stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Bytecode MEX: get M_ and options_ through input arguments rather than as global variables

kalman-mex
Sébastien Villemot 2023-10-11 16:50:18 -04:00
parent b59dc2cf1a
commit e4a4d2d8e6
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
14 changed files with 71 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/dyn_ramsey_static.m
View File

@ -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

4
matlab/evaluate_static_model.m
View File

@ -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

20
matlab/evaluate_steady_state.m
View File

@ -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

8
matlab/model_diagnostics.m
View File

@ -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

7
matlab/partial_information/add_auxiliary_variables_to_steadystate.m
View File

@ -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 <https://www.gnu.org/licenses/>.
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

12
matlab/perfect-foresight-models/det_cond_forecast.m
View File

@ -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);

2
matlab/perfect-foresight-models/perfect_foresight_solver.m
View File

@ -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;

6
matlab/perfect-foresight-models/perfect_foresight_solver_core.m
View File

@ -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

4
matlab/print_bytecode_dynamic_model.m
View File

@ -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');

4
matlab/print_bytecode_static_model.m
View File

@ -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');

2
matlab/stochastic_solvers.m
View File

@ -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

8
mex/sources/bytecode/BasicSymbolTable.cc
View File

@ -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<string> r;

6
mex/sources/bytecode/BasicSymbolTable.hh
View File

@ -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 <map>
#include <utility>
#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<SymbolType, int> getIDAndType(const string &name) const;
size_t maxEndoNameLength() const;

50
mex/sources/bytecode/bytecode.cc
View File

@ -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<int>(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<string> dates;
if (extended_path)