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

Bytecode w/ block decomposition: no longer output derivatives w.r.t. exogenous and endogenous outside the block

silicon
Sébastien Villemot 2023-01-17 13:50:16 +01:00
parent 2bdc185f1d
commit 0c5c96e724
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 12 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
mex/sources/bytecode/Evaluate.cc
View File

@ -1,5 +1,5 @@
/*
* Copyright © 2013-2022 Dynare Team
* Copyright © 2013-2023 Dynare Team
*
* This file is part of Dynare.
*
@ -54,7 +54,6 @@ Evaluate::error_location(it_code_type expr_begin, it_code_type faulty_op, bool s
Error_loc << "equation";
break;
case ExpressionType::FirstEndoDerivative:
case ExpressionType::FirstOtherEndoDerivative:
case ExpressionType::FirstExoDerivative:
case ExpressionType::FirstExodetDerivative:
Error_loc << "first order derivative of equation";
@ -71,9 +70,6 @@ Evaluate::error_location(it_code_type expr_begin, it_code_type faulty_op, bool s
case ExpressionType::FirstEndoDerivative:
Error_loc << " with respect to endogenous variable " << symbol_table.getName(SymbolType::endogenous, EQN_dvar1);
break;
case ExpressionType::FirstOtherEndoDerivative:
Error_loc << " with respect to other endogenous variable " << symbol_table.getName(SymbolType::endogenous, EQN_dvar1);
break;
case ExpressionType::FirstExoDerivative:
Error_loc << " with respect to exogenous variable " << symbol_table.getName(SymbolType::exogenous, EQN_dvar1);
break;
@ -180,9 +176,6 @@ Evaluate::print_expression(const it_code_type &expr_begin, const optional<it_cod
case ExpressionType::FirstEndoDerivative:
equation_type = ExpressionType::FirstEndoDerivative;
break;
case ExpressionType::FirstOtherEndoDerivative:
equation_type = ExpressionType::FirstOtherEndoDerivative;
break;
case ExpressionType::FirstExoDerivative:
equation_type = ExpressionType::FirstExoDerivative;
break;
@ -345,8 +338,6 @@ Evaluate::print_expression(const it_code_type &expr_begin, const optional<it_cod
{
case ExpressionType::FirstEndoDerivative:
return "jacob";
case ExpressionType::FirstOtherEndoDerivative:
return "jacob_other_endo";
case ExpressionType::FirstExoDerivative:
return "jacob_exo";
case ExpressionType::FirstExodetDerivative:
@ -805,7 +796,7 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
bool go_on = true;
double ll;
double rr;
double *jacob = nullptr, *jacob_other_endo = nullptr, *jacob_exo = nullptr, *jacob_exo_det = nullptr;
double *jacob = nullptr, *jacob_exo = nullptr, *jacob_exo_det = nullptr;
EQN_block = block_num;
stack<double> Stack;
ExternalFunctionCallType call_type{ExternalFunctionCallType::levelWithoutDerivative};
@ -819,7 +810,6 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
jacob = mxGetPr(jacobian_block[block_num]);
if (!steady_state)
{
jacob_other_endo = mxGetPr(jacobian_other_endo_block[block_num]);
jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
jacob_exo_det = mxGetPr(jacobian_det_exo_block[block_num]);
}
@ -866,15 +856,6 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
EQN_dvar1 = static_cast<FNUMEXPR_ *>(*it_code)->get_dvariable1();
EQN_lag1 = static_cast<FNUMEXPR_ *>(*it_code)->get_lag1();
break;
case ExpressionType::FirstOtherEndoDerivative:
#ifdef DEBUG
mexPrintf("FirstOtherEndoDerivative\n");
#endif
EQN_type = ExpressionType::FirstOtherEndoDerivative;
EQN_equation = static_cast<FNUMEXPR_ *>(*it_code)->get_equation();
EQN_dvar1 = static_cast<FNUMEXPR_ *>(*it_code)->get_dvariable1();
EQN_lag1 = static_cast<FNUMEXPR_ *>(*it_code)->get_lag1();
break;
case ExpressionType::FirstExoDerivative:
#ifdef DEBUG
mexPrintf("FirstExoDerivative\n");
@ -1309,18 +1290,6 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
#endif
jacob[eq + size*pos_col] = rr;
break;
case ExpressionType::FirstOtherEndoDerivative:
//eq = static_cast<FSTPG3_ *>(*it_code)->get_row();
eq = EQN_equation;
var = static_cast<FSTPG3_ *>(*it_code)->get_col();
lag = static_cast<FSTPG3_ *>(*it_code)->get_lag();
pos_col = static_cast<FSTPG3_ *>(*it_code)->get_col_pos();
#ifdef DEBUG
mexPrintf("other_endo eq=%d, pos_col=%d, size=%d\n", eq, pos_col, size);
mexEvalString("drawnow;");
#endif
jacob_other_endo[eq + size*pos_col] = rr;
break;
case ExpressionType::FirstExoDerivative:
//eq = static_cast<FSTPG3_ *>(*it_code)->get_row();
eq = EQN_equation;

4
mex/sources/bytecode/Evaluate.hh
View File

@ -1,5 +1,5 @@
/*
* Copyright © 2007-2022 Dynare Team
* Copyright © 2007-2023 Dynare Team
*
* This file is part of Dynare.
*
@ -56,7 +56,7 @@ protected:
double *u;
double *steady_y, *steady_x;
double *g1, *r, *res;
vector<mxArray *> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
vector<mxArray *> jacobian_block, jacobian_exo_block, jacobian_det_exo_block;
mxArray *GlobalTemporaryTerms;
it_code_type start_code, end_code;
double pow1(double a, double b);

12
mex/sources/bytecode/Interpreter.cc
View File

@ -652,9 +652,8 @@ Interpreter::MainLoop(const string &bin_basename, const CodeLoad &code, bool eva
mexPrintf("(0) Allocating Jacobian\n");
#endif
jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_exo_jacob(), mxREAL));
jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_det_exo_jacob(), mxREAL));
jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(), mxREAL));
jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_exo_size(), mxREAL));
jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_det_exo_size(), mxREAL));
}
if (block >= 0)
{
@ -677,9 +676,8 @@ Interpreter::MainLoop(const string &bin_basename, const CodeLoad &code, bool eva
//mexPrintf("(1) Allocating Jacobian fb->get_size()=%d fb->get_nb_col_jacob()=%d\n", fb->get_size(), fb->get_nb_col_jacob());
jacobian_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_jacob(), mxREAL));
//mexPrintf("mxGetPr(jacobian_block[block_num])=%x\n",mxGetPr(jacobian_block[0]));
jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_exo_jacob(), mxREAL));
jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_det_exo_jacob(), mxREAL));
jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(), mxREAL));
jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_exo_size(), mxREAL));
jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_det_exo_size(), mxREAL));
residual = vector<double>(fb->get_size()*(periods+y_kmin));
result = simulate_a_block(vector_table_conditional_local);
@ -689,8 +687,6 @@ Interpreter::MainLoop(const string &bin_basename, const CodeLoad &code, bool eva
jacobian_exo_block.pop_back();
mxDestroyArray(jacobian_det_exo_block.back());
jacobian_det_exo_block.pop_back();
mxDestroyArray(jacobian_other_endo_block.back());
jacobian_other_endo_block.pop_back();
}
else
result = simulate_a_block(vector_table_conditional_local);

5
mex/sources/bytecode/Interpreter.hh
View File

@ -73,11 +73,6 @@ public:
{
return jacobian_det_exo_block[block_num];
}
inline mxArray *
get_jacob_other_endo(int block_num) const
{
return jacobian_other_endo_block[block_num];
}
inline vector<double>
get_residual() const
{

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

@ -799,16 +799,15 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (evaluate)
{
int jacob_field_number = 0, jacob_exo_field_number = 0,
jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
jacob_exo_det_field_number = 0;
if (!block_structur)
{
const char *field_names[] = {"g1", "g1_x", "g1_xd", "g1_o"};
const char *field_names[] = {"g1", "g1_x", "g1_xd"};
jacob_field_number = 0;
jacob_exo_field_number = 1;
jacob_exo_det_field_number = 2;
jacob_other_endo_field_number = 3;
mwSize dims[1] = { static_cast<mwSize>(nb_blocks) };
plhs[1] = mxCreateStructArray(1, dims, 4, field_names);
plhs[1] = mxCreateStructArray(1, dims, 3, field_names);
}
else if (!mxIsStruct(block_structur))
{
@ -827,9 +826,6 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
jacob_exo_det_field_number = mxAddField(plhs[1], "g1_xd");
if (jacob_exo_det_field_number == -1)
mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob_exo_det to the structArray");
jacob_other_endo_field_number = mxAddField(plhs[1], "g1_o");
if (jacob_other_endo_field_number == -1)
mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray");
}
if (!dont_store_a_structure)
for (int i = 0; i < nb_blocks; i++)
@ -839,7 +835,6 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mxSetFieldByNumber(plhs[1], i, jacob_exo_field_number, interprete.get_jacob_exo(i));
mxSetFieldByNumber(plhs[1], i, jacob_exo_det_field_number, interprete.get_jacob_exo_det(i));
mxSetFieldByNumber(plhs[1], i, jacob_other_endo_field_number, interprete.get_jacob_other_endo(i));
}
}
}

2
preprocessor

@ -1 +1 @@
Subproject commit bbdbd0807bad5d17c534782b9fec52d2f407ec5c
Subproject commit 2e09df90e70e8e1bae94e19143a312d036e83543