From 479c2c029f2fdaa7d6608ad482365691abd3c3de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Mon, 22 Jun 2020 12:46:33 +0200 Subject: [PATCH] Block decomposition: move core of the routine for writing per-block files in separate function This is a preparatory step to allow use_dll with block decomposition. --- src/DynamicModel.cc | 459 ++++++++++++++++++++++++-------------------- src/DynamicModel.hh | 2 + src/StaticModel.cc | 206 +++++++++++--------- src/StaticModel.hh | 3 + 4 files changed, 368 insertions(+), 302 deletions(-) diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index d56a00ce..5abc8f43 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -223,11 +223,232 @@ DynamicModel::additionalBlockTemporaryTerms(int blk, d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count); } +void +DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const +{ + BlockSimulationType simulation_type = blocks[blk].simulation_type; + int block_size = blocks[blk].size; + int block_mfs_size = blocks[blk].mfs_size; + int block_recursive_size = blocks[blk].getRecursiveSize(); + + deriv_node_temp_terms_t tef_terms; + + auto write_eq_tt = [&](int eq) + { + for (auto it : blocks_temporary_terms[blk][eq]) + { + if (dynamic_cast(it)) + it->writeExternalFunctionOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); + + output << " "; + it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms); + output << '='; + it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); + temporary_terms.insert(it); + output << ';' << endl; + } + }; + + // The equations + for (int eq = 0; eq < block_size; eq++) + { + write_eq_tt(eq); + + EquationType equ_type = getBlockEquationType(blk, eq); + BinaryOpNode *e = getBlockEquationExpr(blk, eq); + expr_t lhs = e->arg1, rhs = e->arg2; + switch (simulation_type) + { + case BlockSimulationType::evaluateBackward: + case BlockSimulationType::evaluateForward: + evaluation: + if (equ_type == EquationType::evaluateRenormalized) + { + e = getBlockEquationRenormalizedExpr(blk, eq); + lhs = e->arg1; + rhs = e->arg2; + } + else if (equ_type != EquationType::evaluate) + { + cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl; + exit(EXIT_FAILURE); + } + output << " "; + lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << '='; + rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << ';' << endl; + break; + case BlockSimulationType::solveBackwardSimple: + case BlockSimulationType::solveForwardSimple: + case BlockSimulationType::solveBackwardComplete: + case BlockSimulationType::solveForwardComplete: + case BlockSimulationType::solveTwoBoundariesComplete: + case BlockSimulationType::solveTwoBoundariesSimple: + if (eq < block_recursive_size) + goto evaluation; + output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type) + << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=("; + goto end; + default: + end: + lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << ")-("; + rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << ");" << endl; + } + } + + // The Jacobian if we have to solve the block + + // Write temporary terms for derivatives + write_eq_tt(blocks[blk].size); + + 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]) + { + 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+1); + 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]) + { + 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+1); + 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]) + { + 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+1); + 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]) + { + 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+1); + output << i_output.str() << j_output.str() << v_output.str(); + + // Deterministic mode + if (simulation_type != BlockSimulationType::evaluateForward + && simulation_type != BlockSimulationType::evaluateBackward) + { + 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]) + { + auto [eq, var, lag] = indices; + if (lag == 0) + { + 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) << '=' + << 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]) + { + 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+1); + output << i_output.str() << j_output.str() << v_output.str(); + } + output << " end" << endl; +} + void DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const { temporary_terms_t temporary_terms; // Temp terms written so far - constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicModel; for (int blk = 0; blk < static_cast(blocks.size()); blk++) { @@ -282,234 +503,58 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const << BlockSim(simulation_type) << " //" << endl << " % ////////////////////////////////////////////////////////////////////////" << endl; - if (simulation_type == BlockSimulationType::solveBackwardSimple - || simulation_type == BlockSimulationType::solveForwardSimple - || simulation_type == BlockSimulationType::solveBackwardComplete - || simulation_type == BlockSimulationType::solveForwardComplete - || simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) + if (simulation_type != BlockSimulationType::evaluateForward + && simulation_type != BlockSimulationType::evaluateBackward) output << " residual=zeros(" << block_mfs_size << ",1);" << endl; - deriv_node_temp_terms_t tef_terms; + 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; - auto write_eq_tt = [&](int eq) - { - for (auto it : blocks_temporary_terms[blk][eq]) - { - if (dynamic_cast(it)) - it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); + writeDynamicPerBlockHelper(blk, output, ExprNodeOutputType::matlabDynamicModel, temporary_terms, + nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo); - output << " "; - it->writeOutput(output, local_output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms); - output << " = "; - it->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); - temporary_terms.insert(it); - output << ";" << endl; - } - }; - - // The equations - for (int eq = 0; eq < block_size; eq++) - { - write_eq_tt(eq); - - EquationType equ_type = getBlockEquationType(blk, eq); - BinaryOpNode *e = getBlockEquationExpr(blk, eq); - expr_t lhs = e->arg1, rhs = e->arg2; - switch (simulation_type) - { - case BlockSimulationType::evaluateBackward: - case BlockSimulationType::evaluateForward: - evaluation: - if (equ_type == EquationType::evaluateRenormalized) - { - e = getBlockEquationRenormalizedExpr(blk, eq); - lhs = e->arg1; - rhs = e->arg2; - } - else if (equ_type != EquationType::evaluate) - { - cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl; - exit(EXIT_FAILURE); - } - output << " "; - lhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << " = "; - rhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << ";" << endl; - break; - case BlockSimulationType::solveBackwardSimple: - case BlockSimulationType::solveForwardSimple: - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - case BlockSimulationType::solveTwoBoundariesComplete: - case BlockSimulationType::solveTwoBoundariesSimple: - if (eq < block_recursive_size) - goto evaluation; - output << " residual(" << eq+1-block_recursive_size << ") = ("; - goto end; - default: - end: - lhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << ") - ("; - rhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << ");" << endl; - } - } - - // The Jacobian if we have to solve the block - - // Write temporary terms for derivatives - write_eq_tt(blocks[blk].size); - - output << " % Jacobian" << endl + output << endl << " 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; - ostringstream i_output, j_output, v_output; - int line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives[blk]) - { - auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag }); - i_output << " g1_i(" << line_counter << ")="<< eq+1 << ";" << endl; - j_output << " g1_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl; - v_output << " g1_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - assert(line_counter == nze_stochastic+1); - output << i_output.str() << j_output.str() << v_output.str(); - i_output.str(""); - j_output.str(""); - v_output.str(""); - - line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives_exo[blk]) - { - auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag }); - i_output << " g1_x_i(" << line_counter << ")="<< eq+1 << ";" << endl; - j_output << " g1_x_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl; - v_output << " g1_x_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - assert(line_counter == nze_exo+1); - output << i_output.str() << j_output.str() << v_output.str(); - i_output.str(""); - j_output.str(""); - v_output.str(""); - - line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives_exo_det[blk]) - { - auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag }); - i_output << " g1_xd_i(" << line_counter << ")="<< eq+1 << ";" << endl; - j_output << " g1_xd_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl; - v_output << " g1_xd_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - assert(line_counter == nze_exo_det+1); - output << i_output.str() << j_output.str() << v_output.str(); - i_output.str(""); - j_output.str(""); - v_output.str(""); - - line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives_other_endo[blk]) - { - auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag }); - i_output << " g1_o_i(" << line_counter << ")="<< eq+1 << ";" << endl; - j_output << " g1_o_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl; - v_output << " g1_o_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - assert(line_counter == nze_other_endo+1); - output << i_output.str() << j_output.str() << v_output.str(); - i_output.str(""); - j_output.str(""); - v_output.str(""); - - output << " g1=sparse(g1_i, g1_j, g1_v, " << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ");" << 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; + output << " g1=[];" << endl; break; case BlockSimulationType::solveBackwardSimple: case BlockSimulationType::solveForwardSimple: case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: - output << " g1_i = zeros(" << nze_deterministic << ", 1);" << endl - << " g1_j = zeros(" << nze_deterministic << ", 1);" << endl - << " g1_v = zeros(" << nze_deterministic << ", 1);" << endl; - line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives[blk]) - { - auto [eq, var, lag] = indices; - if (lag == 0) - { - i_output << " g1_i(" << line_counter << ")="<< eq+1 << ";" << endl; - j_output << " g1_j(" << line_counter << ")="<< var+1-block_recursive_size << ";" << endl; - v_output << " g1_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - } - assert(line_counter == nze_deterministic+1); - output << i_output.str() << j_output.str() << v_output.str() - << " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size + 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_i = zeros(" << nze_deterministic << ", 1);" << endl - << " g1_j = zeros(" << nze_deterministic << ", 1);" << endl - << " g1_v = zeros(" << nze_deterministic << ", 1);" << endl; - line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives[blk]) - { - auto [eq, var, lag] = indices; - assert(lag >= -1 && lag <= 1); - if (eq >= block_recursive_size && var >= block_recursive_size) - { - i_output << " g1_i(" << line_counter << ")="<< eq+1-block_recursive_size << ";" << endl; - j_output << " g1_j(" << line_counter << ")="<< var+1-block_recursive_size+block_mfs_size*(lag+1) << ";" << endl; - v_output << " g1_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - } - assert(line_counter == nze_deterministic+1); - output << i_output.str() << j_output.str() << v_output.str() - << " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size + output << " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size << ", " << 3*block_mfs_size << ");" << endl; break; default: diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index d8345b70..42548b74 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -133,6 +133,8 @@ private: void writeDynamicModel(const string &basename, ostream &DynamicOutput, bool use_dll, bool julia) const; //! Writes the main dynamic function of block decomposed model (MATLAB version) void writeDynamicBlockMFile(const string &basename) const; + //! Helper for writing the per-block dynamic files of block decomposed models + void writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, 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 version) void writeDynamicPerBlockMFiles(const string &basename) const; //! Writes the code of the block-decomposed model in virtual machine bytecode diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 22679e27..87e09b09 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -127,17 +127,115 @@ StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instr } } +void +StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms) const +{ + BlockSimulationType simulation_type = blocks[blk].simulation_type; + int block_recursive_size = blocks[blk].getRecursiveSize(); + + // The equations + deriv_node_temp_terms_t tef_terms; + + auto write_eq_tt = [&](int eq) + { + for (auto it : blocks_temporary_terms[blk][eq]) + { + if (dynamic_cast(it)) + it->writeExternalFunctionOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); + + output << " "; + it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms); + output << '='; + it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); + temporary_terms.insert(it); + output << ';' << endl; + } + }; + + for (int eq = 0; eq < blocks[blk].size; eq++) + { + write_eq_tt(eq); + + EquationType equ_type = getBlockEquationType(blk, eq); + BinaryOpNode *e = getBlockEquationExpr(blk, eq); + expr_t lhs = e->arg1, rhs = e->arg2; + switch (simulation_type) + { + case BlockSimulationType::evaluateBackward: + case BlockSimulationType::evaluateForward: + evaluation: + if (equ_type == EquationType::evaluateRenormalized) + { + e = getBlockEquationRenormalizedExpr(blk, eq); + lhs = e->arg1; + rhs = e->arg2; + } + else if (equ_type != EquationType::evaluate) + { + cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl; + exit(EXIT_FAILURE); + } + output << " "; + lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << '='; + rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << ';' << endl; + break; + case BlockSimulationType::solveBackwardSimple: + case BlockSimulationType::solveForwardSimple: + case BlockSimulationType::solveBackwardComplete: + case BlockSimulationType::solveForwardComplete: + if (eq < block_recursive_size) + goto evaluation; + output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type) + << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=("; + lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << ")-("; + rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs); + output << ");" << endl; + break; + default: + cerr << "Incorrect type for block " << blk+1 << endl; + exit(EXIT_FAILURE); + } + } + // The Jacobian if we have to solve the block + if (simulation_type != BlockSimulationType::evaluateBackward + && simulation_type != BlockSimulationType::evaluateForward) + { + // Write temporary terms for derivatives + write_eq_tt(blocks[blk].size); + + ostringstream i_output, j_output, v_output; + int line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type); + for (const auto &[indices, d] : blocks_derivatives[blk]) + { + 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(); + } +} + void StaticModel::writeStaticPerBlockMFiles(const string &basename) const { temporary_terms_t temporary_terms; // Temp terms written so far - constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabStaticModel; for (int blk = 0; blk < static_cast(blocks.size()); blk++) { - // For a block composed of a single equation determines wether we have to evaluate or to solve the equation BlockSimulationType simulation_type = blocks[blk].simulation_type; - int block_recursive_size = blocks[blk].getRecursiveSize(); string filename = packageDir(basename + ".block") + "/static_" + to_string(blk+1) + ".m"; ofstream output; @@ -163,101 +261,19 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const if (simulation_type != BlockSimulationType::evaluateBackward && simulation_type != BlockSimulationType::evaluateForward) - output << " residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl; + 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; - // The equations - deriv_node_temp_terms_t tef_terms; + writeStaticPerBlockHelper(blk, output, ExprNodeOutputType::matlabStaticModel, temporary_terms); - auto write_eq_tt = [&](int eq) - { - for (auto it : blocks_temporary_terms[blk][eq]) - { - if (dynamic_cast(it)) - it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs, tef_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 << " "; - it->writeOutput(output, local_output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms); - output << " = "; - it->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms); - temporary_terms.insert(it); - output << ";" << endl; - } - }; - - for (int eq = 0; eq < blocks[blk].size; eq++) - { - write_eq_tt(eq); - - EquationType equ_type = getBlockEquationType(blk, eq); - BinaryOpNode *e = getBlockEquationExpr(blk, eq); - expr_t lhs = e->arg1, rhs = e->arg2; - switch (simulation_type) - { - case BlockSimulationType::evaluateBackward: - case BlockSimulationType::evaluateForward: - evaluation: - if (equ_type == EquationType::evaluateRenormalized) - { - e = getBlockEquationRenormalizedExpr(blk, eq); - lhs = e->arg1; - rhs = e->arg2; - } - else if (equ_type != EquationType::evaluate) - { - cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl; - exit(EXIT_FAILURE); - } - output << " "; - lhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << " = "; - rhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << ";" << endl; - break; - case BlockSimulationType::solveBackwardSimple: - case BlockSimulationType::solveForwardSimple: - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - if (eq < block_recursive_size) - goto evaluation; - output << " residual(" << eq+1-block_recursive_size << ") = ("; - lhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << ") - ("; - rhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - output << ");" << endl; - break; - default: - cerr << "Incorrect type for block " << blk+1 << endl; - exit(EXIT_FAILURE); - } - } - // The Jacobian if we have to solve the block - if (simulation_type == BlockSimulationType::solveBackwardSimple - || simulation_type == BlockSimulationType::solveForwardSimple - || simulation_type == BlockSimulationType::solveBackwardComplete - || simulation_type == BlockSimulationType::solveForwardComplete) - { - // Write temporary terms for derivatives - write_eq_tt(blocks[blk].size); - - output << " 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; - - ostringstream i_output, j_output, v_output; - int line_counter = 1; - for (const auto &[indices, d] : blocks_derivatives[blk]) - { - auto [eq, var, ignore] = indices; - i_output << " g1_i(" << line_counter << ")=" << eq+1-block_recursive_size << ";" << endl; - j_output << " g1_j(" << line_counter << ")=" << var+1-block_recursive_size << ";"; - v_output << " g1_v(" << line_counter << ")="; - d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); - v_output << ";" << endl; - line_counter++; - } - output << i_output.str() << j_output.str() << v_output.str() - << " g1=sparse(g1_i, g1_j, g1_v, " << blocks[blk].mfs_size << "," << blocks[blk].mfs_size << ");" << endl; - } output << "end" << endl; output.close(); } diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 3ebf93a9..12327ca7 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -48,6 +48,9 @@ private: //! Writes the main static function of block decomposed model (MATLAB version) void writeStaticBlockMFile(const string &basename) const; + //! Helper for writing a per-block static file of block decomposed model + void writeStaticPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms) const; + //! Writes the per-block static files of block decomposed model (MATLAB version) void writeStaticPerBlockMFiles(const string &basename) const;