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

Bytecode: in evaluate+dynamic mode, return residuals with as many columns as periods 

Previously it would also include initial and terminal conditions (i.e.
residuals would have periods+maximum_lag+maximum_lead columns). But we do not
care about residuals at the initial and terminal conditions.

This change is for consistency with the perfect_foresight_problem MEX.
remove-submodule
Sébastien Villemot 2023-06-12 19:12:12 +02:00
parent 8cf8731c65
commit 08d378c244
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
4 changed files with 25 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -228,10 +228,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(0, true, false);
if (single_block)
for (int j = 0; j < size; j++)
residual[it_*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
residual[(it_-y_kmin)*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
else
for (int j = 0; j < size; j++)
residual[it_*size+Block_Contain[j].Equation] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
}
}
break;
@ -258,10 +258,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(0, true, false);
if (!single_block)
for (int j = 0; j < size; j++)
residual[Per_y_+Block_Contain[j].Equation] = r[j];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = r[j];
else
for (int j = 0; j < size; j++)
residual[it_*size+j] = r[j];
residual[(it_-y_kmin)*size+j] = r[j];
}
}
mxFree(g1);
@ -296,10 +296,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(0, true, false);
if (!single_block)
for (int j = 0; j < size; j++)
residual[it_*y_size+Block_Contain[j].Equation] = r[j];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = r[j];
else
for (int j = 0; j < size; j++)
residual[it_*size+j] = r[j];
residual[(it_-y_kmin)*size+j] = r[j];
}
}
mxFree(r);
@ -323,10 +323,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(0, true, false);
if (single_block)
for (int j = 0; j < size; j++)
residual[it_*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
residual[(it_-y_kmin)*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
else
for (int j = 0; j < size; j++)
residual[it_*size+Block_Contain[j].Equation] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
}
}
break;
@ -353,10 +353,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(0, true, false);
if (!single_block)
for (int j = 0; j < size; j++)
residual[Per_y_+Block_Contain[j].Equation] = r[j];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = r[j];
else
for (int j = 0; j < size; j++)
residual[it_*size+j] = r[j];
residual[(it_-y_kmin)*size+j] = r[j];
}
}
mxFree(g1);
@ -388,10 +388,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(0, true, false);
if (!single_block)
for (int j = 0; j < size; j++)
residual[Per_y_+Block_Contain[j].Equation] = r[j];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = r[j];
else
for (int j = 0; j < size; j++)
residual[it_*size+j] = r[j];
residual[(it_-y_kmin)*size+j] = r[j];
}
}
mxFree(r);
@ -413,10 +413,10 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
compute_block_time(Per_u_, true, false);
if (!single_block)
for (int j = 0; j < size; j++)
residual[it_*y_size+Block_Contain[j].Equation] = r[j];
residual[(it_-y_kmin)*y_size+Block_Contain[j].Equation] = r[j];
else
for (int j = 0; j < size; j++)
residual[it_*size+j] = r[j];
residual[(it_-y_kmin)*size+j] = r[j];
}
mxFree(r);
break;
@ -677,7 +677,7 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool
if (steady_state)
residual = vector<double>(y_size);
else
residual = vector<double>(y_size*(periods+y_kmin));
residual = vector<double>(y_size*periods);
}
for (int current_block : blocks)
@ -697,7 +697,7 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool
if (steady_state)
residual = vector<double>(size);
else
residual = vector<double>(size*(periods+y_kmin));
residual = vector<double>(size*periods);
evaluator.printCurrentBlock();
}
else if (evaluate)
@ -721,7 +721,7 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool
if (steady_state)
residual = vector<double>(size);
else
residual = vector<double>(size*(periods+y_kmin));
residual = vector<double>(size*periods);
}
evaluate_a_block(true, block >= 0, bin_basename);
}
@ -737,7 +737,7 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool
jacobian_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockNbColJacob(), mxREAL);
jacobian_exo_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockExoSize(), mxREAL);
jacobian_det_exo_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockExoDetSize(), mxREAL);
residual = vector<double>(size*(periods+y_kmin));
residual = vector<double>(size*periods);
result = simulate_a_block(vector_table_conditional_local, block >= 0, bin_basename);
}
else

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

@ -40,6 +40,7 @@ private:
bool global_temporary_terms;
bool print;
int col_x, col_y;
vector<double> residual;
void evaluate_over_periods(bool forward);
void solve_simple_one_periods();
void solve_simple_over_periods(bool forward);

1
mex/sources/bytecode/SparseMatrix.hh
View File

@ -178,7 +178,6 @@ protected:
double res1a;
long int nop_all, nop1, nop2;
map<tuple<int, int, int>, int> IM_i;
vector<double> residual;
int u_count_alloc, u_count_alloc_save;
vector<double *> jac;
double *jcb;

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

@ -746,40 +746,18 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
bool dont_store_a_structure = false;
if (nlhs > 0)
{
if (block >= 0)
if (evaluate)
{
if (evaluate)
{
vector<double> residual = interprete.get_residual();
plhs[0] = mxCreateDoubleMatrix(static_cast<int>(residual.size()/static_cast<double>(col_y)),
static_cast<int>(col_y), mxREAL);
pind = mxGetPr(plhs[0]);
for (i = 0; i < residual.size(); i++)
pind[i] = residual[i];
}
else
{
int out_periods = extended_path ? max_periods + y_kmin : row_y;
plhs[0] = mxCreateDoubleMatrix(out_periods, static_cast<int>(col_y), mxREAL);
pind = mxGetPr(plhs[0]);
for (i = 0; i < out_periods*col_y; i++)
pind[i] = y[i];
}
vector<double> residual = interprete.get_residual();
plhs[0] = mxCreateDoubleMatrix(static_cast<int>(residual.size()/periods),
static_cast<int>(periods), mxREAL);
std::copy(residual.begin(), residual.end(), mxGetPr(plhs[0]));
}
else
{
int out_periods = extended_path ? max_periods + y_kmin : col_y;
plhs[0] = mxCreateDoubleMatrix(static_cast<int>(row_y), out_periods, mxREAL);
pind = mxGetPr(plhs[0]);
if (evaluate)
{
vector<double> residual = interprete.get_residual();
for (i = 0; i < residual.size(); i++)
pind[i] = residual[i];
}
else
for (i = 0; i < row_y*out_periods; i++)
pind[i] = y[i];
std::copy_n(y, row_y*out_periods, mxGetPr(plhs[0]));
}
if (nlhs > 1)
{