From 752c5aff30dd54ea9c7138bb03c00cf21e088fc3 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Thu, 26 Jun 2025 14:30:01 -0400 Subject: [PATCH 1/2] initial early stop work --- R/engine_functions.R | 23 +++++ R/enum.R | 1 + R/mp_tmb_model_spec.R | 1 + inst/starter_models/nfds/README.md | 4 +- man/engine_functions.Rd | 29 ++++++ misc/dev/dev.cpp | 160 ++++++++++++++++++++++------- src/macpan2.cpp | 160 ++++++++++++++++++++++------- 7 files changed, 304 insertions(+), 74 deletions(-) diff --git a/R/engine_functions.R b/R/engine_functions.R index d694340f..6f9bf812 100644 --- a/R/engine_functions.R +++ b/R/engine_functions.R @@ -859,6 +859,28 @@ #' engine_eval(~ check_finite(1/0)) ## returns nothing and throws an error #' ``` #' +#' ## Control Flow +#' +#' Functions for controlling the flow of simulations. +#' +#' ### Functions +#' +#' * `stop_if_lt(x, y)` : Stop the simulation loop and return +#' simulations if any element of `x` is less than the +#' corresponding element in `y`. The elements of `y` are +#' recycled to match those in `x` if necessary and if possible. +#' +#' ### Arguments +#' +#' * `x` : Matrix to compare with `y`. +#' * `y` : Matrix to compare with `x`, with a default value of +#' `0`. +#' +#' ### Return +#' +#' An empty matrix. This function is returned for its effect of +#' breaking out of the simulation loop. +#' #' ## Assign (deprecated) #' #' Assign values to a subset of the elements in a matrix. @@ -1041,6 +1063,7 @@ #' @aliases pnorm #' @aliases invlogit #' @aliases logit +#' @aliases stop_if_lt #' @aliases assign #' @aliases unpack NULL diff --git a/R/enum.R b/R/enum.R index 6d6ff791..15223d31 100644 --- a/R/enum.R +++ b/R/enum.R @@ -59,6 +59,7 @@ valid_func_sigs = c( , "fwrap: pnorm(q, mean, sd)" , "fwrap: invlogit(x)" , "fwrap: logit(x)" + , "fwrap: stop_if_lt(x, y)" , "fwrap: assign(x, i, j, v)" , "fwrap: unpack(x, ...)" ) diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R index 66960efc..82cdd48d 100644 --- a/R/mp_tmb_model_spec.R +++ b/R/mp_tmb_model_spec.R @@ -42,6 +42,7 @@ TMBModelSpec = function( self$update_method$before() , self$update_method$during() , self$update_method$after() + , .simulate_exprs = self$sim_exprs ) } diff --git a/inst/starter_models/nfds/README.md b/inst/starter_models/nfds/README.md index 3579a333..cfdd5601 100644 --- a/inst/starter_models/nfds/README.md +++ b/inst/starter_models/nfds/README.md @@ -355,14 +355,14 @@ start_time <- Sys.time() old_sim = mp_simulator(old_spec,time_steps = 1L,outputs = "Y") end_time <- Sys.time() end_time - start_time -#> Time difference of 0.1707058 secs +#> Time difference of 0.225244 secs # new spec start_time <- Sys.time() new_sim = mp_simulator(new_spec,time_steps = 1L,outputs = "Y") end_time <- Sys.time() end_time - start_time -#> Time difference of 0.1190519 secs +#> Time difference of 0.187005 secs ``` Creating a simulator with one time step results in the new model diff --git a/man/engine_functions.Rd b/man/engine_functions.Rd index 89d9c79d..1045e51d 100644 --- a/man/engine_functions.Rd +++ b/man/engine_functions.Rd @@ -61,6 +61,7 @@ \alias{pnorm} \alias{invlogit} \alias{logit} +\alias{stop_if_lt} \alias{assign} \alias{unpack} \title{Functions Available in the Simulation Engine} @@ -1047,6 +1048,34 @@ engine_eval(~ check_finite(1/0)) ## returns nothing and throws an error } +\subsection{Control Flow}{ + +Functions for controlling the flow of simulations. +\subsection{Functions}{ +\itemize{ +\item \code{stop_if_lt(x, y)} : Stop the simulation loop and return +simulations if any element of \code{x} is less than the +corresponding element in \code{y}. The elements of \code{y} are +recycled to match those in \code{x} if necessary and if possible. +} +} + +\subsection{Arguments}{ +\itemize{ +\item \code{x} : Matrix to compare with \code{y}. +\item \code{y} : Matrix to compare with \code{x}, with a default value of +\code{0}. +} +} + +\subsection{Return}{ + +An empty matrix. This function is returned for its effect of +breaking out of the simulation loop. +} + +} + \subsection{Assign (deprecated)}{ Assign values to a subset of the elements in a matrix. diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp index 5bf95de7..e362799f 100644 --- a/misc/dev/dev.cpp +++ b/misc/dev/dev.cpp @@ -135,8 +135,9 @@ enum macpan2_func , MP2_PNORM = 57 // fwrap: pnorm(q, mean, sd) , MP2_INVLOGIT = 58 // fwrap: invlogit(x) , MP2_LOGIT = 59 // fwrap: logit(x) - , MP2_ASSIGN = 60 // fwrap: assign(x, i, j, v) - , MP2_UNPACK = 61 // fwrap: unpack(x, ...) + , MP2_STOP_IF_LT = 60 // fwrap: stop_if_lt(x, y) + , MP2_ASSIGN = 61 // fwrap: assign(x, i, j, v) + , MP2_UNPACK = 62 // fwrap: unpack(x, ...) }; enum macpan2_meth @@ -824,7 +825,6 @@ class ArgList // Rcpp::Rcout << "step f" << std::endl; // Rcpp::Rcout << "bad column vector" << std::endl; error_code = 501; - // break; // Exit the loop on error } } else if (mat.cols() == cols) { @@ -838,7 +838,6 @@ class ArgList { // Rcpp::Rcout << "bad row vector" << std::endl; error_code = 501; - // break; // Exit the loop on error } } else @@ -846,7 +845,6 @@ class ArgList // Rcpp::Rcout << "really bad" << std::endl; // Rcpp::Rcout << "step g" << std::endl; error_code = 501; - // break; // Exit the loop on error } if (error_code != 0) @@ -900,11 +898,6 @@ class ArgList 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 -// }; template class ExprEvaluator @@ -919,6 +912,7 @@ class ExprEvaluator ListOfIntVecs meth_int_vecs; ListOfIntVecs valid_int_vecs; vector valid_literals; + bool early_stop; public: // constructor @@ -950,6 +944,10 @@ class ExprEvaluator meth_int_vecs = meth_int_vecs_; valid_int_vecs = valid_int_vecs_; valid_literals = valid_literals_; + // Rcpp::Rcout << "setting early_stop to false" << std::endl; + early_stop = false; + // Rcpp::Rcout << "early_stop now equals: " << early_stop << std::endl; + // Rcpp::Rcout << "but GetEarlyStop() now returns: " << GetEarlyStop() << std::endl; strcpy(error_message, "OK"); }; @@ -963,6 +961,7 @@ class ExprEvaluator std::vector GetArgRows() { return arg_rows; }; std::vector GetArgCols() { return arg_cols; }; std::vector GetArgTypeInts() { return arg_type_ints; }; + bool GetEarlyStop() { return early_stop; }; // 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) { @@ -975,13 +974,17 @@ class ExprEvaluator arg_type_ints = a_type_ints; strcpy(error_message, message); }; + void SetEarlyStop() { + early_stop = true; + }; // evaluator matrix EvalExpr( const vector> &hist, // current simulation history 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 + int row = 0, // current expression parse table row being evaluated + bool simulate_block = false // evaluated within simulate block? ) { // total number of time steps in the simulation loop @@ -1123,7 +1126,7 @@ class ExprEvaluator { // otherwise, recursively descend into the parse tree // to pick out the arguments - args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i)); + args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i, simulate_block)); } // Check here if index2mats actually points at @@ -3090,6 +3093,55 @@ class ExprEvaluator } return m; + + // #' ## Control Flow + // #' + // #' Functions for controlling the flow of simulations. + // #' + // #' ### Functions + // #' + // #' * `stop_if_lt(x, y)` : Stop the simulation loop and return + // #' simulations if any element of `x` is less than the + // #' corresponding element in `y`. The elements of `y` are + // #' recycled to match those in `x` if necessary and if possible. + // #' + // #' ### Arguments + // #' + // #' * `x` : Matrix to compare with `y`. + // #' * `y` : Matrix to compare with `x`, with a default value of + // #' `0`. + // #' + // #' ### Return + // #' + // #' An empty matrix. This function is returned for its effect of + // #' breaking out of the simulation loop. + // #' + case MP2_STOP_IF_LT: + m = args.get_as_mat(0); + rows = m.rows(); + cols = m.cols(); + if (n == 1) { + m1 = matrix::Zero(rows, cols); + } else { + v1.push_back(1); + args = args.recycle_to_shape(v1, rows, cols); + m1 = args.get_as_mat(1); + } + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + if (m.coeff(i, j) < m1.coeff(i, j)) { + //Rcpp::Rcout << "TIME: " << t << std::endl; + //is_finite_mat = GetEarlyStop(); + //Rcpp::Rcout << "Early Stop BEFORE SET: " << is_finite_mat << std::endl; + SetEarlyStop(); + //is_finite_mat = GetEarlyStop(); + //Rcpp::Rcout << "Early Stop AFTER SET: " << is_finite_mat << std::endl; + return m2; // empty matrix + } + } + } + return m2; // empty matrix + // #' ## Assign (deprecated) // #' // #' Assign values to a subset of the elements in a matrix. @@ -3196,7 +3248,7 @@ class ExprEvaluator valid_vars.m_matrices[matIndex].coeffRef(rowIndex, colIndex) = args[3].coeff(k, 0); } return m2; // empty matrix - + // #' ## Unpack (deprecated) // #' // #' Unpack elements of a matrix into smaller matrices. @@ -3475,6 +3527,7 @@ class MatAssigner { int func_int = exprEvaluator.GetFuncInt(); \ int time_int = exprEvaluator.GetTimeInt(); \ const char *err_msg = exprEvaluator.GetErrorMessage(); \ + bool early_stop = exprEvaluator.GetEarlyStop(); \ REPORT(error); \ REPORT(expr_row); \ REPORT(func_int); \ @@ -3489,6 +3542,7 @@ class MatAssigner { logfile << "Expression row = " << expr_row << std::endl; \ logfile << "Function code = " << func_int << std::endl; \ logfile << "Time step = " << time_int << std::endl; \ + logfile << "Early stop = " << early_stop << std::endl; \ logfile.close(); \ } @@ -3665,6 +3719,14 @@ Type objective_function::operator()() for (int i = 0; i < n; i++) mats.m_matrices[r_mat_id[i]].coeffRef(r_row_id[i], r_col_id[i]) = random[r_par_id[i]]; + // in case you need to stop early you can use this to fill the remaining + // iterations of the simulation history + ListOfMatrices empty_mats_list(mats); + matrix empty_matrix; + for (int i = 0; i < empty_mats_list.m_matrices.size(); i++) { + empty_mats_list.m_matrices[i] = empty_matrix; + } + // Simulation history /// each element of this history 'vector' is a list of the matrices /// in the model at a particular point in history @@ -3719,7 +3781,7 @@ Type objective_function::operator()() if (expr_sim_block[i] == 1) { SIMULATE { result = exprEvaluator.EvalExpr( - simulation_history, 0, mats, p_table_row); + simulation_history, 0, mats, p_table_row, true); } } else @@ -3767,23 +3829,35 @@ Type objective_function::operator()() Rcpp::Rcout << "Eval expression --- " << i << std::endl; Rcpp::Rcout << "expr_num_p_table_rows[i] " << expr_num_p_table_rows[expr_index + i] << std::endl; #endif - matrix result; - if (expr_sim_block[i] == 1) { - SIMULATE { + //Rcpp::Rcout << "GetEarlyStop(): " << exprEvaluator.GetEarlyStop() << std::endl; + if (!exprEvaluator.GetEarlyStop()) { + matrix result; + if (expr_sim_block[i] == 1) { + SIMULATE { + result = exprEvaluator.EvalExpr( + simulation_history + , k + 1 + , mats + , p_table_row2 + , true + ); + } + } + else result = exprEvaluator.EvalExpr( - simulation_history, k + 1, mats, p_table_row2); + simulation_history + , k + 1 + , mats + , p_table_row2 + , false + ); + + if (exprEvaluator.GetErrorCode()) { + REPORT_ERROR + return 0.0; } + matAssigner.matAssign(result, mats, a_table_row2); } - else - result = exprEvaluator.EvalExpr( - 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); p_table_row2 += expr_num_p_table_rows[expr_index + i]; a_table_row2 += assign_num_a_table_rows[expr_index + i]; @@ -3794,13 +3868,24 @@ Type objective_function::operator()() Rcpp::Rcout << "mats = " << mats.m_matrices[ii] << std::endl; #endif } - // simulation_history[k+1] = mats; - UpdateSimulationHistory( - simulation_history, - k + 1, - mats, - mats_save_hist, - hist_shape_template); + if (exprEvaluator.GetEarlyStop()) { + UpdateSimulationHistory( + simulation_history, + k + 1, + empty_mats_list, + mats_save_hist, + hist_shape_template + ); + } else { + // simulation_history[k+1] = mats; + UpdateSimulationHistory( + simulation_history, + k + 1, + mats, + mats_save_hist, + hist_shape_template + ); + } } p_table_row = p_table_row2; a_table_row = a_table_row2; @@ -3817,7 +3902,7 @@ Type objective_function::operator()() if (expr_sim_block[i] == 1) { SIMULATE { result = exprEvaluator.EvalExpr( - simulation_history, time_steps + 1, mats, p_table_row); + simulation_history, time_steps + 1, mats, p_table_row, true); } } else { @@ -3844,6 +3929,7 @@ Type objective_function::operator()() mats, mats_save_hist, hist_shape_template); + #ifdef MP_VERBOSE Rcpp::Rcout << "Simulation history ..." << std::endl; @@ -3915,6 +4001,7 @@ Type objective_function::operator()() matrix value_column = values.block(0, 4, values.rows(), 1); ADREPORT(value_column) } + // 7 Calc the return of the objective function matrix ret; @@ -3931,5 +4018,6 @@ Type objective_function::operator()() #ifdef MP_VERBOSE Rcpp::Rcout << "======== end of objective function ========" << std::endl; #endif + return ret.coeff(0, 0); } diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 2f5b74d3..25d5884b 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -136,8 +136,9 @@ enum macpan2_func , MP2_PNORM = 57 // fwrap: pnorm(q, mean, sd) , MP2_INVLOGIT = 58 // fwrap: invlogit(x) , MP2_LOGIT = 59 // fwrap: logit(x) - , MP2_ASSIGN = 60 // fwrap: assign(x, i, j, v) - , MP2_UNPACK = 61 // fwrap: unpack(x, ...) + , MP2_STOP_IF_LT = 60 // fwrap: stop_if_lt(x, y) + , MP2_ASSIGN = 61 // fwrap: assign(x, i, j, v) + , MP2_UNPACK = 62 // fwrap: unpack(x, ...) }; enum macpan2_meth @@ -825,7 +826,6 @@ class ArgList // Rcpp::Rcout << "step f" << std::endl; // Rcpp::Rcout << "bad column vector" << std::endl; error_code = 501; - // break; // Exit the loop on error } } else if (mat.cols() == cols) { @@ -839,7 +839,6 @@ class ArgList { // Rcpp::Rcout << "bad row vector" << std::endl; error_code = 501; - // break; // Exit the loop on error } } else @@ -847,7 +846,6 @@ class ArgList // Rcpp::Rcout << "really bad" << std::endl; // Rcpp::Rcout << "step g" << std::endl; error_code = 501; - // break; // Exit the loop on error } if (error_code != 0) @@ -901,11 +899,6 @@ class ArgList 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 -// }; template class ExprEvaluator @@ -920,6 +913,7 @@ class ExprEvaluator ListOfIntVecs meth_int_vecs; ListOfIntVecs valid_int_vecs; vector valid_literals; + bool early_stop; public: // constructor @@ -951,6 +945,10 @@ class ExprEvaluator meth_int_vecs = meth_int_vecs_; valid_int_vecs = valid_int_vecs_; valid_literals = valid_literals_; + // Rcpp::Rcout << "setting early_stop to false" << std::endl; + early_stop = false; + // Rcpp::Rcout << "early_stop now equals: " << early_stop << std::endl; + // Rcpp::Rcout << "but GetEarlyStop() now returns: " << GetEarlyStop() << std::endl; strcpy(error_message, "OK"); }; @@ -964,6 +962,7 @@ class ExprEvaluator std::vector GetArgRows() { return arg_rows; }; std::vector GetArgCols() { return arg_cols; }; std::vector GetArgTypeInts() { return arg_type_ints; }; + bool GetEarlyStop() { return early_stop; }; // 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) { @@ -976,13 +975,17 @@ class ExprEvaluator arg_type_ints = a_type_ints; strcpy(error_message, message); }; + void SetEarlyStop() { + early_stop = true; + }; // evaluator matrix EvalExpr( const vector> &hist, // current simulation history 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 + int row = 0, // current expression parse table row being evaluated + bool simulate_block = false // evaluated within simulate block? ) { // total number of time steps in the simulation loop @@ -1124,7 +1127,7 @@ class ExprEvaluator { // otherwise, recursively descend into the parse tree // to pick out the arguments - args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i)); + args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i, simulate_block)); } // Check here if index2mats actually points at @@ -3091,6 +3094,55 @@ class ExprEvaluator } return m; + + // #' ## Control Flow + // #' + // #' Functions for controlling the flow of simulations. + // #' + // #' ### Functions + // #' + // #' * `stop_if_lt(x, y)` : Stop the simulation loop and return + // #' simulations if any element of `x` is less than the + // #' corresponding element in `y`. The elements of `y` are + // #' recycled to match those in `x` if necessary and if possible. + // #' + // #' ### Arguments + // #' + // #' * `x` : Matrix to compare with `y`. + // #' * `y` : Matrix to compare with `x`, with a default value of + // #' `0`. + // #' + // #' ### Return + // #' + // #' An empty matrix. This function is returned for its effect of + // #' breaking out of the simulation loop. + // #' + case MP2_STOP_IF_LT: + m = args.get_as_mat(0); + rows = m.rows(); + cols = m.cols(); + if (n == 1) { + m1 = matrix::Zero(rows, cols); + } else { + v1.push_back(1); + args = args.recycle_to_shape(v1, rows, cols); + m1 = args.get_as_mat(1); + } + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + if (m.coeff(i, j) < m1.coeff(i, j)) { + //Rcpp::Rcout << "TIME: " << t << std::endl; + //is_finite_mat = GetEarlyStop(); + //Rcpp::Rcout << "Early Stop BEFORE SET: " << is_finite_mat << std::endl; + SetEarlyStop(); + //is_finite_mat = GetEarlyStop(); + //Rcpp::Rcout << "Early Stop AFTER SET: " << is_finite_mat << std::endl; + return m2; // empty matrix + } + } + } + return m2; // empty matrix + // #' ## Assign (deprecated) // #' // #' Assign values to a subset of the elements in a matrix. @@ -3197,7 +3249,7 @@ class ExprEvaluator valid_vars.m_matrices[matIndex].coeffRef(rowIndex, colIndex) = args[3].coeff(k, 0); } return m2; // empty matrix - + // #' ## Unpack (deprecated) // #' // #' Unpack elements of a matrix into smaller matrices. @@ -3476,6 +3528,7 @@ class MatAssigner { int func_int = exprEvaluator.GetFuncInt(); \ int time_int = exprEvaluator.GetTimeInt(); \ const char *err_msg = exprEvaluator.GetErrorMessage(); \ + bool early_stop = exprEvaluator.GetEarlyStop(); \ REPORT(error); \ REPORT(expr_row); \ REPORT(func_int); \ @@ -3490,6 +3543,7 @@ class MatAssigner { logfile << "Expression row = " << expr_row << std::endl; \ logfile << "Function code = " << func_int << std::endl; \ logfile << "Time step = " << time_int << std::endl; \ + logfile << "Early stop = " << early_stop << std::endl; \ logfile.close(); \ } @@ -3666,6 +3720,14 @@ Type objective_function::operator()() for (int i = 0; i < n; i++) mats.m_matrices[r_mat_id[i]].coeffRef(r_row_id[i], r_col_id[i]) = random[r_par_id[i]]; + // in case you need to stop early you can use this to fill the remaining + // iterations of the simulation history + ListOfMatrices empty_mats_list(mats); + matrix empty_matrix; + for (int i = 0; i < empty_mats_list.m_matrices.size(); i++) { + empty_mats_list.m_matrices[i] = empty_matrix; + } + // Simulation history /// each element of this history 'vector' is a list of the matrices /// in the model at a particular point in history @@ -3720,7 +3782,7 @@ Type objective_function::operator()() if (expr_sim_block[i] == 1) { SIMULATE { result = exprEvaluator.EvalExpr( - simulation_history, 0, mats, p_table_row); + simulation_history, 0, mats, p_table_row, true); } } else @@ -3768,23 +3830,35 @@ Type objective_function::operator()() Rcpp::Rcout << "Eval expression --- " << i << std::endl; Rcpp::Rcout << "expr_num_p_table_rows[i] " << expr_num_p_table_rows[expr_index + i] << std::endl; #endif - matrix result; - if (expr_sim_block[i] == 1) { - SIMULATE { + //Rcpp::Rcout << "GetEarlyStop(): " << exprEvaluator.GetEarlyStop() << std::endl; + if (!exprEvaluator.GetEarlyStop()) { + matrix result; + if (expr_sim_block[i] == 1) { + SIMULATE { + result = exprEvaluator.EvalExpr( + simulation_history + , k + 1 + , mats + , p_table_row2 + , true + ); + } + } + else result = exprEvaluator.EvalExpr( - simulation_history, k + 1, mats, p_table_row2); + simulation_history + , k + 1 + , mats + , p_table_row2 + , false + ); + + if (exprEvaluator.GetErrorCode()) { + REPORT_ERROR + return 0.0; } + matAssigner.matAssign(result, mats, a_table_row2); } - else - result = exprEvaluator.EvalExpr( - 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); p_table_row2 += expr_num_p_table_rows[expr_index + i]; a_table_row2 += assign_num_a_table_rows[expr_index + i]; @@ -3795,13 +3869,24 @@ Type objective_function::operator()() Rcpp::Rcout << "mats = " << mats.m_matrices[ii] << std::endl; #endif } - // simulation_history[k+1] = mats; - UpdateSimulationHistory( - simulation_history, - k + 1, - mats, - mats_save_hist, - hist_shape_template); + if (exprEvaluator.GetEarlyStop()) { + UpdateSimulationHistory( + simulation_history, + k + 1, + empty_mats_list, + mats_save_hist, + hist_shape_template + ); + } else { + // simulation_history[k+1] = mats; + UpdateSimulationHistory( + simulation_history, + k + 1, + mats, + mats_save_hist, + hist_shape_template + ); + } } p_table_row = p_table_row2; a_table_row = a_table_row2; @@ -3818,7 +3903,7 @@ Type objective_function::operator()() if (expr_sim_block[i] == 1) { SIMULATE { result = exprEvaluator.EvalExpr( - simulation_history, time_steps + 1, mats, p_table_row); + simulation_history, time_steps + 1, mats, p_table_row, true); } } else { @@ -3845,6 +3930,7 @@ Type objective_function::operator()() mats, mats_save_hist, hist_shape_template); + #ifdef MP_VERBOSE Rcpp::Rcout << "Simulation history ..." << std::endl; @@ -3916,6 +4002,7 @@ Type objective_function::operator()() matrix value_column = values.block(0, 4, values.rows(), 1); ADREPORT(value_column) } + // 7 Calc the return of the objective function matrix ret; @@ -3932,5 +4019,6 @@ Type objective_function::operator()() #ifdef MP_VERBOSE Rcpp::Rcout << "======== end of objective function ========" << std::endl; #endif + return ret.coeff(0, 0); } From c455f159a4ee2a73a36cb2031619a5279c624d92 Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Thu, 26 Jun 2025 15:44:38 -0400 Subject: [PATCH 2/2] testing --- misc/dev/dev.cpp | 6 ------ src/macpan2.cpp | 6 ------ tests/testthat/test-stop-if.R | 16 ++++++++++++++++ 3 files changed, 16 insertions(+), 12 deletions(-) create mode 100644 tests/testthat/test-stop-if.R diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp index e362799f..e63f40d1 100644 --- a/misc/dev/dev.cpp +++ b/misc/dev/dev.cpp @@ -3130,12 +3130,7 @@ class ExprEvaluator for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { if (m.coeff(i, j) < m1.coeff(i, j)) { - //Rcpp::Rcout << "TIME: " << t << std::endl; - //is_finite_mat = GetEarlyStop(); - //Rcpp::Rcout << "Early Stop BEFORE SET: " << is_finite_mat << std::endl; SetEarlyStop(); - //is_finite_mat = GetEarlyStop(); - //Rcpp::Rcout << "Early Stop AFTER SET: " << is_finite_mat << std::endl; return m2; // empty matrix } } @@ -3829,7 +3824,6 @@ Type objective_function::operator()() Rcpp::Rcout << "Eval expression --- " << i << std::endl; Rcpp::Rcout << "expr_num_p_table_rows[i] " << expr_num_p_table_rows[expr_index + i] << std::endl; #endif - //Rcpp::Rcout << "GetEarlyStop(): " << exprEvaluator.GetEarlyStop() << std::endl; if (!exprEvaluator.GetEarlyStop()) { matrix result; if (expr_sim_block[i] == 1) { diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 25d5884b..da28aabe 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -3131,12 +3131,7 @@ class ExprEvaluator for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { if (m.coeff(i, j) < m1.coeff(i, j)) { - //Rcpp::Rcout << "TIME: " << t << std::endl; - //is_finite_mat = GetEarlyStop(); - //Rcpp::Rcout << "Early Stop BEFORE SET: " << is_finite_mat << std::endl; SetEarlyStop(); - //is_finite_mat = GetEarlyStop(); - //Rcpp::Rcout << "Early Stop AFTER SET: " << is_finite_mat << std::endl; return m2; // empty matrix } } @@ -3830,7 +3825,6 @@ Type objective_function::operator()() Rcpp::Rcout << "Eval expression --- " << i << std::endl; Rcpp::Rcout << "expr_num_p_table_rows[i] " << expr_num_p_table_rows[expr_index + i] << std::endl; #endif - //Rcpp::Rcout << "GetEarlyStop(): " << exprEvaluator.GetEarlyStop() << std::endl; if (!exprEvaluator.GetEarlyStop()) { matrix result; if (expr_sim_block[i] == 1) { diff --git a/tests/testthat/test-stop-if.R b/tests/testthat/test-stop-if.R new file mode 100644 index 00000000..8e80c858 --- /dev/null +++ b/tests/testthat/test-stop-if.R @@ -0,0 +1,16 @@ +test_that("simulations can be stopped early", { + exprs = list(x ~ x - 2, sims = dummy ~ stop_if_lt(x), y ~ y - 2) + def = list(x = 10, y = 10) + ts = 10 + spec = mp_tmb_model_spec(during = exprs, default = def) + sim = mp_simulator(spec, ts, c("x", "y")) + + expect_equal( + mp_final(sim)$value + , c(-2, 0) + ) + expect_equal( + mp_trajectory(sim)$value + , rep(seq(8, 0, by = -2), each = 2) + ) +})