diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 19c2ea4aea165f0f12d9e6db892f275beb56b46a..4c55c28232ccfcc97f0e46fe6ff59ce1db6af9b4 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -188,10 +188,8 @@ if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_fores
yT = NaN(ny, 1);
end
yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
- model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']);
- [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
- residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll, options_.threads.perfect_foresight_problem);
+ residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
if max(abs(residuals))< options_.dynatol.f
oo_.deterministic_simulation.status = 1;
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index b336c4d0356fc576888583721bf2c8d05d66e0a1..a81aa300e62536c152069cb683d3641718537d95 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -141,10 +141,8 @@ if nargout>1
yT = NaN(ny, 1);
end
yy = oo_.endo_simul(:,M_.maximum_lag+(1:periods));
- model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']);
- [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
- residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll, options_.threads.perfect_foresight_problem);
+ residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
end
maxerror = max(max(abs(residuals)));
end
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 58efa011bba930439c8436185390f71d52ac7dc0..b79c8a4b45d3bbc5c56347f80bfcd76662246cbb 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -59,15 +59,12 @@ if verbose
skipline()
end
-model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']);
-[nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
-
h1 = clock;
for iter = 1:options.simul.maxit
h2 = clock;
- [res, A] = perfect_foresight_problem(y, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll, options.threads.perfect_foresight_problem);
+ [res, A] = perfect_foresight_problem(y, y0, yT, exogenousvariables, M.params, steadystate, periods, M, options);
if options.endogenous_terminal_period && iter > 1
for it = 1:periods
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index 07d1804b35ad616e2066ede561d79c9cd4ae6af1..e52a833092abb71dbc77a8b3c73569c26684d505 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -53,10 +53,7 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
eq_index);
else
- model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']);
- [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
-
- [y, check] = dynare_solve(@perfect_foresight_problem,z(:), options, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, options.periods, M.endo_nbr, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll, options.threads.perfect_foresight_problem);
+ [y, check] = dynare_solve(@perfect_foresight_problem,z(:), options, y0, yT, exogenousvariables, M.params, steadystate, options.periods, M, options);
end
if all(imag(y)<.1*options.dynatol.x)
diff --git a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
index 5a1447295ddc54db0d4d1ff0f8263028447a2ea9..107dc3caaf1d31719e2d803af98d5b45342edf46 100644
--- a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
+++ b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
@@ -28,86 +28,78 @@
void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
- if (nlhs < 1 || nlhs > 2 || nrhs != 19)
- mexErrMsgTxt("Must have 19 input arguments and 1 or 2 output arguments");
+ if (nlhs < 1 || nlhs > 2 || nrhs != 9)
+ mexErrMsgTxt("Must have 9 input arguments and 1 or 2 output arguments");
bool compute_jacobian = nlhs == 2;
// Give explicit names to input arguments
const mxArray *y_mx = prhs[0];
- const mxArray *basename_mx = prhs[1];
- const mxArray *ntt_mx = prhs[2];
- const mxArray *y0_mx = prhs[3];
- const mxArray *yT_mx = prhs[4];
- const mxArray *exo_path_mx = prhs[5];
- const mxArray *params_mx = prhs[6];
- const mxArray *steady_state_mx = prhs[7];
- const mxArray *periods_mx = prhs[8];
- const mxArray *ny_mx = prhs[9];
- const mxArray *maximum_lag_mx = prhs[10];
- const mxArray *maximum_endo_lag_mx = prhs[11];
- const mxArray *lead_lag_incidence_mx = prhs[12];
- const mxArray *nzij_pred_mx = prhs[13];
- const mxArray *nzij_current_mx = prhs[14];
- const mxArray *nzij_fwrd_mx = prhs[15];
- const mxArray *has_external_function_mx = prhs[16];
- const mxArray *use_dll_mx = prhs[17];
- const mxArray *num_threads_mx = prhs[18];
-
- // Check input and map it to local variables
- if (!(mxIsChar(basename_mx) && mxGetM(basename_mx) == 1))
- mexErrMsgTxt("basename should be a character string");
- char *basename = mxArrayToString(basename_mx);
-
- if (!(mxIsScalar(ntt_mx) && mxIsNumeric(ntt_mx)))
- mexErrMsgTxt("ntt should be a numeric scalar");
- size_t ntt = mxGetScalar(ntt_mx);
-
- if (!(mxIsScalar(periods_mx) && mxIsNumeric(periods_mx)))
- mexErrMsgTxt("periods should be a numeric scalar");
- mwIndex periods = static_cast<mwIndex>(mxGetScalar(periods_mx));
-
- if (!(mxIsScalar(ny_mx) && mxIsNumeric(ny_mx)))
- mexErrMsgTxt("ny should be a numeric scalar");
- mwIndex ny = static_cast<mwIndex>(mxGetScalar(ny_mx));
-
- if (!(mxIsScalar(maximum_lag_mx) && mxIsNumeric(maximum_lag_mx)))
- mexErrMsgTxt("maximum_lag should be a numeric scalar");
+ const mxArray *y0_mx = prhs[1];
+ const mxArray *yT_mx = prhs[2];
+ const mxArray *exo_path_mx = prhs[3];
+ const mxArray *params_mx = prhs[4];
+ const mxArray *steady_state_mx = prhs[5];
+ const mxArray *periods_mx = prhs[6];
+ const mxArray *M_mx = prhs[7];
+ const mxArray *options_mx = prhs[8];
+
+ // Extract various fields from M_
+ const mxArray *basename_mx = mxGetField(M_mx, 0, "fname");
+ if (!(basename_mx && mxIsChar(basename_mx) && mxGetM(basename_mx) == 1))
+ mexErrMsgTxt("M_.fname should be a character string");
+ std::string basename{mxArrayToString(basename_mx)};
+
+ const mxArray *endo_nbr_mx = mxGetField(M_mx, 0, "endo_nbr");
+ if (!(endo_nbr_mx && mxIsScalar(endo_nbr_mx) && mxIsNumeric(endo_nbr_mx)))
+ mexErrMsgTxt("M_.endo_nbr should be a numeric scalar");
+ mwIndex ny = static_cast<mwIndex>(mxGetScalar(endo_nbr_mx));
+
+ const mxArray *maximum_lag_mx = mxGetField(M_mx, 0, "maximum_lag");
+ if (!(maximum_lag_mx && mxIsScalar(maximum_lag_mx) && mxIsNumeric(maximum_lag_mx)))
+ mexErrMsgTxt("M_.maximum_lag should be a numeric scalar");
mwIndex maximum_lag = static_cast<mwIndex>(mxGetScalar(maximum_lag_mx));
- if (!(mxIsScalar(maximum_endo_lag_mx) && mxIsNumeric(maximum_endo_lag_mx)))
- mexErrMsgTxt("maximum_endo_lag should be a numeric scalar");
+ const mxArray *maximum_endo_lag_mx = mxGetField(M_mx, 0, "maximum_endo_lag");
+ if (!(maximum_endo_lag_mx && mxIsScalar(maximum_endo_lag_mx) && mxIsNumeric(maximum_endo_lag_mx)))
+ mexErrMsgTxt("M_.maximum_endo_lag should be a numeric scalar");
mwIndex maximum_endo_lag = static_cast<mwIndex>(mxGetScalar(maximum_endo_lag_mx));
- if (!(mxIsDouble(y_mx) && mxGetM(y_mx) == static_cast<size_t>(ny*periods) && mxGetN(y_mx) == 1))
- mexErrMsgTxt("y should be a double precision column-vector of ny*periods elements");
- const double *y = mxGetPr(y_mx);
-
- if (!(mxIsDouble(y0_mx) && mxGetM(y0_mx) == static_cast<size_t>(ny) && mxGetN(y0_mx) == 1))
- mexErrMsgTxt("y0 should be a double precision column-vector of ny elements");
- const double *y0 = mxGetPr(y0_mx);
+ const mxArray *dynamic_tmp_nbr_mx = mxGetField(M_mx, 0, "dynamic_tmp_nbr");
+ if (!(dynamic_tmp_nbr_mx && mxIsDouble(dynamic_tmp_nbr_mx) && mxGetNumberOfElements(dynamic_tmp_nbr_mx) >= 2))
+ mexErrMsgTxt("M_.dynamic_tmp_nbr should be a double array of at least 2 elements");
+ size_t ntt = mxGetPr(dynamic_tmp_nbr_mx)[0] + mxGetPr(dynamic_tmp_nbr_mx)[1];
- if (!(mxIsDouble(yT_mx) && mxGetM(yT_mx) == static_cast<size_t>(ny) && mxGetN(yT_mx) == 1))
- mexErrMsgTxt("yT should be a double precision column-vector of ny elements");
- const double *yT = mxGetPr(yT_mx);
+ const mxArray *lead_lag_incidence_mx = mxGetField(M_mx, 0, "lead_lag_incidence");
+ if (!(lead_lag_incidence_mx && mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == static_cast<size_t>(2+maximum_endo_lag)
+ && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(ny)))
+ mexErrMsgTxt("M_.lead_lag_incidence should be a double precision matrix with 2+M_.maximum_endo_lag rows and M_.endo_nbr columns");
+ const double *lead_lag_incidence = mxGetPr(lead_lag_incidence_mx);
- if (!(mxIsDouble(exo_path_mx) && mxGetM(exo_path_mx) >= static_cast<size_t>(periods+maximum_lag)))
- mexErrMsgTxt("exo_path should be a double precision matrix with at least periods+maximum_lag rows");
- mwIndex nx = static_cast<mwIndex>(mxGetN(exo_path_mx));
- size_t nb_row_x = mxGetM(exo_path_mx);
- const double *exo_path = mxGetPr(exo_path_mx);
+ const mxArray *has_external_function_mx = mxGetField(M_mx, 0, "has_external_function");
+ if (!(has_external_function_mx && mxIsLogicalScalar(has_external_function_mx)))
+ mexErrMsgTxt("M_.has_external_function should be a logical scalar");
+ bool has_external_function = static_cast<bool>(mxGetScalar(has_external_function_mx));
- if (!(mxIsDouble(params_mx) && mxGetN(params_mx) == 1))
- mexErrMsgTxt("params should be a double precision column-vector");
- const double *params = mxGetPr(params_mx);
+ // Extract various fields from options_
+ const mxArray *use_dll_mx = mxGetField(options_mx, 0, "use_dll");
+ if (!(use_dll_mx && mxIsLogicalScalar(use_dll_mx)))
+ mexErrMsgTxt("options_.use_dll should be a logical scalar");
+ bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
- if (!(mxIsDouble(steady_state_mx) && mxGetN(steady_state_mx) == 1))
- mexErrMsgTxt("steady_state should be a double precision column-vector");
- const double *steady_state = mxGetPr(steady_state_mx);
+ const mxArray *threads_mx = mxGetField(options_mx, 0, "threads");
+ if (!threads_mx)
+ mexErrMsgTxt("Can't find field options_.threads");
+ const mxArray *num_threads_mx = mxGetField(threads_mx, 0, "perfect_foresight_problem");
+ if (!(num_threads_mx && mxIsScalar(num_threads_mx) && mxIsNumeric(num_threads_mx)))
+ mexErrMsgTxt("options_.threads.perfect_foresight_problem should be a numeric scalar");
+ int num_threads = static_cast<int>(mxGetScalar(num_threads_mx));
- if (!(mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == static_cast<size_t>(2+maximum_endo_lag)
- && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(ny)))
- mexErrMsgTxt("lead_lag_incidence should be a double precision matrix with 2+maximum_endo_lag rows and endo_nbr columns");
- const double *lead_lag_incidence = mxGetPr(lead_lag_incidence_mx);
+ // Call <model>.dynamic_g1_nz
+ mxArray *g1_nz_plhs[3];
+ mexCallMATLAB(3, g1_nz_plhs, 0, nullptr, (basename + ".dynamic_g1_nz").c_str());
+ const mxArray *nzij_pred_mx = g1_nz_plhs[0];
+ const mxArray *nzij_current_mx = g1_nz_plhs[1];
+ const mxArray *nzij_fwrd_mx = g1_nz_plhs[2];
if (!(mxIsInt32(nzij_pred_mx) && mxGetN(nzij_pred_mx) == 2))
mexErrMsgTxt("nzij_pred should be an int32 matrix with 2 columns");
@@ -136,17 +128,36 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
const int32_T *nzij_fwrd = static_cast<const int32_T *>(mxGetData(nzij_fwrd_mx));
#endif
- if (!(mxIsLogicalScalar(has_external_function_mx)))
- mexErrMsgTxt("has_external_function should be a logical scalar");
- bool has_external_function = static_cast<bool>(mxGetScalar(has_external_function_mx));
+ // Check other input and map it to local variables
+ if (!(mxIsScalar(periods_mx) && mxIsNumeric(periods_mx)))
+ mexErrMsgTxt("periods should be a numeric scalar");
+ mwIndex periods = static_cast<mwIndex>(mxGetScalar(periods_mx));
- if (!(mxIsLogicalScalar(use_dll_mx)))
- mexErrMsgTxt("use_dll should be a logical scalar");
- bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
+ if (!(mxIsDouble(y_mx) && mxGetM(y_mx) == static_cast<size_t>(ny*periods) && mxGetN(y_mx) == 1))
+ mexErrMsgTxt("y should be a double precision column-vector of M_.endo_nbr*periods elements");
+ const double *y = mxGetPr(y_mx);
- if (!(mxIsScalar(num_threads_mx) && mxIsNumeric(num_threads_mx)))
- mexErrMsgTxt("num_threads should be a numeric scalar");
- int num_threads = static_cast<int>(mxGetScalar(num_threads_mx));
+ if (!(mxIsDouble(y0_mx) && mxGetM(y0_mx) == static_cast<size_t>(ny) && mxGetN(y0_mx) == 1))
+ mexErrMsgTxt("y0 should be a double precision column-vector of M_.endo_nbr elements");
+ const double *y0 = mxGetPr(y0_mx);
+
+ if (!(mxIsDouble(yT_mx) && mxGetM(yT_mx) == static_cast<size_t>(ny) && mxGetN(yT_mx) == 1))
+ mexErrMsgTxt("yT should be a double precision column-vector of M_.endo_nbr elements");
+ const double *yT = mxGetPr(yT_mx);
+
+ if (!(mxIsDouble(exo_path_mx) && mxGetM(exo_path_mx) >= static_cast<size_t>(periods+maximum_lag)))
+ mexErrMsgTxt("exo_path should be a double precision matrix with at least periods+M_.maximum_lag rows");
+ mwIndex nx = static_cast<mwIndex>(mxGetN(exo_path_mx));
+ size_t nb_row_x = mxGetM(exo_path_mx);
+ const double *exo_path = mxGetPr(exo_path_mx);
+
+ if (!(mxIsDouble(params_mx) && mxGetN(params_mx) == 1))
+ mexErrMsgTxt("params should be a double precision column-vector");
+ const double *params = mxGetPr(params_mx);
+
+ if (!(mxIsDouble(steady_state_mx) && mxGetN(steady_state_mx) == 1))
+ mexErrMsgTxt("steady_state should be a double precision column-vector");
+ const double *steady_state = mxGetPr(steady_state_mx);
// Allocate output matrices
plhs[0] = mxCreateDoubleMatrix(periods*ny, 1, mxREAL);