diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp index 179d1d4a..b0043379 100644 --- a/misc/dev/dev.cpp +++ b/misc/dev/dev.cpp @@ -154,17 +154,17 @@ enum macpan2_meth // functions that can only take numerical matrices -- no integer vectors std::vector mp_math = { - MP2_ADD, MP2_SUBTRACT, MP2_MULTIPLY, MP2_DIVIDE, MP2_POWER, MP2_EXP, MP2_LOG - , MP2_MATRIX, MP2_MATRIX_MULTIPLY - , MP2_SUM, MP2_ROWSUMS, MP2_COLSUMS, MP2_TRANSPOSE - , MP2_CONVOLUTION, MP2_CBIND, MP2_RBIND, MP2_RECYCLE, MP2_CLAMP - , MP2_POISSON_DENSITY, MP2_NEGBIN_DENSITY, MP2_NORMAL_DENSITY - , MP2_BINOM_DENSITY - , MP2_POISSON_SIM, MP2_NEGBIN_SIM, MP2_NORMAL_SIM, MP2_KRONECKER - , MP2_TO_DIAG, MP2_FROM_DIAG, MP2_COS, MP2_SIN, MP2_COS - , MP2_BINOM_SIM, MP2_EULER_MULTINOM_SIM - , MP2_ROUND, MP2_PGAMMA, MP2_PNORM - , MP2_MEAN, MP2_SD, MP2_INVLOGIT, MP2_LOGIT + MP2_ADD, MP2_SUBTRACT, MP2_MULTIPLY, MP2_DIVIDE, MP2_POWER, MP2_EXP, MP2_LOG + , MP2_MATRIX, MP2_MATRIX_MULTIPLY + , MP2_SUM, MP2_ROWSUMS, MP2_COLSUMS, MP2_TRANSPOSE + , MP2_CONVOLUTION, MP2_CBIND, MP2_RBIND, MP2_RECYCLE, MP2_CLAMP + , MP2_POISSON_DENSITY, MP2_NEGBIN_DENSITY, MP2_NORMAL_DENSITY + , MP2_BINOM_DENSITY + , MP2_POISSON_SIM, MP2_NEGBIN_SIM, MP2_NORMAL_SIM, MP2_KRONECKER + , MP2_TO_DIAG, MP2_FROM_DIAG, MP2_COS, MP2_SIN, MP2_COS + , MP2_BINOM_SIM, MP2_EULER_MULTINOM_SIM + , MP2_ROUND, MP2_PGAMMA, MP2_PNORM + , MP2_MEAN, MP2_SD, MP2_INVLOGIT, MP2_LOGIT }; std::string bail_out_log_file = ".macpan2/bail-out/log.txt"; @@ -218,6 +218,40 @@ Type FUN(const Type x, const int &index) { \ #define MP2_ERR(CODE, MSG, FUN) SetError(CODE, MSG, row, FUN, args.all_rows(), args.all_cols(), args.all_type_ints(), t); \ +#define MP2_REPORT_ERROR_FINAL(evaluator) { \ + int error = evaluator.GetErrorCode(); \ + int expr_row = evaluator.GetExprRow(); \ + vector arg_rows = evaluator.GetArgRows(); \ + vector arg_cols = evaluator.GetArgCols(); \ + vector arg_type_ints = evaluator.GetArgTypeInts();\ + int func_int = evaluator.GetFuncInt(); \ + int time_int = evaluator.GetTimeInt(); \ + const char *err_msg = evaluator.GetErrorMessage(); \ + REPORT(error); \ + REPORT(expr_row); \ + REPORT(func_int); \ + REPORT(time_int); \ + REPORT(arg_rows); \ + REPORT(arg_cols); \ + REPORT(arg_type_ints); \ + \ + logfile.open(log_file, std::ios_base::app); \ + logfile << "Error code = " << error << std::endl; \ + logfile << "Error message = " << err_msg << std::endl; \ + logfile << "Expression row = " << expr_row << std::endl;\ + logfile << "Function code = " << func_int << std::endl;\ + logfile << "Time step = " << time_int << std::endl; \ + logfile.close(); \ + return 0.0; \ +} \ + +#define MP2_REPORT_ERROR(evaluator) \ + { \ + if (evaluator.GetErrorCode()) { \ + MP2_REPORT_ERROR_FINAL(evaluator) \ + } \ + } \ + // ENGINE FUNCTIONS USED BY OTHER ENGINE FUNCTIONS @@ -361,9 +395,8 @@ struct ListOfMatrices // containing the concatenation of the vectors and // another containing the lengths of each concatenated // vector -class ListOfIntVecs -{ -public: +class ListOfIntVecs { +public: // ListOfIntVecs // this nestedVector will contain examples of unflattened // integer vectors std::vector> nestedVector; @@ -569,9 +602,8 @@ matrix getNthMat( } template -class ArgList -{ -public: +class ArgList { +public: // ArgList enum class ItemType { Matrix, IntVector @@ -889,28 +921,20 @@ class ArgList return error_code_; } -private: - struct Item - { +private: // ArgList + struct Item { ItemType type; matrix mat; std::vector intVec; }; - std::vector items_; int size_; int error_code_ = 0; // Initialize the error code to 0 (no error) by default -}; -// private: -// std::vector items_; -// int size_; -// int error_code_ = 0; // Initialize the error code to 0 (no error) by default -// }; +}; // end ArgList template -class ExprEvaluator -{ -private: +class ExprEvaluator { +private: // ExprEvaluator vector mats_save_hist; vector table_x; vector table_n; @@ -921,7 +945,7 @@ class ExprEvaluator ListOfIntVecs valid_int_vecs; vector valid_literals; -public: +public: // ExprEvaluator // constructor ExprEvaluator( vector &mats_save_hist_, @@ -1014,8 +1038,7 @@ class ExprEvaluator // Check if error has already happened at some point // of the recursive call of EvalExpr. - if (GetErrorCode()) - return m; + if (GetErrorCode()) return m; switch (table_n[row]) { case -2: // methods (pre-processed matrices) @@ -1114,14 +1137,12 @@ class ExprEvaluator vector index2mats(n); vector index2what(n); for (int i = 0; i < n; i++) { - if (table_n[table_i[row] + i] == -3) - { + if (table_n[table_i[row] + i] == -3) { // -3 in the 'number of arguments' column of the // parse table means 'integer vector' args.set(i, valid_int_vecs[table_x[table_i[row] + i]]); } - else - { + else { // otherwise, recursively descend into the parse tree // to pick out the arguments args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i)); @@ -1136,23 +1157,19 @@ class ExprEvaluator // or something. // // index2what = 0 (for a matrix), 1 (for an int vec), -1 (for something else) - if (table_n[table_i[row] + i] == 0) - { + if (table_n[table_i[row] + i] == 0) { index2mats[i] = table_x[table_i[row] + i]; index2what[i] = 0; // pointing at matrix } - else if (table_n[table_i[row] + i] == -3) - { + else if (table_n[table_i[row] + i] == -3) { index2mats[i] = table_x[table_i[row] + i]; index2what[i] = 1; // pointing at integer vector } - else - { + else { index2mats[i] = -1; index2what[i] = -1; } - if (GetErrorCode()) - return m; + if (GetErrorCode()) return m; } if (is_int_in(table_x[row] + 1, mp_math)) { @@ -3335,7 +3352,7 @@ class ExprEvaluator } // switch (table_n[row]) }; -private: +private: // ExprEvaluator unsigned char error_code; int expr_row; int func_int; @@ -3344,33 +3361,68 @@ class ExprEvaluator std::vector arg_cols; std::vector arg_type_ints; char error_message[256]; -}; +}; // end ExprEvaluator + + template class MatAssigner { -private: +private: // MatAssigner vector table_x; vector table_n; vector table_i; ListOfIntVecs valid_int_vecs; vector valid_literals; -public: +public: // MatAssigner MatAssigner( vector &table_x_, vector &table_n_, vector &table_i_, ListOfIntVecs &valid_int_vecs_, - vector &valid_literals_) { + vector &valid_literals_ + ) { + error_code = 0; // non-zero means error has occurred; otherwise, no error + expr_row = 0; + func_int = -99; // assume no function information is available + time_int = 0; + arg_rows = {0}; + arg_cols = {0}; + arg_type_ints = {0}; table_x = table_x_; table_n = table_n_; table_i = table_i_; valid_int_vecs = valid_int_vecs_; valid_literals = valid_literals_; + + strcpy(error_message, "OK"); + }; + + // getters + unsigned char GetErrorCode() { return error_code; }; + const char *GetErrorMessage() { return error_message; }; + int GetExprRow() { return expr_row; }; + int GetFuncInt() { return func_int; }; + int GetTimeInt() { return time_int; }; + std::vector GetArgRows() { return arg_rows; }; + std::vector GetArgCols() { return arg_cols; }; + std::vector GetArgTypeInts() { return arg_type_ints; }; + + // setters + void SetError(unsigned char code, const char *message, int e_row, int f_int, std::vector a_rows, std::vector a_cols, std::vector a_type_ints, int t_int) { + error_code = code; + expr_row = e_row; + func_int = f_int; + time_int = t_int; + arg_rows = a_rows; + arg_cols = a_cols; + arg_type_ints = a_type_ints; + strcpy(error_message, message); }; void matAssign( matrix assignment_value, + int t, // current time step ListOfMatrices &valid_vars, // current list of values of each matrix int row = 0 // current expression parse table row being evaluated ) { @@ -3388,9 +3440,13 @@ class MatAssigner { std::vector v1; std::vector v2; int err_code; + matrix for_literals; + for_literals = matrix::Zero(1, 1); // Rcpp::Rcout << "---- assignment ----" << std::endl; // Rcpp::Rcout << "n: " << n << std::endl; // Rcpp::Rcout << "x: " << n << std::endl; + + if (GetErrorCode()) return; switch (n) { case 0: valid_vars.m_matrices[x] = assignment_value; @@ -3408,6 +3464,24 @@ class MatAssigner { // to switch on the particular function being used. at most one // function can be used on the left-hand-side, and only particular ones // can be used as the switch statement shows + // + + ArgList args(n); + for (int i = 0; i < n; i++) { + if (table_n[table_i[row] + i] == -3) { + // -3 in the 'number of arguments' column of the + // parse table means 'integer vector' + args.set(i, valid_int_vecs[table_x[table_i[row] + i]]); + } else if (table_n[table_i[row] + i] == -1) { + // -1 in the 'number of arguments' column of the + // parse table means 'literal' + for_literals.coeffRef(0, 0) = valid_literals[table_x[table_i[row] + i]]; + args.set(i, for_literals); + } else if (table_n[table_i[row] + i] == 0) { + args.set(i, valid_vars.m_matrices[table_x[table_i[row] + i]]); + } + } + // // #' ## Assignment // #' @@ -3442,9 +3516,9 @@ class MatAssigner { if (n == 3) { // two index vectors if (table_n[table_i[row] + 2] == -1) { // second index vector is a literal v2.push_back(CppAD::Integer(valid_literals[table_x[table_i[row] + 2]])); - //Rf_error("indexing on the left-hand-side cannot be done using literals"); } else if (table_n[table_i[row] + 2] != -3) { // second index vector is not an integer vector - Rf_error("indexing on the left-hand-side needs to be done using integer vectors or literals"); + MP2_ERR(MP2_SQUARE_BRACKET, "indexing on the left-hand-side needs to be done using integer vectors or literals", MP2_SQUARE_BRACKET); + return; } } else if (n == 2){ // one index vector v2.push_back(0); // assume the second index vector is length-1 with a 0 (i.e. points to the first column) @@ -3454,7 +3528,8 @@ class MatAssigner { if (table_n[table_i[row] + 1] == -1) { // first index vector is a literal v1.push_back(CppAD::Integer(valid_literals[table_x[table_i[row] + 1]])); } else if (table_n[table_i[row] + 1] != -3) { // first index vector is not an integer vector or literal - Rf_error("indexing on the left-hand-side needs to be done using integer vectors or literals"); + MP2_ERR(MP2_SQUARE_BRACKET, "indexing on the left-hand-side needs to be done using integer vectors or literals", MP2_SQUARE_BRACKET); + return; } else { // first index is an integer vector v1 = valid_int_vecs[table_x[table_i[row] + 1]]; } @@ -3517,34 +3592,17 @@ class MatAssigner { Rf_error("square bracket (e.g. x[i, j]) and concatenation (e.g. c(x, y, z)) are the only functions allowed on the left-hand-side"); }; -}; +private: // MatAssigner + unsigned char error_code; + int expr_row; + int func_int; + int time_int; + std::vector arg_rows; + std::vector arg_cols; + std::vector arg_type_ints; + char error_message[256]; +}; // end MatAssigner -#define REPORT_ERROR \ - { \ - int error = exprEvaluator.GetErrorCode(); \ - int expr_row = exprEvaluator.GetExprRow(); \ - vector arg_rows = exprEvaluator.GetArgRows(); \ - vector arg_cols = exprEvaluator.GetArgCols(); \ - vector arg_type_ints = exprEvaluator.GetArgTypeInts(); \ - int func_int = exprEvaluator.GetFuncInt(); \ - int time_int = exprEvaluator.GetTimeInt(); \ - const char *err_msg = exprEvaluator.GetErrorMessage(); \ - REPORT(error); \ - REPORT(expr_row); \ - REPORT(func_int); \ - REPORT(time_int); \ - REPORT(arg_rows); \ - REPORT(arg_cols); \ - REPORT(arg_type_ints); \ - \ - logfile.open(log_file, std::ios_base::app); \ - logfile << "Error code = " << error << std::endl; \ - logfile << "Error message = " << err_msg << std::endl; \ - logfile << "Expression row = " << expr_row << std::endl; \ - logfile << "Function code = " << func_int << std::endl; \ - logfile << "Time step = " << time_int << std::endl; \ - logfile.close(); \ - } template vector> MakeSimulationHistory( @@ -3556,7 +3614,6 @@ vector> MakeSimulationHistory( vector> simulation_history(time_steps + 2); matrix empty_matrix; for (unsigned int i = 0; i < mats_save_hist.size(); i++) - //Rcpp::Rcout << "matrix: " << size << std::endl; if (mats_save_hist[i] == 0) hist_shape_template.m_matrices[i] = empty_matrix; @@ -3572,8 +3629,6 @@ void UpdateSimulationHistory( const vector &mats_save_hist, ListOfMatrices &hist_shape_template) { - // matrix emptyMat; - // ListOfMatrices ms(mats); // if the history of the matrix is not to be saved, // just save a 1-by-1 with a zero instead to save space for (unsigned int i = 0; i < mats_save_hist.size(); i++) @@ -3741,7 +3796,8 @@ Type objective_function::operator()() meth_mats, meth_int_vecs, const_int_vecs, - literals); + literals + ); ExprEvaluator objFunEvaluator( mats_save_hist, // this seems odd given that objective functions can't access history o_table_x, @@ -3751,13 +3807,15 @@ Type objective_function::operator()() meth_mats, meth_int_vecs, const_int_vecs, - literals); + literals + ); MatAssigner matAssigner( a_table_x, a_table_n, a_table_i, const_int_vecs, - literals); + literals + ); // 3 Pre-simulation (the 'before' step) int expr_index = 0; @@ -3780,13 +3838,9 @@ Type objective_function::operator()() result = exprEvaluator.EvalExpr( simulation_history, 0, mats, p_table_row); - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR - return 0.0; - } + MP2_REPORT_ERROR(exprEvaluator) - // mats.m_matrices[expr_output_id[expr_index+i]] = result; - matAssigner.matAssign(result, mats, a_table_row); + matAssigner.matAssign(result, 0, mats, a_table_row); p_table_row += expr_num_p_table_rows[i]; a_table_row += assign_num_a_table_rows[i]; @@ -3795,7 +3849,6 @@ Type objective_function::operator()() // to the start of that loop } // p_table_row is fine here - // simulation_history[0] = mats; UpdateSimulationHistory( simulation_history, 0, @@ -3825,19 +3878,19 @@ Type objective_function::operator()() if (expr_sim_block[i] == 1) { SIMULATE { result = exprEvaluator.EvalExpr( - simulation_history, k + 1, mats, p_table_row2); + simulation_history, k + 1, mats, p_table_row2 + ); } } else result = exprEvaluator.EvalExpr( - simulation_history, k + 1, mats, p_table_row2); + simulation_history, k + 1, mats, p_table_row2 + ); - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR - return 0.0; - } - // mats.m_matrices[expr_output_id[expr_index+i]] = result; - matAssigner.matAssign(result, mats, a_table_row2); + MP2_REPORT_ERROR(exprEvaluator) + matAssigner.matAssign(result, k + 1, mats, a_table_row2); + Rcpp::Rcout << matAssigner.GetErrorCode() << std::endl; + MP2_REPORT_ERROR(matAssigner) p_table_row2 += expr_num_p_table_rows[expr_index + i]; a_table_row2 += assign_num_a_table_rows[expr_index + i]; @@ -3848,7 +3901,6 @@ Type objective_function::operator()() Rcpp::Rcout << "mats = " << mats.m_matrices[ii] << std::endl; #endif } - // simulation_history[k+1] = mats; UpdateSimulationHistory( simulation_history, k + 1, @@ -3879,19 +3931,13 @@ Type objective_function::operator()() simulation_history, time_steps + 1, mats, p_table_row); } - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR - return 0.0; - } - - // mats.m_matrices[expr_output_id[expr_index+i]] = result; - matAssigner.matAssign(result, mats, a_table_row); + MP2_REPORT_ERROR(exprEvaluator) + matAssigner.matAssign(result, time_steps + 1, mats, a_table_row); p_table_row += expr_num_p_table_rows[expr_index + i]; a_table_row += assign_num_a_table_rows[expr_index + i]; } - // simulation_history[time_steps+1] = mats; UpdateSimulationHistory( simulation_history, time_steps + 1, @@ -3937,7 +3983,10 @@ Type objective_function::operator()() int cur = 0; for (unsigned int i = 0; i < mats_return.size(); i++) { if (mats_return[i] == 1) { - if (mats_save_hist[i] == 0) { // Report the last one + + // Report the last one + // (just 'after' phases) + if (mats_save_hist[i] == 0) { for (int jj = 0; jj < mats.m_matrices[i].cols(); jj++) for (int ii = 0; ii < mats.m_matrices[i].rows(); ii++) { values(cur, 0) = i; @@ -3948,7 +3997,10 @@ Type objective_function::operator()() cur++; } } - else { // Report the whole simulation history + + // Report the whole simulation history + // (including 'before', 'during' and 'after', phasese) + else { int hist_len = time_steps + 2; for (int k = 0; k < hist_len; k++) for (int jj = 0; jj < simulation_history[k].m_matrices[i].cols(); jj++) @@ -3973,14 +4025,12 @@ Type objective_function::operator()() // 7 Calc the return of the objective function matrix ret; ret = objFunEvaluator.EvalExpr(simulation_history, time_steps + 2, mats, 0); - if (ret.size() != 1) Rf_error("Objective function did not return a scalar."); - - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR; - return 0.0; - } - REPORT_ERROR + MP2_REPORT_ERROR(exprEvaluator); + MP2_REPORT_ERROR(objFunEvaluator); + if (ret.size() != 1) Rf_error("Objective function did not return a scalar."); + + MP2_REPORT_ERROR_FINAL(exprEvaluator); #ifdef MP_VERBOSE Rcpp::Rcout << "======== end of objective function ========" << std::endl; diff --git a/misc/notes b/misc/notes index e5cb557d..63d397a0 100644 --- a/misc/notes +++ b/misc/notes @@ -1,2 +1,33 @@ Julia Gog, Viggo Andreasen -https://www.overleaf.com/login +https://link.springer.com/article/10.1007/s00285-015-0873-4 + +Looking more at odin-monty. + +This design philosophy is different from the one I have in mind for macpan2 : https://mrc-ide.github.io/odin-monty/slides/fitting.html#/design-philosophy + + + +This roadmap is focused on things I don’t understand and things I do understand that are about computational efficiency https://mrc-ide.github.io/odin-monty/slides/fitting.html#/roadmap + +* Update on the Plague project +* Dev work + * Quiet pkg check + * Document adding of engine functions in CONTRIBUTING.md, with some simplification + * New engine functions: `invlogit`, `logit`, `sin`, `sqrt` + * pkgdown (nicer formatting on pkgdown topic titles, Fix rendering of math equations) + * Complete engine function documentation + * Using `mp_state_vars`/`mp_flow_vars` more consistently + * Documented how assignment works (i.e., functions on the left hand side) +* GUI fun +* Differences between us and odin2: + * state space model and MCMC is the generic approach + * design philosophy is provide modular pieces with 'minimal handholding' +* Summer workshop +* Roadmap + * Vector valued FOI + * Log linear model of time variation + * Distribution specification documentation and hierarchical and binomial + * + + +https://github.com/canmod/macpan2/compare/7ef2c8053fa475d00b14cd2c87b2392d7f5df297..main \ No newline at end of file diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 54ae64d5..2ebf3a2e 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -155,17 +155,17 @@ enum macpan2_meth // functions that can only take numerical matrices -- no integer vectors std::vector mp_math = { - MP2_ADD, MP2_SUBTRACT, MP2_MULTIPLY, MP2_DIVIDE, MP2_POWER, MP2_EXP, MP2_LOG - , MP2_MATRIX, MP2_MATRIX_MULTIPLY - , MP2_SUM, MP2_ROWSUMS, MP2_COLSUMS, MP2_TRANSPOSE - , MP2_CONVOLUTION, MP2_CBIND, MP2_RBIND, MP2_RECYCLE, MP2_CLAMP - , MP2_POISSON_DENSITY, MP2_NEGBIN_DENSITY, MP2_NORMAL_DENSITY - , MP2_BINOM_DENSITY - , MP2_POISSON_SIM, MP2_NEGBIN_SIM, MP2_NORMAL_SIM, MP2_KRONECKER - , MP2_TO_DIAG, MP2_FROM_DIAG, MP2_COS, MP2_SIN, MP2_COS - , MP2_BINOM_SIM, MP2_EULER_MULTINOM_SIM - , MP2_ROUND, MP2_PGAMMA, MP2_PNORM - , MP2_MEAN, MP2_SD, MP2_INVLOGIT, MP2_LOGIT + MP2_ADD, MP2_SUBTRACT, MP2_MULTIPLY, MP2_DIVIDE, MP2_POWER, MP2_EXP, MP2_LOG + , MP2_MATRIX, MP2_MATRIX_MULTIPLY + , MP2_SUM, MP2_ROWSUMS, MP2_COLSUMS, MP2_TRANSPOSE + , MP2_CONVOLUTION, MP2_CBIND, MP2_RBIND, MP2_RECYCLE, MP2_CLAMP + , MP2_POISSON_DENSITY, MP2_NEGBIN_DENSITY, MP2_NORMAL_DENSITY + , MP2_BINOM_DENSITY + , MP2_POISSON_SIM, MP2_NEGBIN_SIM, MP2_NORMAL_SIM, MP2_KRONECKER + , MP2_TO_DIAG, MP2_FROM_DIAG, MP2_COS, MP2_SIN, MP2_COS + , MP2_BINOM_SIM, MP2_EULER_MULTINOM_SIM + , MP2_ROUND, MP2_PGAMMA, MP2_PNORM + , MP2_MEAN, MP2_SD, MP2_INVLOGIT, MP2_LOGIT }; std::string bail_out_log_file = ".macpan2/bail-out/log.txt"; @@ -219,6 +219,40 @@ Type FUN(const Type x, const int &index) { \ #define MP2_ERR(CODE, MSG, FUN) SetError(CODE, MSG, row, FUN, args.all_rows(), args.all_cols(), args.all_type_ints(), t); \ +#define MP2_REPORT_ERROR_FINAL(evaluator) { \ + int error = evaluator.GetErrorCode(); \ + int expr_row = evaluator.GetExprRow(); \ + vector arg_rows = evaluator.GetArgRows(); \ + vector arg_cols = evaluator.GetArgCols(); \ + vector arg_type_ints = evaluator.GetArgTypeInts(); \ + int func_int = evaluator.GetFuncInt(); \ + int time_int = evaluator.GetTimeInt(); \ + const char *err_msg = evaluator.GetErrorMessage(); \ + REPORT(error); \ + REPORT(expr_row); \ + REPORT(func_int); \ + REPORT(time_int); \ + REPORT(arg_rows); \ + REPORT(arg_cols); \ + REPORT(arg_type_ints); \ + \ + logfile.open(log_file, std::ios_base::app); \ + logfile << "Error code = " << error << std::endl; \ + logfile << "Error message = " << err_msg << std::endl; \ + logfile << "Expression row = " << expr_row << std::endl; \ + logfile << "Function code = " << func_int << std::endl; \ + logfile << "Time step = " << time_int << std::endl; \ + logfile.close(); \ + return 0.0; \ +} \ + +#define MP2_REPORT_ERROR(evaluator) \ + { \ + if (evaluator.GetErrorCode()) { \ + MP2_REPORT_ERROR_FINAL(evaluator) \ + } \ + } \ + // ENGINE FUNCTIONS USED BY OTHER ENGINE FUNCTIONS @@ -362,9 +396,8 @@ struct ListOfMatrices // containing the concatenation of the vectors and // another containing the lengths of each concatenated // vector -class ListOfIntVecs -{ -public: +class ListOfIntVecs { +public: // ListOfIntVecs // this nestedVector will contain examples of unflattened // integer vectors std::vector> nestedVector; @@ -570,9 +603,8 @@ matrix getNthMat( } template -class ArgList -{ -public: +class ArgList { +public: // ArgList enum class ItemType { Matrix, IntVector @@ -890,28 +922,20 @@ class ArgList return error_code_; } -private: - struct Item - { +private: // ArgList + struct Item { ItemType type; matrix mat; std::vector intVec; }; - std::vector items_; int size_; int error_code_ = 0; // Initialize the error code to 0 (no error) by default -}; -// private: -// std::vector items_; -// int size_; -// int error_code_ = 0; // Initialize the error code to 0 (no error) by default -// }; +}; // end ArgList template -class ExprEvaluator -{ -private: +class ExprEvaluator { +private: // ExprEvaluator vector mats_save_hist; vector table_x; vector table_n; @@ -922,7 +946,7 @@ class ExprEvaluator ListOfIntVecs valid_int_vecs; vector valid_literals; -public: +public: // ExprEvaluator // constructor ExprEvaluator( vector &mats_save_hist_, @@ -1015,8 +1039,7 @@ class ExprEvaluator // Check if error has already happened at some point // of the recursive call of EvalExpr. - if (GetErrorCode()) - return m; + if (GetErrorCode()) return m; switch (table_n[row]) { case -2: // methods (pre-processed matrices) @@ -1115,14 +1138,12 @@ class ExprEvaluator vector index2mats(n); vector index2what(n); for (int i = 0; i < n; i++) { - if (table_n[table_i[row] + i] == -3) - { + if (table_n[table_i[row] + i] == -3) { // -3 in the 'number of arguments' column of the // parse table means 'integer vector' args.set(i, valid_int_vecs[table_x[table_i[row] + i]]); } - else - { + else { // otherwise, recursively descend into the parse tree // to pick out the arguments args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i)); @@ -1137,23 +1158,19 @@ class ExprEvaluator // or something. // // index2what = 0 (for a matrix), 1 (for an int vec), -1 (for something else) - if (table_n[table_i[row] + i] == 0) - { + if (table_n[table_i[row] + i] == 0) { index2mats[i] = table_x[table_i[row] + i]; index2what[i] = 0; // pointing at matrix } - else if (table_n[table_i[row] + i] == -3) - { + else if (table_n[table_i[row] + i] == -3) { index2mats[i] = table_x[table_i[row] + i]; index2what[i] = 1; // pointing at integer vector } - else - { + else { index2mats[i] = -1; index2what[i] = -1; } - if (GetErrorCode()) - return m; + if (GetErrorCode()) return m; } if (is_int_in(table_x[row] + 1, mp_math)) { @@ -3336,7 +3353,7 @@ class ExprEvaluator } // switch (table_n[row]) }; -private: +private: // ExprEvaluator unsigned char error_code; int expr_row; int func_int; @@ -3345,33 +3362,68 @@ class ExprEvaluator std::vector arg_cols; std::vector arg_type_ints; char error_message[256]; -}; +}; // end ExprEvaluator + + template class MatAssigner { -private: +private: // MatAssigner vector table_x; vector table_n; vector table_i; ListOfIntVecs valid_int_vecs; vector valid_literals; -public: +public: // MatAssigner MatAssigner( vector &table_x_, vector &table_n_, vector &table_i_, ListOfIntVecs &valid_int_vecs_, - vector &valid_literals_) { + vector &valid_literals_ + ) { + error_code = 0; // non-zero means error has occurred; otherwise, no error + expr_row = 0; + func_int = -99; // assume no function information is available + time_int = 0; + arg_rows = {0}; + arg_cols = {0}; + arg_type_ints = {0}; table_x = table_x_; table_n = table_n_; table_i = table_i_; valid_int_vecs = valid_int_vecs_; valid_literals = valid_literals_; + + strcpy(error_message, "OK"); + }; + + // getters + unsigned char GetErrorCode() { return error_code; }; + const char *GetErrorMessage() { return error_message; }; + int GetExprRow() { return expr_row; }; + int GetFuncInt() { return func_int; }; + int GetTimeInt() { return time_int; }; + std::vector GetArgRows() { return arg_rows; }; + std::vector GetArgCols() { return arg_cols; }; + std::vector GetArgTypeInts() { return arg_type_ints; }; + + // setters + void SetError(unsigned char code, const char *message, int e_row, int f_int, std::vector a_rows, std::vector a_cols, std::vector a_type_ints, int t_int) { + error_code = code; + expr_row = e_row; + func_int = f_int; + time_int = t_int; + arg_rows = a_rows; + arg_cols = a_cols; + arg_type_ints = a_type_ints; + strcpy(error_message, message); }; void matAssign( matrix assignment_value, + int t, // current time step ListOfMatrices &valid_vars, // current list of values of each matrix int row = 0 // current expression parse table row being evaluated ) { @@ -3389,9 +3441,13 @@ class MatAssigner { std::vector v1; std::vector v2; int err_code; + matrix for_literals; + for_literals = matrix::Zero(1, 1); // Rcpp::Rcout << "---- assignment ----" << std::endl; // Rcpp::Rcout << "n: " << n << std::endl; // Rcpp::Rcout << "x: " << n << std::endl; + + if (GetErrorCode()) return; switch (n) { case 0: valid_vars.m_matrices[x] = assignment_value; @@ -3409,6 +3465,24 @@ class MatAssigner { // to switch on the particular function being used. at most one // function can be used on the left-hand-side, and only particular ones // can be used as the switch statement shows + // + + ArgList args(n); + for (int i = 0; i < n; i++) { + if (table_n[table_i[row] + i] == -3) { + // -3 in the 'number of arguments' column of the + // parse table means 'integer vector' + args.set(i, valid_int_vecs[table_x[table_i[row] + i]]); + } else if (table_n[table_i[row] + i] == -1) { + // -1 in the 'number of arguments' column of the + // parse table means 'literal' + for_literals.coeffRef(0, 0) = valid_literals[table_x[table_i[row] + i]]; + args.set(i, for_literals); + } else if (table_n[table_i[row] + i] == 0) { + args.set(i, valid_vars.m_matrices[table_x[table_i[row] + i]]); + } + } + // // #' ## Assignment // #' @@ -3443,9 +3517,9 @@ class MatAssigner { if (n == 3) { // two index vectors if (table_n[table_i[row] + 2] == -1) { // second index vector is a literal v2.push_back(CppAD::Integer(valid_literals[table_x[table_i[row] + 2]])); - //Rf_error("indexing on the left-hand-side cannot be done using literals"); } else if (table_n[table_i[row] + 2] != -3) { // second index vector is not an integer vector - Rf_error("indexing on the left-hand-side needs to be done using integer vectors or literals"); + MP2_ERR(MP2_SQUARE_BRACKET, "indexing on the left-hand-side needs to be done using integer vectors or literals", MP2_SQUARE_BRACKET); + return; } } else if (n == 2){ // one index vector v2.push_back(0); // assume the second index vector is length-1 with a 0 (i.e. points to the first column) @@ -3455,7 +3529,8 @@ class MatAssigner { if (table_n[table_i[row] + 1] == -1) { // first index vector is a literal v1.push_back(CppAD::Integer(valid_literals[table_x[table_i[row] + 1]])); } else if (table_n[table_i[row] + 1] != -3) { // first index vector is not an integer vector or literal - Rf_error("indexing on the left-hand-side needs to be done using integer vectors or literals"); + MP2_ERR(MP2_SQUARE_BRACKET, "indexing on the left-hand-side needs to be done using integer vectors or literals", MP2_SQUARE_BRACKET); + return; } else { // first index is an integer vector v1 = valid_int_vecs[table_x[table_i[row] + 1]]; } @@ -3518,34 +3593,17 @@ class MatAssigner { Rf_error("square bracket (e.g. x[i, j]) and concatenation (e.g. c(x, y, z)) are the only functions allowed on the left-hand-side"); }; -}; +private: // MatAssigner + unsigned char error_code; + int expr_row; + int func_int; + int time_int; + std::vector arg_rows; + std::vector arg_cols; + std::vector arg_type_ints; + char error_message[256]; +}; // end MatAssigner -#define REPORT_ERROR \ - { \ - int error = exprEvaluator.GetErrorCode(); \ - int expr_row = exprEvaluator.GetExprRow(); \ - vector arg_rows = exprEvaluator.GetArgRows(); \ - vector arg_cols = exprEvaluator.GetArgCols(); \ - vector arg_type_ints = exprEvaluator.GetArgTypeInts(); \ - int func_int = exprEvaluator.GetFuncInt(); \ - int time_int = exprEvaluator.GetTimeInt(); \ - const char *err_msg = exprEvaluator.GetErrorMessage(); \ - REPORT(error); \ - REPORT(expr_row); \ - REPORT(func_int); \ - REPORT(time_int); \ - REPORT(arg_rows); \ - REPORT(arg_cols); \ - REPORT(arg_type_ints); \ - \ - logfile.open(log_file, std::ios_base::app); \ - logfile << "Error code = " << error << std::endl; \ - logfile << "Error message = " << err_msg << std::endl; \ - logfile << "Expression row = " << expr_row << std::endl; \ - logfile << "Function code = " << func_int << std::endl; \ - logfile << "Time step = " << time_int << std::endl; \ - logfile.close(); \ - } template vector> MakeSimulationHistory( @@ -3557,7 +3615,6 @@ vector> MakeSimulationHistory( vector> simulation_history(time_steps + 2); matrix empty_matrix; for (unsigned int i = 0; i < mats_save_hist.size(); i++) - //Rcpp::Rcout << "matrix: " << size << std::endl; if (mats_save_hist[i] == 0) hist_shape_template.m_matrices[i] = empty_matrix; @@ -3573,8 +3630,6 @@ void UpdateSimulationHistory( const vector &mats_save_hist, ListOfMatrices &hist_shape_template) { - // matrix emptyMat; - // ListOfMatrices ms(mats); // if the history of the matrix is not to be saved, // just save a 1-by-1 with a zero instead to save space for (unsigned int i = 0; i < mats_save_hist.size(); i++) @@ -3742,7 +3797,8 @@ Type objective_function::operator()() meth_mats, meth_int_vecs, const_int_vecs, - literals); + literals + ); ExprEvaluator objFunEvaluator( mats_save_hist, // this seems odd given that objective functions can't access history o_table_x, @@ -3752,13 +3808,15 @@ Type objective_function::operator()() meth_mats, meth_int_vecs, const_int_vecs, - literals); + literals + ); MatAssigner matAssigner( a_table_x, a_table_n, a_table_i, const_int_vecs, - literals); + literals + ); // 3 Pre-simulation (the 'before' step) int expr_index = 0; @@ -3781,13 +3839,9 @@ Type objective_function::operator()() result = exprEvaluator.EvalExpr( simulation_history, 0, mats, p_table_row); - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR - return 0.0; - } + MP2_REPORT_ERROR(exprEvaluator) - // mats.m_matrices[expr_output_id[expr_index+i]] = result; - matAssigner.matAssign(result, mats, a_table_row); + matAssigner.matAssign(result, 0, mats, a_table_row); p_table_row += expr_num_p_table_rows[i]; a_table_row += assign_num_a_table_rows[i]; @@ -3796,7 +3850,6 @@ Type objective_function::operator()() // to the start of that loop } // p_table_row is fine here - // simulation_history[0] = mats; UpdateSimulationHistory( simulation_history, 0, @@ -3826,19 +3879,19 @@ Type objective_function::operator()() if (expr_sim_block[i] == 1) { SIMULATE { result = exprEvaluator.EvalExpr( - simulation_history, k + 1, mats, p_table_row2); + simulation_history, k + 1, mats, p_table_row2 + ); } } else result = exprEvaluator.EvalExpr( - simulation_history, k + 1, mats, p_table_row2); + simulation_history, k + 1, mats, p_table_row2 + ); - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR - return 0.0; - } - // mats.m_matrices[expr_output_id[expr_index+i]] = result; - matAssigner.matAssign(result, mats, a_table_row2); + MP2_REPORT_ERROR(exprEvaluator) + matAssigner.matAssign(result, k + 1, mats, a_table_row2); + Rcpp::Rcout << matAssigner.GetErrorCode() << std::endl; + MP2_REPORT_ERROR(matAssigner) p_table_row2 += expr_num_p_table_rows[expr_index + i]; a_table_row2 += assign_num_a_table_rows[expr_index + i]; @@ -3849,7 +3902,6 @@ Type objective_function::operator()() Rcpp::Rcout << "mats = " << mats.m_matrices[ii] << std::endl; #endif } - // simulation_history[k+1] = mats; UpdateSimulationHistory( simulation_history, k + 1, @@ -3880,19 +3932,13 @@ Type objective_function::operator()() simulation_history, time_steps + 1, mats, p_table_row); } - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR - return 0.0; - } - - // mats.m_matrices[expr_output_id[expr_index+i]] = result; - matAssigner.matAssign(result, mats, a_table_row); + MP2_REPORT_ERROR(exprEvaluator) + matAssigner.matAssign(result, time_steps + 1, mats, a_table_row); p_table_row += expr_num_p_table_rows[expr_index + i]; a_table_row += assign_num_a_table_rows[expr_index + i]; } - // simulation_history[time_steps+1] = mats; UpdateSimulationHistory( simulation_history, time_steps + 1, @@ -3938,7 +3984,10 @@ Type objective_function::operator()() int cur = 0; for (unsigned int i = 0; i < mats_return.size(); i++) { if (mats_return[i] == 1) { - if (mats_save_hist[i] == 0) { // Report the last one + + // Report the last one + // (just 'after' phases) + if (mats_save_hist[i] == 0) { for (int jj = 0; jj < mats.m_matrices[i].cols(); jj++) for (int ii = 0; ii < mats.m_matrices[i].rows(); ii++) { values(cur, 0) = i; @@ -3949,7 +3998,10 @@ Type objective_function::operator()() cur++; } } - else { // Report the whole simulation history + + // Report the whole simulation history + // (including 'before', 'during' and 'after', phasese) + else { int hist_len = time_steps + 2; for (int k = 0; k < hist_len; k++) for (int jj = 0; jj < simulation_history[k].m_matrices[i].cols(); jj++) @@ -3974,14 +4026,12 @@ Type objective_function::operator()() // 7 Calc the return of the objective function matrix ret; ret = objFunEvaluator.EvalExpr(simulation_history, time_steps + 2, mats, 0); - if (ret.size() != 1) Rf_error("Objective function did not return a scalar."); - if (exprEvaluator.GetErrorCode()) { - REPORT_ERROR; - return 0.0; - } - - REPORT_ERROR + MP2_REPORT_ERROR(exprEvaluator); + MP2_REPORT_ERROR(objFunEvaluator); + if (ret.size() != 1) Rf_error("Objective function did not return a scalar."); + + MP2_REPORT_ERROR_FINAL(exprEvaluator); #ifdef MP_VERBOSE Rcpp::Rcout << "======== end of objective function ========" << std::endl; diff --git a/tests/testthat/test-convolution.R b/tests/testthat/test-convolution.R index 9288a024..87ee756b 100644 --- a/tests/testthat/test-convolution.R +++ b/tests/testthat/test-convolution.R @@ -36,14 +36,6 @@ X_sim = (sim_data manual_convolution = k[1] * X_sim$value + k[2] * X_sim$lag_1 + k[3] * X_sim$lag_2 macpan2_convolution = sim_data %>% filter(matrix == "Y") %>% pull(value) -# TODO: add this formula to engine docs -# - `l` : index for the time lag, starting with `l = 0` -# - `t` : index for time, starting at `t = 1` -# - `k_l` : `l`th value of the kernel -# - `x_t` : value of the variable at time `t` -# k_0 * x_t + k_1 * x_{t-1} + k_2 * x_{t-2} -# \sum_{l = 0}^{min(m,t)-1} k_l x_{t-l} - ## view differences cbind(manual_convolution, macpan2_convolution) plot(manual_convolution, macpan2_convolution)