Drop the legacy representation of the block-decomposed model
As a consequence, and as a temporary measure, always output the non-block-decomposed legacy representation. Also drop the block kalman filter output, and drop now useless variables in M_.block_structure.master
Sébastien Villemot
parent
35ac73fad8
commit
93b9ed6957
|
|
@ -179,140 +179,6 @@ DynamicModel::additionalBlockTemporaryTerms(int blk,
|
|||
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
|
||||
}
|
||||
|
||||
int
|
||||
DynamicModel::nzeDeterministicJacobianForBlock(int blk) const
|
||||
{
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
int block_recursive_size = blocks[blk].getRecursiveSize();
|
||||
|
||||
int nze_deterministic = 0;
|
||||
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|
||||
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
|
||||
nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(),
|
||||
[=](const auto &kv) {
|
||||
auto [eq, var, lag] = kv.first;
|
||||
return eq >= block_recursive_size && var >= block_recursive_size;
|
||||
});
|
||||
else if (simulation_type == BlockSimulationType::solveBackwardSimple
|
||||
|| simulation_type == BlockSimulationType::solveForwardSimple
|
||||
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
||||
|| simulation_type == BlockSimulationType::solveForwardComplete)
|
||||
nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(),
|
||||
[=](const auto &kv) {
|
||||
auto [eq, var, lag] = kv.first;
|
||||
return lag == 0 && eq >= block_recursive_size && var >= block_recursive_size;
|
||||
});
|
||||
return nze_deterministic;
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
|
||||
{
|
||||
temporary_terms_t temporary_terms; // Temp terms written so far
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
int block_size = blocks[blk].size;
|
||||
int block_mfs_size = blocks[blk].mfs_size;
|
||||
|
||||
// Number of nonzero derivatives for the various Jacobians
|
||||
int nze_stochastic = blocks_derivatives[blk].size();
|
||||
int nze_deterministic = nzeDeterministicJacobianForBlock(blk);
|
||||
int nze_other_endo = blocks_derivatives_other_endo[blk].size();
|
||||
int nze_exo = blocks_derivatives_exo[blk].size();
|
||||
int nze_exo_det = blocks_derivatives_exo_det[blk].size();
|
||||
|
||||
filesystem::path filename {packageDir(basename) / "+block" / ("dynamic_" + to_string(blk+1) + ".m")};
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
output << "%" << endl
|
||||
<< "% " << filename.string() << " : Computes dynamic version of one block" << endl
|
||||
<< "%" << endl
|
||||
<< "% Warning : this file is generated automatically by Dynare" << endl
|
||||
<< "% from model file (.mod)" << endl << endl
|
||||
<< "%" << endl;
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << "function [y, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode)" << endl;
|
||||
else
|
||||
output << "function [residual, y, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode)" << endl;
|
||||
|
||||
output << " % ////////////////////////////////////////////////////////////////////////" << endl
|
||||
<< " % //" << " Block "sv.substr(static_cast<int>(log10(blk + 1))) << blk+1
|
||||
<< " //" << endl
|
||||
<< " % // Simulation type "
|
||||
<< BlockSim(simulation_type) << " //" << endl
|
||||
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
|
||||
|
||||
if (simulation_type != BlockSimulationType::evaluateForward
|
||||
&& simulation_type != BlockSimulationType::evaluateBackward)
|
||||
output << " residual=zeros(" << block_mfs_size << ",1);" << endl;
|
||||
|
||||
output << " if stochastic_mode" << endl
|
||||
<< " g1_i=zeros(" << nze_stochastic << ",1);" << endl
|
||||
<< " g1_j=zeros(" << nze_stochastic << ",1);" << endl
|
||||
<< " g1_v=zeros(" << nze_stochastic << ",1);" << endl
|
||||
<< " g1_x_i=zeros(" << nze_exo << ",1);" << endl
|
||||
<< " g1_x_j=zeros(" << nze_exo << ",1);" << endl
|
||||
<< " g1_x_v=zeros(" << nze_exo << ",1);" << endl
|
||||
<< " g1_xd_i=zeros(" << nze_exo_det << ",1);" << endl
|
||||
<< " g1_xd_j=zeros(" << nze_exo_det << ",1);" << endl
|
||||
<< " g1_xd_v=zeros(" << nze_exo_det << ",1);" << endl
|
||||
<< " g1_o_i=zeros(" << nze_other_endo << ",1);" << endl
|
||||
<< " g1_o_j=zeros(" << nze_other_endo << ",1);" << endl
|
||||
<< " g1_o_v=zeros(" << nze_other_endo << ",1);" << endl;
|
||||
if (simulation_type != BlockSimulationType::evaluateForward
|
||||
&& simulation_type != BlockSimulationType::evaluateBackward)
|
||||
output << " else" << endl
|
||||
<< " g1_i=zeros(" << nze_deterministic << ",1);" << endl
|
||||
<< " g1_j=zeros(" << nze_deterministic << ",1);" << endl
|
||||
<< " g1_v=zeros(" << nze_deterministic << ",1);" << endl;
|
||||
output << " end" << endl
|
||||
<< endl;
|
||||
|
||||
writeDynamicPerBlockHelper<ExprNodeOutputType::matlabDynamicModel>(blk, output, temporary_terms,
|
||||
nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);
|
||||
|
||||
output << endl
|
||||
<< " if stochastic_mode" << endl
|
||||
<< " g1=sparse(g1_i, g1_j, g1_v, " << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ");" << endl
|
||||
<< " varargout{1}=sparse(g1_x_i, g1_x_j, g1_x_v, " << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ");" << endl
|
||||
<< " varargout{2}=sparse(g1_xd_i, g1_xd_j, g1_xd_v, " << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
|
||||
<< " varargout{3}=sparse(g1_o_i, g1_o_j, g1_o_v, " << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ");" << endl
|
||||
<< " else" << endl;
|
||||
switch (simulation_type)
|
||||
{
|
||||
case BlockSimulationType::evaluateForward:
|
||||
case BlockSimulationType::evaluateBackward:
|
||||
output << " g1=[];" << endl;
|
||||
break;
|
||||
case BlockSimulationType::solveBackwardSimple:
|
||||
case BlockSimulationType::solveForwardSimple:
|
||||
case BlockSimulationType::solveBackwardComplete:
|
||||
case BlockSimulationType::solveForwardComplete:
|
||||
output << " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
|
||||
<< ", " << block_mfs_size << ");" << endl;
|
||||
break;
|
||||
case BlockSimulationType::solveTwoBoundariesSimple:
|
||||
case BlockSimulationType::solveTwoBoundariesComplete:
|
||||
output << " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
|
||||
<< ", " << 3*block_mfs_size << ");" << endl;
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
output << " end" << endl
|
||||
<< "end" << endl;
|
||||
output.close();
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
|
||||
const temporary_terms_t &temporary_terms_union,
|
||||
|
|
@ -346,212 +212,6 @@ DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file,
|
|||
}
|
||||
}
|
||||
|
||||
vector<filesystem::path>
|
||||
DynamicModel::writeDynamicPerBlockCFiles(const string &basename, const string &mexext,
|
||||
const filesystem::path &matlabroot,
|
||||
const filesystem::path &dynareroot) const
|
||||
{
|
||||
temporary_terms_t temporary_terms; // Temp terms written so far
|
||||
const filesystem::path model_src_dir { filesystem::path{basename} / "model" / "src" };
|
||||
vector<filesystem::path> compiled_object_files;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
int block_size = blocks[blk].size;
|
||||
int block_mfs_size = blocks[blk].mfs_size;
|
||||
|
||||
// Number of nonzero derivatives for the various Jacobians
|
||||
int nze_stochastic = blocks_derivatives[blk].size();
|
||||
int nze_deterministic = nzeDeterministicJacobianForBlock(blk);
|
||||
int nze_other_endo = blocks_derivatives_other_endo[blk].size();
|
||||
int nze_exo = blocks_derivatives_exo[blk].size();
|
||||
int nze_exo_det = blocks_derivatives_exo_det[blk].size();
|
||||
|
||||
filesystem::path filename { model_src_dir / ("dynamic_" + to_string(blk+1) + ".c") };
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
output << "/* Block " << blk+1 << endl
|
||||
<< " " << BlockSim(simulation_type) << " */" << endl
|
||||
<< endl
|
||||
<< "#include <math.h>" << endl
|
||||
<< "#include <stdlib.h>" << endl
|
||||
<< "#include <stdbool.h>" << endl
|
||||
<< R"(#include "mex.h")" << endl
|
||||
<< endl;
|
||||
|
||||
// Write function definition if BinaryOpcode::powerDeriv is used
|
||||
writePowerDerivHeader(output);
|
||||
|
||||
output << endl;
|
||||
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
|
||||
else
|
||||
output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
|
||||
output << '{' << endl;
|
||||
|
||||
writeDynamicPerBlockHelper<ExprNodeOutputType::CDynamicModel>(blk, output, temporary_terms,
|
||||
nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);
|
||||
|
||||
output << '}' << endl
|
||||
<< endl;
|
||||
|
||||
ostringstream header;
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)";
|
||||
else
|
||||
header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **residual, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)";
|
||||
output << header.str() << endl
|
||||
<< '{' << endl
|
||||
<< " int nb_row_x = mxGetM(x);" << endl;
|
||||
|
||||
if (simulation_type != BlockSimulationType::evaluateForward
|
||||
&& simulation_type != BlockSimulationType::evaluateBackward)
|
||||
output << " *residual = mxCreateDoubleMatrix(" << block_mfs_size << ",1,mxREAL);" << endl;
|
||||
|
||||
output << " mxArray *g1_i = NULL, *g1_j = NULL, *g1_v = NULL;" << endl
|
||||
<< " mxArray *g1_x_i = NULL, *g1_x_j = NULL, *g1_x_v = NULL;" << endl
|
||||
<< " mxArray *g1_xd_i = NULL, *g1_xd_j = NULL, *g1_xd_v = NULL;" << endl
|
||||
<< " mxArray *g1_o_i = NULL, *g1_o_j = NULL, *g1_o_v = NULL;" << endl
|
||||
<< " if (mxGetScalar(stochastic_mode)) {" << endl
|
||||
<< " g1_i=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
|
||||
<< " g1_j=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
|
||||
<< " g1_v=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
|
||||
<< " g1_x_i=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
|
||||
<< " g1_x_j=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
|
||||
<< " g1_x_v=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
|
||||
<< " g1_xd_i=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
|
||||
<< " g1_xd_j=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
|
||||
<< " g1_xd_v=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
|
||||
<< " g1_o_i=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl
|
||||
<< " g1_o_j=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl
|
||||
<< " g1_o_v=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl;
|
||||
if (simulation_type != BlockSimulationType::evaluateForward
|
||||
&& simulation_type != BlockSimulationType::evaluateBackward)
|
||||
output << " } else {" << endl
|
||||
<< " g1_i=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl
|
||||
<< " g1_j=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl
|
||||
<< " g1_v=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl;
|
||||
output << " }" << endl
|
||||
<< endl;
|
||||
|
||||
// N.B.: In the following, it_ is decreased by 1, to follow C convention
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << " dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_)-1, mxGetScalar(stochastic_mode), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl;
|
||||
else
|
||||
output << " dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_)-1, mxGetScalar(stochastic_mode), mxGetPr(*residual), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl;
|
||||
|
||||
output << endl
|
||||
<< " if (mxGetScalar(stochastic_mode)) {" << endl
|
||||
<< " mxArray *m = mxCreateDoubleScalar(" << block_size << ");" << endl
|
||||
<< " mxArray *n = mxCreateDoubleScalar(" << blocks_jacob_cols_endo[blk].size() << ");" << endl
|
||||
<< " mxArray *plhs[1];" << endl
|
||||
<< " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
|
||||
<< " *g1=plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_i);" << endl
|
||||
<< " mxDestroyArray(g1_j);" << endl
|
||||
<< " mxDestroyArray(g1_v);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< " n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo[blk].size() << ");" << endl
|
||||
<< " mxArray *prhs_x[5] = { g1_x_i, g1_x_j, g1_x_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs_x, "sparse");)" << endl
|
||||
<< " *g1_x=plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_x_i);" << endl
|
||||
<< " mxDestroyArray(g1_x_j);" << endl
|
||||
<< " mxDestroyArray(g1_x_v);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< " n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
|
||||
<< " mxArray *prhs_xd[5] = { g1_xd_i, g1_xd_j, g1_xd_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs_xd, "sparse");)" << endl
|
||||
<< " *g1_xd=plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_xd_i);" << endl
|
||||
<< " mxDestroyArray(g1_xd_j);" << endl
|
||||
<< " mxDestroyArray(g1_xd_v);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< " n = mxCreateDoubleScalar(" << blocks_jacob_cols_other_endo[blk].size() << ");" << endl
|
||||
<< " mxArray *prhs_o[5] = { g1_o_i, g1_o_j, g1_o_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs_o, "sparse");)" << endl
|
||||
<< " *g1_o=plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_o_i);" << endl
|
||||
<< " mxDestroyArray(g1_o_j);" << endl
|
||||
<< " mxDestroyArray(g1_o_v);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< " mxDestroyArray(m);" << endl
|
||||
<< " } else {" << endl;
|
||||
switch (simulation_type)
|
||||
{
|
||||
case BlockSimulationType::evaluateForward:
|
||||
case BlockSimulationType::evaluateBackward:
|
||||
output << " *g1=mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
|
||||
break;
|
||||
case BlockSimulationType::solveBackwardSimple:
|
||||
case BlockSimulationType::solveForwardSimple:
|
||||
case BlockSimulationType::solveBackwardComplete:
|
||||
case BlockSimulationType::solveForwardComplete:
|
||||
output << " mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
|
||||
<< " mxArray *n = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
|
||||
<< " mxArray *plhs[1];" << endl
|
||||
<< " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
|
||||
<< " *g1=plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_i);" << endl
|
||||
<< " mxDestroyArray(g1_j);" << endl
|
||||
<< " mxDestroyArray(g1_v);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< " mxDestroyArray(m);" << endl;
|
||||
break;
|
||||
case BlockSimulationType::solveTwoBoundariesSimple:
|
||||
case BlockSimulationType::solveTwoBoundariesComplete:
|
||||
output << " mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
|
||||
<< " mxArray *n = mxCreateDoubleScalar(" << 3*block_mfs_size << ");" << endl
|
||||
<< " mxArray *plhs[1];" << endl
|
||||
<< " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
|
||||
<< " *g1=plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_i);" << endl
|
||||
<< " mxDestroyArray(g1_j);" << endl
|
||||
<< " mxDestroyArray(g1_v);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< " mxDestroyArray(m);" << endl;
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
output << " *g1_x=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
|
||||
<< " *g1_xd=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
|
||||
<< " *g1_o=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
|
||||
<< " }" << endl
|
||||
<< "}" << endl;
|
||||
output.close();
|
||||
|
||||
// Compile intermediary object under <MODFILE>/model/src/
|
||||
compiled_object_files.emplace_back(compileMEX(model_src_dir, "dynamic_" + to_string(blk+1),
|
||||
mexext, { filename }, matlabroot, dynareroot,
|
||||
false));
|
||||
|
||||
filename = model_src_dir / ("dynamic_" + to_string(blk+1) + ".h");
|
||||
ofstream header_output{filename, ios::out | ios::binary};
|
||||
if (!header_output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
header_output << header.str() << ';' << endl;
|
||||
header_output.close();
|
||||
}
|
||||
return compiled_object_files;
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeDynamicBytecode(const string &basename) const
|
||||
{
|
||||
|
|
@ -756,127 +416,6 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const
|
|||
output << "]";
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeDynamicBlockMFile(const string &basename) const
|
||||
{
|
||||
filesystem::path filename {packageDir(basename) / "dynamic.m"};
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
output << "function [residual, y, T, g1, varargout] = dynamic(nblock, y, x, params, steady_state, T, it_, stochastic_mode)" << endl
|
||||
<< " switch nblock" << endl;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
output << " case " << blk+1 << endl;
|
||||
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << " [y, T, g1, varargout{1:nargout-4}] = " << basename << ".block.dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode);" << endl
|
||||
<< " residual = [];" << endl;
|
||||
else
|
||||
output << " [residual, y, T, g1, varargout{1:nargout-4}] = " << basename << ".block.dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode);" << endl;
|
||||
}
|
||||
output << " end" << endl
|
||||
<< "end" << endl;
|
||||
|
||||
output.close();
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeDynamicBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
|
||||
{
|
||||
const filesystem::path filename {basename + "/model/src/dynamic.c"};
|
||||
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "Error: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
output << "#include <math.h>" << endl
|
||||
<< R"(#include "mex.h")" << endl;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
output << R"(#include "dynamic_)" << blk+1 << R"(.h")" << endl;
|
||||
|
||||
output << endl;
|
||||
writePowerDeriv(output);
|
||||
|
||||
output << endl
|
||||
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
|
||||
<< "{" << endl
|
||||
<< " if (nrhs != 8)" << endl
|
||||
<< R"( mexErrMsgTxt("Requires exactly 8 input arguments");)" << endl
|
||||
<< " if (nlhs > 7)" << endl
|
||||
<< R"( mexErrMsgTxt("Accepts at most 7 output arguments");)" << endl
|
||||
<< " int nblock = (int) mxGetScalar(prhs[0]);" << endl
|
||||
<< " const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *steady_state = prhs[4], *T = prhs[5], *it_ = prhs[6], *stochastic_mode = prhs[7];" << endl
|
||||
<< " mxArray *T_new = mxDuplicateArray(T);" << endl
|
||||
<< " mxArray *y_new = mxDuplicateArray(y);" << endl
|
||||
<< " mxArray *residual, *g1, *g1_x, *g1_xd, *g1_o;" << endl
|
||||
<< " switch (nblock)" << endl
|
||||
<< " {" << endl;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
output << " case " << blk+1 << ':' << endl;
|
||||
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << " dynamic_" << blk+1 << "_mx(y_new, x, params, steady_state, T_new, it_, stochastic_mode, &g1, &g1_x, &g1_xd, &g1_o);" << endl
|
||||
<< " residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
|
||||
else
|
||||
output << " dynamic_" << blk+1 << "_mx(y_new, x, params, steady_state, T_new, it_, stochastic_mode, &residual, &g1, &g1_x, &g1_xd, &g1_o);" << endl;
|
||||
output << " break;" << endl;
|
||||
}
|
||||
output << " }" << endl
|
||||
<< endl
|
||||
<< " if (nlhs >= 1)" << endl
|
||||
<< " plhs[0] = residual;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(residual);" << endl
|
||||
<< " if (nlhs >= 2)" << endl
|
||||
<< " plhs[1] = y_new;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(y_new);" << endl
|
||||
<< " if (nlhs >= 3)" << endl
|
||||
<< " plhs[2] = T_new;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(T_new);" << endl
|
||||
<< " if (nlhs >= 4)" << endl
|
||||
<< " plhs[3] = g1;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(g1);" << endl
|
||||
<< " if (nlhs >= 5)" << endl
|
||||
<< " plhs[4] = g1_x;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(g1_x);" << endl
|
||||
<< " if (nlhs >= 6)" << endl
|
||||
<< " plhs[5] = g1_xd;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(g1_xd);" << endl
|
||||
<< " if (nlhs >= 7)" << endl
|
||||
<< " plhs[6] = g1_o;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(g1_o);" << endl
|
||||
<< "}" << endl;
|
||||
|
||||
output.close();
|
||||
|
||||
per_block_object_files.push_back(filename);
|
||||
compileMEX(packageDir(basename), "dynamic", mexext, per_block_object_files, matlabroot, dynareroot);
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeDynamicMWrapperFunction(const string &basename, const string &ending) const
|
||||
{
|
||||
|
|
@ -1353,8 +892,7 @@ DynamicModel::includeExcludeEquations(const string &inc_exc_option_value, bool e
|
|||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
|
||||
const vector<int> &state_var, bool estimation_present) const
|
||||
DynamicModel::writeBlockDriverOutput(ostream &output) const
|
||||
{
|
||||
output << "M_.block_structure.time_recursive = " << boolalpha << time_recursive_block_decomposition << ";" << endl;
|
||||
|
||||
|
|
@ -1362,14 +900,6 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
|
|||
{
|
||||
int block_size = blocks[blk].size;
|
||||
output << "M_.block_structure.block(" << blk+1 << ").Simulation_Type = " << static_cast<int>(blocks[blk].simulation_type) << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_lag = " << blocks[blk].max_lag << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_lead = " << blocks[blk].max_lead << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_endo_lag = " << blocks[blk].max_endo_lag << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_endo_lead = " << blocks[blk].max_endo_lead << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_exo_lag = " << blocks[blk].max_exo_lag << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_exo_lead = " << blocks[blk].max_exo_lead << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_exo_det_lag = " << blocks[blk].max_exo_det_lag << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").maximum_exo_det_lead = " << blocks[blk].max_exo_det_lead << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").endo_nbr = " << block_size << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").mfs = " << blocks[blk].mfs_size << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").equation = [";
|
||||
|
|
@ -1380,91 +910,6 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
|
|||
for (int var = 0; var < block_size; var++)
|
||||
output << " " << getBlockVariableID(blk, var)+1;
|
||||
output << "];" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").exogenous = [";
|
||||
for (int exo : blocks_exo[blk])
|
||||
output << " " << exo+1;
|
||||
output << "];" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").exo_nbr = " << blocks_exo[blk].size() << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").exogenous_det = [";
|
||||
for (int exo_det : blocks_exo_det[blk])
|
||||
output << " " << exo_det+1;
|
||||
output << "];" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").exo_det_nbr = " << blocks_exo_det[blk].size() << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").other_endogenous = [";
|
||||
for (int other_endo : blocks_other_endo[blk])
|
||||
output << " " << other_endo+1;
|
||||
output << "];" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").other_endogenous_block = [";
|
||||
for (int other_endo : blocks_other_endo[blk])
|
||||
output << " " << endo2block[other_endo]+1;
|
||||
output << "];" << endl;
|
||||
|
||||
output << "M_.block_structure.block(" << blk+1 << ").tm1 = zeros(" << blocks_other_endo[blk].size() << ", " << state_var.size() << ");" << endl;
|
||||
for (int line{1};
|
||||
auto other_endo : blocks_other_endo[blk])
|
||||
{
|
||||
if (auto it = find(state_var.begin(), state_var.end(), other_endo);
|
||||
it != state_var.end())
|
||||
output << "M_.block_structure.block(" << blk+1 << ").tm1("
|
||||
<< line << ", "
|
||||
<< distance(state_var.begin(), it)+1 << ") = 1;" << endl;
|
||||
line++;
|
||||
}
|
||||
|
||||
output << "M_.block_structure.block(" << blk+1 << ").other_endo_nbr = " << blocks_other_endo[blk].size() << ";" << endl;
|
||||
|
||||
int count_lead_lag_incidence = 0;
|
||||
vector<int> local_state_var;
|
||||
output << "M_.block_structure.block(" << blk+1 << ").lead_lag_incidence = [" << endl;
|
||||
for (int lag = -1; lag <= 1; lag++)
|
||||
{
|
||||
for (int var = 0; var < block_size; var++)
|
||||
{
|
||||
for (int eq = 0; eq < block_size; eq++)
|
||||
if (blocks_derivatives[blk].contains({ eq, var, lag }))
|
||||
{
|
||||
if (lag == -1)
|
||||
local_state_var.push_back(getBlockVariableID(blk, var));
|
||||
output << " " << ++count_lead_lag_incidence;
|
||||
goto var_found;
|
||||
}
|
||||
output << " 0";
|
||||
var_found:
|
||||
;
|
||||
}
|
||||
output << ";" << endl;
|
||||
}
|
||||
output << "];" << endl;
|
||||
|
||||
output << "M_.block_structure.block(" << blk+1 << ").sorted_col_dr_ghx = [";
|
||||
for (int lsv : local_state_var)
|
||||
output << distance(state_var.begin(), find(state_var.begin(), state_var.end(), lsv))+1 << " ";
|
||||
output << "];" << endl;
|
||||
|
||||
count_lead_lag_incidence = 0;
|
||||
output << "M_.block_structure.block(" << blk+1 << ").lead_lag_incidence_other = [" << endl;
|
||||
for (int lag = -1; lag <= 1; lag++)
|
||||
{
|
||||
for (int other_endo : blocks_other_endo[blk])
|
||||
{
|
||||
for (int eq = 0; eq < block_size; eq++)
|
||||
if (blocks_derivatives_other_endo[blk].contains({ eq, other_endo, lag }))
|
||||
{
|
||||
output << " " << ++count_lead_lag_incidence;
|
||||
goto other_endo_found;
|
||||
}
|
||||
output << " 0";
|
||||
other_endo_found:
|
||||
;
|
||||
}
|
||||
output << ";" << endl;
|
||||
}
|
||||
output << "];" << endl;
|
||||
|
||||
output << "M_.block_structure.block(" << blk+1 << ").n_static = " << blocks[blk].n_static << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").n_forward = " << blocks[blk].n_forward << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").n_backward = " << blocks[blk].n_backward << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").n_mixed = " << blocks[blk].n_mixed << ";" << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").is_linear = " << boolalpha << blocks[blk].linear << ';' << endl
|
||||
<< "M_.block_structure.block(" << blk+1 << ").NNZDerivatives = " << blocks_derivatives[blk].size() << ';' << endl;
|
||||
}
|
||||
|
|
@ -1499,138 +944,10 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
|
|||
output << "];" << endl;
|
||||
}
|
||||
output << "M_.block_structure.dyn_tmp_nbr = " << blocks_temporary_terms_idxs.size() << ';' << endl;
|
||||
|
||||
if (estimation_present)
|
||||
{
|
||||
filesystem::create_directories(basename + "/model/bytecode");
|
||||
const filesystem::path main_name {basename + "/model/bytecode/kfi"};
|
||||
ofstream KF_index_file{main_name, ios::out | ios::binary | ios::ate};
|
||||
if (!KF_index_file.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << main_name.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
int n_obs = symbol_table.observedVariablesNbr();
|
||||
int n_state = state_var.size();
|
||||
for (int it : state_var)
|
||||
if (symbol_table.isObservedVariable(symbol_table.getID(SymbolType::endogenous, it)))
|
||||
n_obs--;
|
||||
|
||||
int n = n_obs + n_state;
|
||||
output << "M_.nobs_non_statevar = " << n_obs << ";" << endl;
|
||||
int nb_diag = 0;
|
||||
|
||||
vector<int> i_nz_state_var(n);
|
||||
for (int i = 0; i < n_obs; i++)
|
||||
i_nz_state_var[i] = n;
|
||||
int lp = n_obs;
|
||||
|
||||
vector<int> state_equ;
|
||||
for (int it : state_var)
|
||||
state_equ.push_back(eq_idx_block2orig[endo_idx_orig2block[it]]);
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
int nze = 0;
|
||||
for (int i = 0; i < blocks[blk].size; i++)
|
||||
if (int var = getBlockVariableID(blk, i);
|
||||
find(state_var.begin(), state_var.end(), var) != state_var.end())
|
||||
nze++;
|
||||
|
||||
if (blk == 0)
|
||||
{
|
||||
set<pair<int, int>> row_state_var_incidence;
|
||||
for (const auto &[idx, ignore] : blocks_derivatives[blk])
|
||||
if (auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(blk, get<1>(idx)));
|
||||
it_state_var != state_var.end())
|
||||
if (auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(blk, get<0>(idx)));
|
||||
it_state_equ != state_equ.end())
|
||||
row_state_var_incidence.emplace(it_state_equ - state_equ.begin(), it_state_var - state_var.begin());
|
||||
auto row_state_var_incidence_it = row_state_var_incidence.begin();
|
||||
bool diag = true;
|
||||
int nb_diag_r = 0;
|
||||
while (row_state_var_incidence_it != row_state_var_incidence.end() && diag)
|
||||
{
|
||||
diag = (row_state_var_incidence_it->first == row_state_var_incidence_it->second);
|
||||
if (diag)
|
||||
{
|
||||
int equ = row_state_var_incidence_it->first;
|
||||
row_state_var_incidence_it++;
|
||||
if (equ != row_state_var_incidence_it->first)
|
||||
nb_diag_r++;
|
||||
}
|
||||
|
||||
}
|
||||
set<pair<int, int>> col_state_var_incidence;
|
||||
for (auto [equ, var] : row_state_var_incidence)
|
||||
col_state_var_incidence.emplace(var, equ);
|
||||
auto col_state_var_incidence_it = col_state_var_incidence.begin();
|
||||
diag = true;
|
||||
int nb_diag_c = 0;
|
||||
while (col_state_var_incidence_it != col_state_var_incidence.end() && diag)
|
||||
{
|
||||
diag = (col_state_var_incidence_it->first == col_state_var_incidence_it->second);
|
||||
if (diag)
|
||||
{
|
||||
int var = col_state_var_incidence_it->first;
|
||||
col_state_var_incidence_it++;
|
||||
if (var != col_state_var_incidence_it->first)
|
||||
nb_diag_c++;
|
||||
}
|
||||
}
|
||||
nb_diag = min(nb_diag_r, nb_diag_c);
|
||||
row_state_var_incidence.clear();
|
||||
col_state_var_incidence.clear();
|
||||
}
|
||||
for (int i = 0; i < nze; i++)
|
||||
i_nz_state_var[lp + i] = lp + nze;
|
||||
lp += nze;
|
||||
}
|
||||
output << "M_.nz_state_var = [";
|
||||
for (int i = 0; i < lp; i++)
|
||||
output << i_nz_state_var[i] << " ";
|
||||
output << "];" << endl
|
||||
<< "M_.n_diag = " << nb_diag << ";" << endl;
|
||||
KF_index_file.write(reinterpret_cast<char *>(&nb_diag), sizeof(nb_diag));
|
||||
|
||||
using index_KF = pair<int, pair<int, int >>;
|
||||
vector<index_KF> v_index_KF;
|
||||
for (int i = 0; i < n; i++)
|
||||
for (int j = n_obs; j < n; j++)
|
||||
{
|
||||
int j1 = j - n_obs;
|
||||
int j1_n_state = j1 * n_state - n_obs;
|
||||
if ((i < n_obs) || (i >= nb_diag + n_obs) || (j1 >= nb_diag))
|
||||
for (int k = n_obs; k < i_nz_state_var[i]; k++)
|
||||
v_index_KF.emplace_back(i + j1 * n, pair(i + k * n, k + j1_n_state));
|
||||
}
|
||||
int size_v_index_KF = v_index_KF.size();
|
||||
|
||||
KF_index_file.write(reinterpret_cast<char *>(&size_v_index_KF), sizeof(size_v_index_KF));
|
||||
for (auto &it : v_index_KF)
|
||||
KF_index_file.write(reinterpret_cast<char *>(&it), sizeof(index_KF));
|
||||
|
||||
vector<index_KF> v_index_KF_2;
|
||||
int n_n_obs = n * n_obs;
|
||||
for (int i = 0; i < n; i++)
|
||||
for (int j = i; j < n; j++)
|
||||
if ((i < n_obs) || (i >= nb_diag + n_obs) || (j < n_obs) || (j >= nb_diag + n_obs))
|
||||
for (int k = n_obs; k < i_nz_state_var[j]; k++)
|
||||
{
|
||||
int k_n = k * n;
|
||||
v_index_KF_2.emplace_back(i * n + j, pair(i + k_n - n_n_obs, j + k_n));
|
||||
}
|
||||
int size_v_index_KF_2 = v_index_KF_2.size();
|
||||
|
||||
KF_index_file.write(reinterpret_cast<char *>(&size_v_index_KF_2), sizeof(size_v_index_KF_2));
|
||||
for (auto &it : v_index_KF_2)
|
||||
KF_index_file.write(reinterpret_cast<char *>(&it), sizeof(index_KF));
|
||||
KF_index_file.close();
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool estimation_present, bool compute_xrefs) const
|
||||
DynamicModel::writeDriverOutput(ostream &output, bool compute_xrefs) const
|
||||
{
|
||||
/* Writing initialisation for M_.lead_lag_incidence matrix
|
||||
M_.lead_lag_incidence is a matrix with as many columns as there are
|
||||
|
|
@ -1760,7 +1077,7 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool es
|
|||
|
||||
// Write the block structure of the model
|
||||
if (block_decomposed)
|
||||
writeBlockDriverOutput(output, basename, state_var, estimation_present);
|
||||
writeBlockDriverOutput(output);
|
||||
|
||||
output << "M_.state_var = [";
|
||||
for (int it : state_var)
|
||||
|
|
@ -3321,32 +2638,12 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
|
|||
create_directories(model_dir / "bytecode" / "block");
|
||||
|
||||
// Legacy representation
|
||||
if (block)
|
||||
{
|
||||
if (use_dll)
|
||||
{
|
||||
auto per_block_object_files { writeDynamicPerBlockCFiles(basename, mexext, matlabroot, dynareroot) };
|
||||
writeDynamicBlockCFile(basename, move(per_block_object_files), mexext, matlabroot, dynareroot);
|
||||
}
|
||||
else if (julia)
|
||||
{
|
||||
cerr << "'block' option is not available with Julia" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
else // M-files
|
||||
{
|
||||
writeDynamicPerBlockMFiles(basename);
|
||||
writeDynamicBlockMFile(basename);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (use_dll)
|
||||
writeModelCFile<true>(basename, mexext, matlabroot, dynareroot);
|
||||
else if (!julia) // M-files
|
||||
writeDynamicMFile(basename);
|
||||
// The legacy representation is no longer produced for Julia
|
||||
}
|
||||
if (use_dll)
|
||||
writeModelCFile<true>(basename, mexext, matlabroot, dynareroot);
|
||||
else if (!julia) // M-files
|
||||
writeDynamicMFile(basename);
|
||||
// The legacy representation is no longer produced for Julia
|
||||
|
||||
writeDynamicBytecode(basename);
|
||||
if (block_decomposed)
|
||||
writeDynamicBlockBytecode(basename);
|
||||
|
|
|
|||
|
|
@ -120,26 +120,6 @@ private:
|
|||
|
||||
// Writes dynamic model file (MATLAB/Octave version, legacy representation)
|
||||
void writeDynamicMFile(const string &basename) const;
|
||||
/* Writes the main dynamic function of block decomposed model (MATLAB/Octave
|
||||
version, legacy representation) */
|
||||
void writeDynamicBlockMFile(const string &basename) const;
|
||||
/* Writes the main dynamic functions of block decomposed model (C version,
|
||||
legacy representation), then compiles it with the per-block functions into
|
||||
a single MEX */
|
||||
void writeDynamicBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
|
||||
/* Computes the number of nonzero elements in deterministic Jacobian of
|
||||
block-decomposed model */
|
||||
int nzeDeterministicJacobianForBlock(int blk) const;
|
||||
// Helper for writing the per-block dynamic files of block decomposed models (legacy representation)
|
||||
template<ExprNodeOutputType output_type>
|
||||
void writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const;
|
||||
/* Writes the per-block dynamic files of block decomposed model (MATLAB/Octave
|
||||
version, legacy representation) */
|
||||
void writeDynamicPerBlockMFiles(const string &basename) const;
|
||||
/* Writes the per-block dynamic files of block decomposed model (C version,
|
||||
legacy representation).
|
||||
Returns the list of paths to the generated C source files (not the headers) */
|
||||
vector<filesystem::path> writeDynamicPerBlockCFiles(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
|
||||
//! Writes the code of the block-decomposed model in virtual machine bytecode
|
||||
void writeDynamicBlockBytecode(const string &basename) const;
|
||||
// Writes derivatives w.r.t. exo, exo det and other endogenous
|
||||
|
|
@ -153,8 +133,7 @@ private:
|
|||
void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
|
||||
|
||||
// Write the block structure of the model in the driver file
|
||||
void writeBlockDriverOutput(ostream &output, const string &basename,
|
||||
const vector<int> &state_var, bool estimation_present) const;
|
||||
void writeBlockDriverOutput(ostream &output) const;
|
||||
|
||||
// Used by determineBlockDerivativesType()
|
||||
enum class BlockDerivativeType
|
||||
|
|
@ -334,7 +313,7 @@ public:
|
|||
void computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context,
|
||||
bool no_tmp_terms, bool block, bool use_dll);
|
||||
//! Writes information about the dynamic model to the driver file
|
||||
void writeDriverOutput(ostream &output, const string &basename, bool estimation_present, bool compute_xrefs) const;
|
||||
void writeDriverOutput(ostream &output, bool compute_xrefs) const;
|
||||
|
||||
//! Write JSON AST
|
||||
void writeJsonAST(ostream &output) const;
|
||||
|
|
@ -697,171 +676,6 @@ public:
|
|||
set<int> findPacExpectationEquationNumbers() const;
|
||||
};
|
||||
|
||||
template<ExprNodeOutputType output_type>
|
||||
void
|
||||
DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms,
|
||||
int nze_stochastic, int nze_deterministic, int nze_exo,
|
||||
int nze_exo_det, int nze_other_endo) const
|
||||
{
|
||||
static_assert(!isSparseModelOutput(output_type));
|
||||
|
||||
BlockSimulationType simulation_type { blocks[blk].simulation_type };
|
||||
int block_mfs_size { blocks[blk].mfs_size };
|
||||
int block_recursive_size { blocks[blk].getRecursiveSize() };
|
||||
|
||||
// Write residuals and temporary terms (incl. for derivatives)
|
||||
writePerBlockHelper<output_type>(blk, output, temporary_terms);
|
||||
|
||||
if constexpr(isCOutput(output_type))
|
||||
output << " if (stochastic_mode) {" << endl;
|
||||
else
|
||||
output << " if stochastic_mode" << endl;
|
||||
|
||||
ostringstream i_output, j_output, v_output;
|
||||
int line_counter { ARRAY_SUBSCRIPT_OFFSET(output_type) };
|
||||
for (const auto &[indices, d] : blocks_derivatives[blk])
|
||||
{
|
||||
const auto &[eq, var, lag] {indices};
|
||||
int jacob_col { blocks_jacob_cols_endo[blk].at({ var, lag }) };
|
||||
i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
|
||||
j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
|
||||
v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
assert(line_counter == nze_stochastic+ARRAY_SUBSCRIPT_OFFSET(output_type));
|
||||
output << i_output.str() << j_output.str() << v_output.str();
|
||||
|
||||
i_output.str("");
|
||||
j_output.str("");
|
||||
v_output.str("");
|
||||
line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
|
||||
for (const auto &[indices, d] : blocks_derivatives_exo[blk])
|
||||
{
|
||||
const auto &[eq, var, lag] {indices};
|
||||
int jacob_col { blocks_jacob_cols_exo[blk].at({ var, lag }) };
|
||||
i_output << " g1_x_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
|
||||
j_output << " g1_x_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
|
||||
v_output << " g1_x_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
assert(line_counter == nze_exo+ARRAY_SUBSCRIPT_OFFSET(output_type));
|
||||
output << i_output.str() << j_output.str() << v_output.str();
|
||||
|
||||
i_output.str("");
|
||||
j_output.str("");
|
||||
v_output.str("");
|
||||
line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
|
||||
for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
|
||||
{
|
||||
const auto &[eq, var, lag] {indices};
|
||||
int jacob_col { blocks_jacob_cols_exo_det[blk].at({ var, lag }) };
|
||||
i_output << " g1_xd_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
|
||||
j_output << " g1_xd_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
|
||||
v_output << " g1_xd_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
assert(line_counter == nze_exo_det+ARRAY_SUBSCRIPT_OFFSET(output_type));
|
||||
output << i_output.str() << j_output.str() << v_output.str();
|
||||
|
||||
i_output.str("");
|
||||
j_output.str("");
|
||||
v_output.str("");
|
||||
line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
|
||||
for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
|
||||
{
|
||||
const auto &[eq, var, lag] {indices};
|
||||
int jacob_col { blocks_jacob_cols_other_endo[blk].at({ var, lag }) };
|
||||
i_output << " g1_o_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
|
||||
j_output << " g1_o_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
|
||||
v_output << " g1_o_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
assert(line_counter == nze_other_endo+ARRAY_SUBSCRIPT_OFFSET(output_type));
|
||||
output << i_output.str() << j_output.str() << v_output.str();
|
||||
|
||||
// Deterministic mode
|
||||
if (simulation_type != BlockSimulationType::evaluateForward
|
||||
&& simulation_type != BlockSimulationType::evaluateBackward)
|
||||
{
|
||||
if constexpr(isCOutput(output_type))
|
||||
output << " } else {" << endl;
|
||||
else
|
||||
output << " else" << endl;
|
||||
i_output.str("");
|
||||
j_output.str("");
|
||||
v_output.str("");
|
||||
line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
|
||||
if (simulation_type == BlockSimulationType::solveBackwardSimple
|
||||
|| simulation_type == BlockSimulationType::solveForwardSimple
|
||||
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
||||
|| simulation_type == BlockSimulationType::solveForwardComplete)
|
||||
for (const auto &[indices, d] : blocks_derivatives[blk])
|
||||
{
|
||||
const auto &[eq, var, lag] {indices};
|
||||
if (lag == 0 && eq >= block_recursive_size && var >= block_recursive_size)
|
||||
{
|
||||
i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
|
||||
<< eq+1-block_recursive_size << ';' << endl;
|
||||
j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
|
||||
<< var+1-block_recursive_size << ';' << endl;
|
||||
v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
}
|
||||
else // solveTwoBoundariesSimple || solveTwoBoundariesComplete
|
||||
for (const auto &[indices, d] : blocks_derivatives[blk])
|
||||
{
|
||||
const auto &[eq, var, lag] {indices};
|
||||
assert(lag >= -1 && lag <= 1);
|
||||
if (eq >= block_recursive_size && var >= block_recursive_size)
|
||||
{
|
||||
i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
|
||||
<< eq+1-block_recursive_size << ';' << endl;
|
||||
j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
|
||||
<< var+1-block_recursive_size+block_mfs_size*(lag+1) << ';' << endl;
|
||||
v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
}
|
||||
assert(line_counter == nze_deterministic+ARRAY_SUBSCRIPT_OFFSET(output_type));
|
||||
output << i_output.str() << j_output.str() << v_output.str();
|
||||
}
|
||||
if constexpr(isCOutput(output_type))
|
||||
output << " }" << endl;
|
||||
else
|
||||
output << " end" << endl;
|
||||
}
|
||||
|
||||
template<bool julia>
|
||||
void
|
||||
DynamicModel::writeParamsDerivativesFile(const string &basename) const
|
||||
|
|
|
|||
|
|
@ -367,16 +367,6 @@ ModFile::checkPass(bool nostrict, bool stochastic)
|
|||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
// See dynare#1726
|
||||
if ((stochastic_statement_present || mod_file_struct.check_present) && block && dynamic_model.mfs > 0)
|
||||
{
|
||||
/* NB: If mfs>0 but “block” is not passed, the block-DR routines will not
|
||||
be called, so do not fail in that case (we may want to use the sparse
|
||||
block representation) */
|
||||
cerr << "ERROR: the `block` option used in conjunction with `mfs > 0` is incompatible with check, stoch_simul, estimation, osr, ramsey_policy, discretionary_policy, calib_smoother, identification, methods_of_moments and sensitivity commands" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
|
|
@ -946,7 +936,7 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
|
|||
|
||||
if (dynamic_model.equation_number() > 0)
|
||||
{
|
||||
dynamic_model.writeDriverOutput(mOutputFile, basename, mod_file_struct.estimation_present, compute_xrefs);
|
||||
dynamic_model.writeDriverOutput(mOutputFile, compute_xrefs);
|
||||
if (!no_static)
|
||||
static_model.writeDriverOutput(mOutputFile);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -327,7 +327,6 @@ protected:
|
|||
void writeModelCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
|
||||
|
||||
// Writes per-block residuals and temporary terms (incl. for derivatives)
|
||||
// This part is common to the legacy and sparse representations
|
||||
template<ExprNodeOutputType output_type>
|
||||
void writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
|
||||
|
||||
|
|
|
|||
|
|
@ -95,160 +95,6 @@ StaticModel::StaticModel(const DynamicModel &m) :
|
|||
user_set_compiler = m.user_set_compiler;
|
||||
}
|
||||
|
||||
void
|
||||
StaticModel::writeStaticPerBlockMFiles(const string &basename) const
|
||||
{
|
||||
temporary_terms_t temporary_terms; // Temp terms written so far
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
|
||||
filesystem::path filename {packageDir(basename) / "+block" / ("static_" + to_string(blk+1) + ".m")};
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
output << "%" << endl
|
||||
<< "% " << filename.string() << " : Computes static version of one block" << endl
|
||||
<< "%" << endl
|
||||
<< "% Warning : this file is generated automatically by Dynare" << endl
|
||||
<< "% from model file (.mod)" << endl << endl
|
||||
<< "%" << endl;
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << "function [y, T] = static_" << blk+1 << "(y, x, params, T)" << endl;
|
||||
else
|
||||
output << "function [residual, y, T, g1] = static_" << blk+1 << "(y, x, params, T)" << endl;
|
||||
|
||||
output << " % ////////////////////////////////////////////////////////////////////////" << endl
|
||||
<< " % //" << " Block "sv.substr(static_cast<int>(log10(blk + 1))) << blk+1
|
||||
<< " //" << endl
|
||||
<< " % // Simulation type "
|
||||
<< BlockSim(simulation_type) << " //" << endl
|
||||
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
|
||||
|
||||
if (simulation_type != BlockSimulationType::evaluateBackward
|
||||
&& simulation_type != BlockSimulationType::evaluateForward)
|
||||
output << " residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl
|
||||
<< " g1_i=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
|
||||
<< " g1_j=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
|
||||
<< " g1_v=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
|
||||
<< endl;
|
||||
|
||||
writeStaticPerBlockHelper<ExprNodeOutputType::matlabStaticModel>(blk, output, temporary_terms);
|
||||
|
||||
if (simulation_type != BlockSimulationType::evaluateBackward
|
||||
&& simulation_type != BlockSimulationType::evaluateForward)
|
||||
output << endl
|
||||
<< " g1=sparse(g1_i, g1_j, g1_v, " << blocks[blk].mfs_size << "," << blocks[blk].mfs_size << ");" << endl;
|
||||
|
||||
output << "end" << endl;
|
||||
output.close();
|
||||
}
|
||||
}
|
||||
|
||||
vector<filesystem::path>
|
||||
StaticModel::writeStaticPerBlockCFiles(const string &basename, const string &mexext,
|
||||
const filesystem::path &matlabroot,
|
||||
const filesystem::path &dynareroot) const
|
||||
{
|
||||
temporary_terms_t temporary_terms; // Temp terms written so far
|
||||
const filesystem::path model_src_dir { filesystem::path{basename} / "model" / "src" };
|
||||
vector<filesystem::path> compiled_object_files;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
|
||||
filesystem::path filename { model_src_dir / ("static_" + to_string(blk+1) + ".c") };
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
output << "/* Block " << blk+1 << endl
|
||||
<< " " << BlockSim(simulation_type) << " */" << endl
|
||||
<< endl
|
||||
<< "#include <math.h>" << endl
|
||||
<< "#include <stdlib.h>" << endl
|
||||
<< R"(#include "mex.h")" << endl
|
||||
<< endl;
|
||||
|
||||
// Write function definition if BinaryOpcode::powerDeriv is used
|
||||
writePowerDerivHeader(output);
|
||||
|
||||
output << endl;
|
||||
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl;
|
||||
else
|
||||
output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v)" << endl;
|
||||
output << '{' << endl;
|
||||
|
||||
writeStaticPerBlockHelper<ExprNodeOutputType::CStaticModel>(blk, output, temporary_terms);
|
||||
|
||||
output << '}' << endl
|
||||
<< endl;
|
||||
|
||||
ostringstream header;
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
{
|
||||
header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T)";
|
||||
output << header.str() << endl
|
||||
<< '{' << endl
|
||||
<< " static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T));" << endl
|
||||
<< '}' << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T, mxArray **residual, mxArray **g1)";
|
||||
output << header.str() << endl
|
||||
<< '{' << endl
|
||||
<< " *residual = mxCreateDoubleMatrix(" << blocks[blk].mfs_size << ",1,mxREAL);" << endl
|
||||
<< " mxArray *g1_i = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
|
||||
<< " mxArray *g1_j = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
|
||||
<< " mxArray *g1_v = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
|
||||
<< " static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T), mxGetPr(*residual), mxGetPr(g1_i), mxGetPr(g1_j), mxGetPr(g1_v));" << endl
|
||||
<< " mxArray *plhs[1];" << endl
|
||||
<< " mxArray *m = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl
|
||||
<< " mxArray *n = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl
|
||||
<< " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
|
||||
<< R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
|
||||
<< " *g1 = plhs[0];" << endl
|
||||
<< " mxDestroyArray(g1_i);" << endl
|
||||
<< " mxDestroyArray(g1_j);" << endl
|
||||
<< " mxDestroyArray(g1_v);" << endl
|
||||
<< " mxDestroyArray(m);" << endl
|
||||
<< " mxDestroyArray(n);" << endl
|
||||
<< '}' << endl;
|
||||
}
|
||||
|
||||
output.close();
|
||||
|
||||
// Compile intermediary object under <MODFILE>/model/src/
|
||||
compiled_object_files.emplace_back(compileMEX(model_src_dir, "static_" + to_string(blk+1),
|
||||
mexext, { filename }, matlabroot, dynareroot,
|
||||
false));
|
||||
|
||||
filename = model_src_dir / ("static_" + to_string(blk+1) + ".h");
|
||||
ofstream header_output{filename, ios::out | ios::binary};
|
||||
if (!header_output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
header_output << header.str() << ';' << endl;
|
||||
header_output.close();
|
||||
}
|
||||
return compiled_object_files;
|
||||
}
|
||||
|
||||
void
|
||||
StaticModel::writeStaticBytecode(const string &basename) const
|
||||
{
|
||||
|
|
@ -645,32 +491,12 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
|
|||
create_directories(model_dir / "bytecode" / "block");
|
||||
|
||||
// Legacy representation
|
||||
if (block)
|
||||
{
|
||||
if (use_dll)
|
||||
{
|
||||
auto per_block_object_files { writeStaticPerBlockCFiles(basename, mexext, matlabroot, dynareroot) };
|
||||
writeStaticBlockCFile(basename, move(per_block_object_files), mexext, matlabroot, dynareroot);
|
||||
}
|
||||
else if (julia)
|
||||
{
|
||||
cerr << "'block' option is not available with Julia" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
else // M-files
|
||||
{
|
||||
writeStaticPerBlockMFiles(basename);
|
||||
writeStaticBlockMFile(basename);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (use_dll)
|
||||
writeModelCFile<false>(basename, mexext, matlabroot, dynareroot);
|
||||
else if (!julia) // M-files
|
||||
writeStaticMFile(basename);
|
||||
// The legacy representation is no longer produced for Julia
|
||||
}
|
||||
if (use_dll)
|
||||
writeModelCFile<false>(basename, mexext, matlabroot, dynareroot);
|
||||
else if (!julia) // M-files
|
||||
writeStaticMFile(basename);
|
||||
// The legacy representation is no longer produced for Julia
|
||||
|
||||
writeStaticBytecode(basename);
|
||||
if (block_decomposed)
|
||||
writeStaticBlockBytecode(basename);
|
||||
|
|
@ -699,117 +525,6 @@ StaticModel::exoPresentInEqs() const
|
|||
return false;
|
||||
}
|
||||
|
||||
void
|
||||
StaticModel::writeStaticBlockMFile(const string &basename) const
|
||||
{
|
||||
filesystem::path filename {packageDir(basename) / "static.m"};
|
||||
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
output << "function [residual, y, T, g1] = static(nblock, y, x, params, T)" << endl
|
||||
<< " switch nblock" << endl;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
output << " case " << blk+1 << endl;
|
||||
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << " [y, T] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl
|
||||
<< " residual = [];" << endl
|
||||
<< " g1 = [];" << endl;
|
||||
else
|
||||
output << " [residual, y, T, g1] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl;
|
||||
|
||||
}
|
||||
output << " end" << endl
|
||||
<< "end" << endl;
|
||||
output.close();
|
||||
}
|
||||
|
||||
void
|
||||
StaticModel::writeStaticBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
|
||||
{
|
||||
const filesystem::path filename {basename + "/model/src/static.c"};
|
||||
|
||||
ofstream output{filename, ios::out | ios::binary};
|
||||
if (!output.is_open())
|
||||
{
|
||||
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
output << "#include <math.h>" << endl
|
||||
<< R"(#include "mex.h")" << endl;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
output << R"(#include "static_)" << blk+1 << R"(.h")" << endl;
|
||||
|
||||
output << endl;
|
||||
writePowerDeriv(output);
|
||||
|
||||
output << endl
|
||||
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
|
||||
<< "{" << endl
|
||||
<< " if (nrhs != 5)" << endl
|
||||
<< R"( mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl
|
||||
<< " if (nlhs > 4)" << endl
|
||||
<< R"( mexErrMsgTxt("Accepts at most 4 output arguments");)" << endl
|
||||
<< " int nblock = (int) mxGetScalar(prhs[0]);" << endl
|
||||
<< " const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *T = prhs[4];" << endl
|
||||
<< " mxArray *T_new = mxDuplicateArray(T);" << endl
|
||||
<< " mxArray *y_new = mxDuplicateArray(y);" << endl
|
||||
<< " mxArray *residual, *g1;" << endl
|
||||
<< " switch (nblock)" << endl
|
||||
<< " {" << endl;
|
||||
|
||||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||||
{
|
||||
output << " case " << blk+1 << ':' << endl;
|
||||
|
||||
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
||||
|
||||
if (simulation_type == BlockSimulationType::evaluateBackward
|
||||
|| simulation_type == BlockSimulationType::evaluateForward)
|
||||
output << " static_" << blk+1 << "_mx(y_new, x, params, T_new);" << endl
|
||||
<< " residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl
|
||||
<< " g1 = mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
|
||||
else
|
||||
output << " static_" << blk+1 << "_mx(y_new, x, params, T_new, &residual, &g1);" << endl;
|
||||
output << " break;" << endl;
|
||||
}
|
||||
output << " }" << endl
|
||||
<< endl
|
||||
<< " if (nlhs >= 1)" << endl
|
||||
<< " plhs[0] = residual;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(residual);" << endl
|
||||
<< " if (nlhs >= 2)" << endl
|
||||
<< " plhs[1] = y_new;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(y_new);" << endl
|
||||
<< " if (nlhs >= 3)" << endl
|
||||
<< " plhs[2] = T_new;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(T_new);" << endl
|
||||
<< " if (nlhs >= 4)" << endl
|
||||
<< " plhs[3] = g1;" << endl
|
||||
<< " else" << endl
|
||||
<< " mxDestroyArray(g1);" << endl
|
||||
<< "}" << endl;
|
||||
output.close();
|
||||
|
||||
per_block_object_files.push_back(filename);
|
||||
compileMEX(packageDir(basename), "static", mexext, per_block_object_files, matlabroot, dynareroot);
|
||||
}
|
||||
|
||||
void
|
||||
StaticModel::writeDriverOutput(ostream &output) const
|
||||
{
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
/*
|
||||
* Copyright © 2003-2022 Dynare Team
|
||||
* Copyright © 2003-2023 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
|
|
@ -37,28 +37,6 @@ private:
|
|||
// Writes static model file (MATLAB/Octave version, legacy representation)
|
||||
void writeStaticMFile(const string &basename) const;
|
||||
|
||||
/* Writes the main static function of block decomposed model (MATLAB/Octave
|
||||
version, legacy representation) */
|
||||
void writeStaticBlockMFile(const string &basename) const;
|
||||
|
||||
/* Writes the main static functions of block decomposed model (C version,
|
||||
legacy representation), then compiles it with the per-block functions into
|
||||
a single MEX */
|
||||
void writeStaticBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
|
||||
|
||||
// Helper for writing a per-block static file of block decomposed model (legacy representation)
|
||||
template<ExprNodeOutputType output_type>
|
||||
void writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
|
||||
|
||||
/* Writes the per-block static files of block decomposed model (MATLAB/Octave
|
||||
version, legacy representation) */
|
||||
void writeStaticPerBlockMFiles(const string &basename) const;
|
||||
|
||||
/* Writes the per-block static files of block decomposed model (C version,
|
||||
legacy representation).
|
||||
Returns the list of paths to the generated C source files (not the headers) */
|
||||
vector<filesystem::path> writeStaticPerBlockCFiles(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
|
||||
|
||||
//! Writes the code of the block-decomposed model in virtual machine bytecode
|
||||
void writeStaticBlockBytecode(const string &basename) const;
|
||||
|
||||
|
|
@ -181,43 +159,6 @@ public:
|
|||
void addAllParamDerivId(set<int> &deriv_id_set) override;
|
||||
};
|
||||
|
||||
template<ExprNodeOutputType output_type>
|
||||
void
|
||||
StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
|
||||
{
|
||||
static_assert(!isSparseModelOutput(output_type));
|
||||
|
||||
BlockSimulationType simulation_type { blocks[blk].simulation_type };
|
||||
int block_recursive_size { blocks[blk].getRecursiveSize() };
|
||||
|
||||
// Write residuals and temporary terms (incl. for derivatives)
|
||||
writePerBlockHelper<output_type>(blk, output, temporary_terms);
|
||||
|
||||
// The Jacobian if we have to solve the block
|
||||
if (simulation_type != BlockSimulationType::evaluateBackward
|
||||
&& simulation_type != BlockSimulationType::evaluateForward)
|
||||
{
|
||||
ostringstream i_output, j_output, v_output;
|
||||
for (int line_counter { ARRAY_SUBSCRIPT_OFFSET(output_type) };
|
||||
const auto &[indices, d] : blocks_derivatives[blk])
|
||||
{
|
||||
const auto &[eq, var, ignore] {indices};
|
||||
i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1-block_recursive_size
|
||||
<< ';' << endl;
|
||||
j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << var+1-block_recursive_size
|
||||
<< ';' << endl;
|
||||
v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
|
||||
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
|
||||
d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
||||
v_output << ';' << endl;
|
||||
line_counter++;
|
||||
}
|
||||
output << i_output.str() << j_output.str() << v_output.str();
|
||||
}
|
||||
}
|
||||
|
||||
template<bool julia>
|
||||
void
|
||||
StaticModel::writeParamsDerivativesFile(const string &basename) const
|
||||
|
|
|
|||
Loading…
Reference in New Issue