From e4a4d2d8e6cc128e9ce3a2382e1b534a12d094aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 11 Oct 2023 16:50:18 -0400 Subject: [PATCH] Bytecode MEX: get M_ and options_ through input arguments rather than as global variables --- matlab/dyn_ramsey_static.m | 4 +- matlab/evaluate_static_model.m | 4 +- matlab/evaluate_steady_state.m | 20 ++++---- matlab/model_diagnostics.m | 8 +-- .../add_auxiliary_variables_to_steadystate.m | 7 ++- .../det_cond_forecast.m | 12 ++--- .../perfect_foresight_solver.m | 2 +- .../perfect_foresight_solver_core.m | 6 +-- matlab/print_bytecode_dynamic_model.m | 4 +- matlab/print_bytecode_static_model.m | 4 +- matlab/stochastic_solvers.m | 2 +- mex/sources/bytecode/BasicSymbolTable.cc | 8 +-- mex/sources/bytecode/BasicSymbolTable.hh | 6 ++- mex/sources/bytecode/bytecode.cc | 50 ++++++++++--------- 14 files changed, 71 insertions(+), 66 deletions(-) diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m index 9409b3ae9..32d399b1d 100644 --- a/matlab/dyn_ramsey_static.m +++ b/matlab/dyn_ramsey_static.m @@ -137,7 +137,7 @@ end % Compute the value of the Lagrange multipliers that minimizes the norm of the % residuals, given the other endogenous if options_.bytecode - res = bytecode('static', xx, exo_ss, M.params, 'evaluate'); + res = bytecode('static', M, options, xx, exo_ss, M.params, 'evaluate'); else res = feval([M.fname '.sparse.static_resid'], xx, exo_ss, M.params); end @@ -167,7 +167,7 @@ end function result = check_static_model(ys,exo_ss,M,options_) result = false; if (options_.bytecode) - [res, ~] = bytecode('static', ys, exo_ss, M.params, 'evaluate'); + [res, ~] = bytecode('static', M, options, ys, exo_ss, M.params, 'evaluate'); else res = feval([M.fname '.sparse.static_resid'], ys, exo_ss, M.params); end diff --git a/matlab/evaluate_static_model.m b/matlab/evaluate_static_model.m index 71d09362f..b05c65d41 100644 --- a/matlab/evaluate_static_model.m +++ b/matlab/evaluate_static_model.m @@ -40,10 +40,10 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt check1 = 0; if options.bytecode if nargout<3 - [residuals]= bytecode('evaluate','static',ys,... + [residuals]= bytecode('evaluate', 'static', M, options, ys, ... exo_ss, params, ys, 1); else - [residuals, junk]= bytecode('evaluate','static',ys,... + [residuals, junk]= bytecode('evaluate', 'static', M, options, ys, ... exo_ss, params, ys, 1); jacob = junk.g1; end diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m index 49a33b2b9..0fa2cb183 100644 --- a/matlab/evaluate_steady_state.m +++ b/matlab/evaluate_steady_state.m @@ -395,9 +395,9 @@ elseif options.bytecode if options.solve_algo >= 5 && options.solve_algo <= 8 try if options.block - ys = bytecode('static', 'block_decomposed', ys_init, exo_ss, params); + ys = bytecode('static', 'block_decomposed', M, options, ys_init, exo_ss, params); else - ys = bytecode('static', ys_init, exo_ss, params); + ys = bytecode('static', M, options, ys_init, exo_ss, params); end catch ME if options.verbosity >= 1 @@ -416,7 +416,7 @@ elseif options.bytecode [ys(mfs_idx), errorflag] = dynare_solve(@block_bytecode_mfs_steadystate, ... ys(mfs_idx), options.steady.maxit, ... options.solve_tolf, options.solve_tolx, ... - options, b, ys, exo_ss, params, T, M); + options, b, ys, exo_ss, params, T, M, options); if errorflag check = 1; break @@ -425,7 +425,7 @@ elseif options.bytecode % Compute endogenous if the block is of type evaluate forward/backward or if there are recursive variables in a solve block. % Also update the temporary terms vector (needed for the dynare_solve case) try - [~, ~, ys, T] = bytecode(ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ... + [~, ~, ys, T] = bytecode(M, options, ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ... 'block_decomposed', ['block = ' int2str(b)]); catch ME if options.verbosity >= 1 @@ -438,7 +438,7 @@ elseif options.bytecode else [ys, check] = dynare_solve(@bytecode_steadystate, ys_init, ... options.steady.maxit, options.solve_tolf, options.solve_tolx, ... - options, exo_ss, params); + options, exo_ss, params, M, options); end end @@ -465,7 +465,7 @@ if M.static_and_dynamic_models_differ if options.bytecode z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1); zx = repmat(exo_ss', M.maximum_lead + M.maximum_lag + 1, 1); - r = bytecode('dynamic','evaluate', z, zx, params, ys, 1); + r = bytecode('dynamic','evaluate', M, options, z, zx, params, ys, 1); else r = feval([M.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys); end @@ -513,13 +513,13 @@ y_all(mfs_idx) = y; M.block_structure_stat.block(b).g1_sparse_colval, ... M.block_structure_stat.block(b).g1_sparse_colptr, T); -function [r, g1] = bytecode_steadystate(y, exo, params) +function [r, g1] = bytecode_steadystate(y, exo, params, M, options) % Wrapper around the static file, for bytecode (without block) -[r, g1] = bytecode(y, exo, params, y, 1, exo, 'evaluate', 'static'); +[r, g1] = bytecode(M, options, y, exo, params, y, 1, exo, 'evaluate', 'static'); -function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M) +function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M, options) % Wrapper around the static files, for block without bytecode mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end); y_all(mfs_idx) = y; -[r, g1] = bytecode(y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', 'block_decomposed', ['block = ' int2str(b) ]); +[r, g1] = bytecode(M, options, y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', 'block_decomposed', ['block = ' int2str(b) ]); g1 = g1(:,end-M.block_structure_stat.block(b).mfs+1:end); % Make Jacobian square if mfs>0 diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m index ca4a0846d..d4065a4bc 100644 --- a/matlab/model_diagnostics.m +++ b/matlab/model_diagnostics.m @@ -147,10 +147,10 @@ exo = [oo_.exo_steady_state; oo_.exo_det_steady_state]; for b=1:nb if options_.bytecode if nb == 1 - [res, jacob] = bytecode(dr.ys, exo, M_.params, dr.ys, 1, exo, ... + [res, jacob] = bytecode(M_, options_, dr.ys, exo, M_.params, dr.ys, 1, exo, ... 'evaluate', 'static'); else - [res, jacob] = bytecode(dr.ys, exo, M_.params, dr.ys, 1, exo, ... + [res, jacob] = bytecode(M_, options_, dr.ys, exo, M_.params, dr.ys, 1, exo, ... 'evaluate', 'static', 'block_decomposed', ['block=' ... int2str(b)]); end @@ -275,7 +275,7 @@ z = repmat(dr.ys,1,klen); if options_.order == 1 if (options_.bytecode) - [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + [~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ... M_.params, dr.ys, 1); jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; else @@ -284,7 +284,7 @@ if options_.order == 1 end elseif options_.order >= 2 if (options_.bytecode) - [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + [~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ... M_.params, dr.ys, 1); jacobia_ = [loc_dr.g1 loc_dr.g1_x]; else diff --git a/matlab/partial_information/add_auxiliary_variables_to_steadystate.m b/matlab/partial_information/add_auxiliary_variables_to_steadystate.m index 3e67613af..0ebbecda3 100644 --- a/matlab/partial_information/add_auxiliary_variables_to_steadystate.m +++ b/matlab/partial_information/add_auxiliary_variables_to_steadystate.m @@ -2,7 +2,7 @@ function ys1 = add_auxiliary_variables_to_steadystate(ys,aux_vars,fname, ... exo_steady_state, exo_det_steady_state,params, byte_code) % Add auxiliary variables to the steady state vector -% Copyright © 2009-2020 Dynare Team +% Copyright © 2009-2023 Dynare Team % % This file is part of Dynare. % @@ -18,12 +18,15 @@ function ys1 = add_auxiliary_variables_to_steadystate(ys,aux_vars,fname, ... % % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . + +global M_ options_ + n = length(aux_vars); ys1 = [ys;zeros(n,1)]; for i=1:n+1 if byte_code - res = bytecode('static','evaluate',ys1,... + res = bytecode('static','evaluate', M_, options_, ys1,... [exo_steady_state; ... exo_det_steady_state],params); else diff --git a/matlab/perfect-foresight-models/det_cond_forecast.m b/matlab/perfect-foresight-models/det_cond_forecast.m index 261150218..9cedec866 100644 --- a/matlab/perfect-foresight-models/det_cond_forecast.m +++ b/matlab/perfect-foresight-models/det_cond_forecast.m @@ -162,9 +162,9 @@ else save_options_dynatol_f = options_.dynatol.f; options_.dynatol.f = 1e-7; if options_.block - [endo, exo] = bytecode('extended_path', plan, 'block_decomposed', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods); + [endo, exo] = bytecode('extended_path', plan, 'block_decomposed', M_, options_, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods); else - [endo, exo] = bytecode('extended_path', plan, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods); + [endo, exo] = bytecode('extended_path', plan, M_, options_, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods); end options_.dynatol.f = save_options_dynatol_f; @@ -473,9 +473,9 @@ if pf && ~surprise data1 = M_; if (options_.bytecode) if options_.block - [~, data1]= bytecode('dynamic','block_decomposed','evaluate', z, zx, M_.params, oo_.steady_state, k, data1); + [~, data1]= bytecode('dynamic','block_decomposed','evaluate', M_, options_, z, zx, M_.params, oo_.steady_state, k, data1); else - [~, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1); + [~, data1]= bytecode('dynamic','evaluate', M_, options_, z, zx, M_.params, oo_.steady_state, k, data1); end else [~, g1b] = feval([M_.fname '.dynamic'], z', zx, M_.params, oo_.steady_state, k); @@ -740,9 +740,9 @@ else data1 = M_; if (options_.bytecode) if options_.block - [~, data1]= bytecode('dynamic','block_decomposed','evaluate', z, zx, M_.params, oo_.steady_state, k, data1); + [~, data1]= bytecode('dynamic','block_decomposed','evaluate', M_, options_, z, zx, M_.params, oo_.steady_state, k, data1); else - [~, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1); + [~, data1]= bytecode('dynamic','evaluate', M_, options_, z, zx, M_.params, oo_.steady_state, k, data1); end else [~, g1b] = feval([M_.fname '.dynamic'], z', zx, M_.params, oo_.steady_state, k); diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index f010e9081..29ab48c66 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -411,7 +411,7 @@ function maxerror = recompute_maxerror(endo_simul, oo_, M_, options_) % Computes ∞-norm of residuals for a given path of endogenous, % given the exogenous path and parameters in oo_ and M_ if options_.bytecode - residuals = bytecode('dynamic', 'evaluate', endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1); + residuals = bytecode('dynamic', 'evaluate', M_, options_, endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1); else ny = size(endo_simul, 1); periods = size(endo_simul, 2) - M_.maximum_lag - M_.maximum_lead; diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m index a77514759..c37b77990 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m @@ -62,7 +62,7 @@ if options_.block end if options_.bytecode try - y = bytecode('dynamic', 'block_decomposed', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods); + y = bytecode('dynamic', 'block_decomposed', M_, options_, oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods); success = true; catch ME if options_.verbosity >= 1 @@ -77,7 +77,7 @@ if options_.block else if options_.bytecode try - y = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods); + y = bytecode('dynamic', M_, options_, oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods); success = true; catch ME if options_.verbosity >= 1 @@ -131,7 +131,7 @@ end % Some solvers do not compute the maximum error, so do it here if needed if nargout > 2 && isempty(maxerror) if options_.bytecode - residuals = bytecode('dynamic', 'evaluate', y, oo_.exo_simul, M_.params, oo_.steady_state, periods); + residuals = bytecode('dynamic', 'evaluate', M_, options_, y, oo_.exo_simul, M_.params, oo_.steady_state, periods); else ny = size(y, 1); if M_.maximum_lag > 0 diff --git a/matlab/print_bytecode_dynamic_model.m b/matlab/print_bytecode_dynamic_model.m index ee4adee9b..eef3b6363 100644 --- a/matlab/print_bytecode_dynamic_model.m +++ b/matlab/print_bytecode_dynamic_model.m @@ -35,9 +35,9 @@ if options_.bytecode zx = repmat([oo_.exo_steady_state; oo_.exo_det_steady_state]', M_.maximum_lead + M_.maximum_lag + 1, 1); if options_.block - bytecode('print', 'dynamic', 'block_decomposed', z, zx, M_.params, oo_.steady_state, 1); + bytecode('print', 'dynamic', 'block_decomposed', M_, options_, z, zx, M_.params, oo_.steady_state, 1); else - bytecode('print', 'dynamic', z, zx, M_.params, oo_.steady_state, 1); + bytecode('print', 'dynamic', M_, options_, z, zx, M_.params, oo_.steady_state, 1); end else disp('You have to use bytecode option in model command to use print_bytecode_dynamic_model'); diff --git a/matlab/print_bytecode_static_model.m b/matlab/print_bytecode_static_model.m index ba22848f5..5e1279bd9 100644 --- a/matlab/print_bytecode_static_model.m +++ b/matlab/print_bytecode_static_model.m @@ -32,9 +32,9 @@ global M_ options_ oo_ if options_.bytecode if options_.block - bytecode('print', 'static', 'block_decomposed', oo_.steady_state, [oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params); + bytecode('print', 'static', 'block_decomposed', M_, options_, oo_.steady_state, [oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params); else - bytecode('print', 'static', oo_.steady_state, [oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params); + bytecode('print', 'static', M_, options_, oo_.steady_state, [oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params); end else disp('You have to use bytecode option in model command to use print_bytecode_static_model'); diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m index 7eb01dc88..c86484df5 100644 --- a/matlab/stochastic_solvers.m +++ b/matlab/stochastic_solvers.m @@ -103,7 +103,7 @@ it_ = M_.maximum_lag + 1; z = repmat(dr.ys,1,klen); if local_order == 1 if (options_.bytecode) - [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... + [~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ... M_.params, dr.ys, 1); jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; else diff --git a/mex/sources/bytecode/BasicSymbolTable.cc b/mex/sources/bytecode/BasicSymbolTable.cc index e60ab96b9..623bdca1e 100644 --- a/mex/sources/bytecode/BasicSymbolTable.cc +++ b/mex/sources/bytecode/BasicSymbolTable.cc @@ -1,5 +1,5 @@ /* - * Copyright © 2007-2022 Dynare Team + * Copyright © 2007-2023 Dynare Team * * This file is part of Dynare. * @@ -21,12 +21,8 @@ #include "dynmex.h" -BasicSymbolTable::BasicSymbolTable() +BasicSymbolTable::BasicSymbolTable(const mxArray *M_) { - mxArray *M_ = mexGetVariable("global", "M_"); - if (!M_) - mexErrMsgTxt("Can't find global variable M_"); - auto get_field_names = [&](const char *field_name, SymbolType type) { vector r; diff --git a/mex/sources/bytecode/BasicSymbolTable.hh b/mex/sources/bytecode/BasicSymbolTable.hh index fa078c5bd..9a2f5a9a7 100644 --- a/mex/sources/bytecode/BasicSymbolTable.hh +++ b/mex/sources/bytecode/BasicSymbolTable.hh @@ -1,5 +1,5 @@ /* - * Copyright © 2007-2022 Dynare Team + * Copyright © 2007-2023 Dynare Team * * This file is part of Dynare. * @@ -25,6 +25,8 @@ #include #include +#include "dynmex.h" + #include "Bytecode.hh" using namespace std; @@ -32,7 +34,7 @@ using namespace std; class BasicSymbolTable { public: - BasicSymbolTable(); + BasicSymbolTable(const mxArray *M_); string getName(SymbolType type, int tsid) const; pair getIDAndType(const string &name) const; size_t maxEndoNameLength() const; diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 232dd6302..aa7ef1bbd 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -56,13 +56,13 @@ Get_Arguments_and_global_variables(int nrhs, double *params[], double *steady_yd[], size_t &steady_row_y, size_t &steady_col_y, unsigned int &periods, - mxArray *block_structur[], + mxArray **block_structur, bool &steady_state, bool &block_decomposed, bool &evaluate, int &block, - mxArray *M_[], mxArray *options_[], bool &global_temporary_terms, + const mxArray **M_, const mxArray **options_, bool &global_temporary_terms, bool &print, - mxArray *GlobalTemporaryTerms[], - bool *extended_path, mxArray *ep_struct[]) + mxArray **GlobalTemporaryTerms, + bool *extended_path, mxArray **ep_struct) { int count_array_argument {0}; size_t pos; @@ -78,30 +78,36 @@ Get_Arguments_and_global_variables(int nrhs, switch (count_array_argument) { case 0: + *M_ = prhs[i]; + break; + case 1: + *options_ = prhs[i]; + break; + case 2: *yd = mxGetPr(prhs[i]); row_y = mxGetM(prhs[i]); col_y = mxGetN(prhs[i]); break; - case 1: + case 3: *xd = mxGetPr(prhs[i]); row_x = mxGetM(prhs[i]); col_x = mxGetN(prhs[i]); break; - case 2: + case 4: *params = mxGetPr(prhs[i]); break; - case 3: + case 5: *steady_yd = mxGetPr(prhs[i]); steady_row_y = mxGetM(prhs[i]); steady_col_y = mxGetN(prhs[i]); break; - case 4: + case 6: periods = static_cast(mxGetScalar(prhs[i])); break; - case 5: + case 7: *block_structur = mxDuplicateArray(prhs[i]); break; - case 6: + case 8: global_temporary_terms = true; *GlobalTemporaryTerms = mxDuplicateArray(prhs[i]); break; @@ -151,27 +157,25 @@ Get_Arguments_and_global_variables(int nrhs, throw FatalException{"In main, unknown argument : " + Get_Argument(prhs[i])}; } } - if (count_array_argument < 5) + if (steady_state) { - if (count_array_argument == 3 && steady_state) + if (count_array_argument < 5) + throw FatalException{"In a static context, the following arguments have to be indicated: M_, options_, y, x, params"}; + if (count_array_argument < 7) periods = 1; - else - throw FatalException{"In main, missing arguments. All the following arguments have to be indicated y, x, params, ys, periods"}; } - *M_ = mexGetVariable("global", "M_"); - if (!*M_) - throw FatalException{"In main, global variable not found: M_"}; - - *options_ = mexGetVariable("global", "options_"); - if (!*options_) - throw FatalException{"In main, global variable not found: options_"}; + else + { + if (count_array_argument < 7) + throw FatalException{"In a dynamic context, the following arguments have to be indicated: M_, options_, y, x, params, steady_state, periods"}; + } } /* The gateway routine */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - mxArray *M_, *options_; + const mxArray *M_, *options_; mxArray *GlobalTemporaryTerms; mxArray *block_structur = nullptr; size_t i, row_y = 0, col_y = 0, row_x = 0, col_x = 0; @@ -226,7 +230,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mexPrintf("**************************************\n"); #endif - BasicSymbolTable symbol_table; + BasicSymbolTable symbol_table{M_}; vector dates; if (extended_path)