diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index ea49debc..428a0cf5 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -224,104 +224,88 @@ DynamicModel::additionalBlockTemporaryTerms(int blk, } void -DynamicModel::writeModelEquationsOrdered_M(const string &basename) const +DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const { - string tmp_s, sps; - ostringstream tmp_output, tmp1_output, global_output; - expr_t lhs = nullptr, rhs = nullptr; - BinaryOpNode *eq_node; - ostringstream Ufoss; - vector Uf(symbol_table.endo_nbr(), ""); - map reference_count; - ofstream output; - int nze, nze_exo, nze_exo_det, nze_other_endo; - vector feedback_variables; - temporary_terms_t temporary_terms; // Temp terms written so far - ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicModelSparse; + constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicModelSparse; - //---------------------------------------------------------------------- - //For each block - for (int block = 0; block < static_cast(blocks.size()); block++) + for (int blk = 0; blk < static_cast(blocks.size()); blk++) { + int nze = blocks_derivatives[blk].size(); + 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(); + 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(); - //recursive_variables.clear(); - feedback_variables.clear(); - //For a block composed of a single equation determines wether we have to evaluate or to solve the equation - nze = blocks_derivatives[block].size(); - nze_other_endo = blocks_derivatives_other_endo[block].size(); - nze_exo = blocks_derivatives_exo[block].size(); - nze_exo_det = blocks_derivatives_exo_det[block].size(); - BlockSimulationType simulation_type = blocks[block].simulation_type; - int block_size = blocks[block].size; - int block_mfs = blocks[block].mfs_size; - int block_recursive = blocks[block].getRecursiveSize(); - local_output_type = ExprNodeOutputType::matlabDynamicModelSparse; + string filename = packageDir(basename + ".block") + "/dynamic_" + to_string(blk+1) + ".m"; + ofstream output; + output.open(filename, ios::out | ios::binary); + + vector Uf(block_size); - tmp1_output.str(""); - tmp1_output << packageDir(basename + ".block") << "/dynamic_" << block+1 << ".m"; - output.open(tmp1_output.str(), ios::out | ios::binary); output << "%" << endl - << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare" << endl + << "% " << filename << " : Computes dynamic version of one block" << endl << "%" << endl << "% Warning : this file is generated automatically by Dynare" << endl << "% from model file (.mod)" << endl << endl - << "%/" << endl; + << "%" << endl; if (simulation_type == BlockSimulationType::evaluateBackward || simulation_type == BlockSimulationType::evaluateForward) - { - output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl; - } + output << "function [y, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl; else if (simulation_type == BlockSimulationType::solveForwardComplete || simulation_type == BlockSimulationType::solveBackwardComplete) - output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl; + output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl; else if (simulation_type == BlockSimulationType::solveBackwardSimple || simulation_type == BlockSimulationType::solveForwardSimple) - output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl; + output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl; else - output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl; + output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl; output << " % ////////////////////////////////////////////////////////////////////////" << endl - << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 + << " % //" << string(" Block ").substr(static_cast(log10(blk + 1))) << blk+1 << " //" << endl << " % // Simulation type " << BlockSim(simulation_type) << " //" << endl << " % ////////////////////////////////////////////////////////////////////////" << endl; - //The Temporary terms + if (simulation_type == BlockSimulationType::evaluateBackward || simulation_type == BlockSimulationType::evaluateForward) { - output << " if(jacobian_eval)" << endl - << " g1 = spalloc(" << block_mfs << ", " << blocks_jacob_cols_endo[block].size() << ", " << nze << ");" << endl - << " g1_x=spalloc(" << block_size << ", " << blocks_jacob_cols_exo[block].size() << ", " << nze_exo << ");" << endl - << " g1_xd=spalloc(" << block_size << ", " << blocks_jacob_cols_exo_det[block].size() << ", " << nze_exo_det << ");" << endl - << " g1_o=spalloc(" << block_size << ", " << blocks_jacob_cols_other_endo[block].size() << ", " << nze_other_endo << ");" << endl + output << " if jacobian_eval" << endl + << " g1 = spalloc(" << block_mfs_size << ", " << blocks_jacob_cols_endo[blk].size() << ", " << nze << ");" << endl + << " g1_x = spalloc(" << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ", " << nze_exo << ");" << endl + << " g1_xd = spalloc(" << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ", " << nze_exo_det << ");" << endl + << " g1_o = spalloc(" << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ", " << nze_other_endo << ");" << endl << " end;" << endl; } else { - output << " if(jacobian_eval)" << endl - << " g1 = spalloc(" << block_size << ", " << blocks_jacob_cols_endo[block].size() << ", " << nze << ");" << endl - << " g1_x=spalloc(" << block_size << ", " << blocks_jacob_cols_exo[block].size() << ", " << nze_exo << ");" << endl - << " g1_xd=spalloc(" << block_size << ", " << blocks_jacob_cols_exo_det[block].size() << ", " << nze_exo_det << ");" << endl - << " g1_o=spalloc(" << block_size << ", " << blocks_jacob_cols_other_endo[block].size() << ", " << nze_other_endo << ");" << endl + output << " if jacobian_eval" << endl + << " g1 = spalloc(" << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ", " << nze << ");" << endl + << " g1_x = spalloc(" << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ", " << nze_exo << ");" << endl + << " g1_xd = spalloc(" << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ", " << nze_exo_det << ");" << endl + << " g1_o = spalloc(" << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ", " << nze_other_endo << ");" << endl << " else" << endl; if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - output << " g1 = spalloc(" << block_mfs << "*Periods, " - << block_mfs << "*(Periods+" << blocks[block].max_lag+blocks[block].max_lead+1 << ")" + output << " g1 = spalloc(" << block_mfs_size << "*Periods, " + << block_mfs_size << "*(Periods+" << blocks[blk].max_lag+blocks[blk].max_lead+1 << ")" << ", " << nze << "*Periods);" << endl; else - output << " g1 = spalloc(" << block_mfs - << ", " << block_mfs << ", " << nze << ");" << endl; - output << " end;" << endl; + output << " g1 = spalloc(" << block_mfs_size + << ", " << block_mfs_size << ", " << nze << ");" << endl; + output << " end" << endl; } - output << " g2=0;g3=0;" << endl; + output << " g2=0;" << endl + << " g3=0;" << endl; // Declare global temp terms from this block and the previous ones bool global_keyword_written = false; - for (int blk2 = 0; blk2 <= block; blk2++) + for (int blk2 = 0; blk2 <= blk; blk2++) for (auto &eq_tt : blocks_temporary_terms[blk2]) for (auto tt : eq_tt) { @@ -341,9 +325,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) { - output << " " << "% //Temporary variables initialization" << endl - << " " << "T_zeros = zeros(y_kmin+periods, 1);" << endl; - for (auto &eq_tt : blocks_temporary_terms[block]) + output << " " << "T_zeros = zeros(y_kmin+periods, 1);" << endl; + for (auto &eq_tt : blocks_temporary_terms[blk]) for (auto tt : eq_tt) { output << " "; @@ -356,10 +339,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const || simulation_type == BlockSimulationType::solveForwardSimple || simulation_type == BlockSimulationType::solveBackwardComplete || simulation_type == BlockSimulationType::solveForwardComplete) - output << " residual=zeros(" << block_mfs << ",1);" << endl; + output << " residual=zeros(" << block_mfs_size << ",1);" << endl; else if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - output << " residual=zeros(" << block_mfs << ",y_kmin+periods);" << endl; + output << " residual=zeros(" << block_mfs_size << ",y_kmin+periods);" << endl; if (simulation_type == BlockSimulationType::evaluateBackward) output << " for it_ = (y_kmin+periods):y_kmin+1" << endl; if (simulation_type == BlockSimulationType::evaluateForward) @@ -367,32 +350,30 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - { - output << " b = zeros(periods*y_size,1);" << endl - << " for it_ = y_kmin+1:(periods+y_kmin)" << endl - << " Per_y_=it_*y_size;" << endl - << " Per_J_=(it_-y_kmin-1)*y_size;" << endl - << " Per_K_=(it_-1)*y_size;" << endl; - sps = " "; - } - else - if (simulation_type == BlockSimulationType::evaluateBackward - || simulation_type == BlockSimulationType::evaluateForward) - sps = " "; - else - sps = ""; + output << " b = zeros(periods*y_size,1);" << endl + << " for it_ = y_kmin+1:(periods+y_kmin)" << endl + << " Per_y_=it_*y_size;" << endl + << " Per_J_=(it_-y_kmin-1)*y_size;" << endl + << " Per_K_=(it_-1)*y_size;" << endl; + + string indent_prefix = " "; + if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete + || simulation_type == BlockSimulationType::solveTwoBoundariesSimple + || simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + indent_prefix += " "; deriv_node_temp_terms_t tef_terms; auto write_eq_tt = [&](int eq) { - for (auto it : blocks_temporary_terms[block][eq]) + for (auto it : blocks_temporary_terms[blk][eq]) { if (dynamic_cast(it)) it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, temporary_terms_idxs, tef_terms); - output << " " << sps; - it->writeOutput(output, local_output_type, blocks_temporary_terms[block][eq], {}, tef_terms); + output << indent_prefix; + it->writeOutput(output, local_output_type, blocks_temporary_terms[blk][eq], {}, tef_terms); output << " = "; it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms); temporary_terms.insert(it); @@ -401,258 +382,176 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const }; // The equations - for (int i = 0; i < block_size; i++) + for (int eq = 0; eq < block_size; eq++) { - write_eq_tt(i); + write_eq_tt(eq); - int variable_ID = getBlockVariableID(block, i); - int equation_ID = getBlockEquationID(block, i); - EquationType equ_type = getBlockEquationType(block, i); - string sModel = symbol_table.getName(symbol_table.getID(SymbolType::endogenous, variable_ID)); - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - tmp_output.str(""); - lhs->writeOutput(tmp_output, local_output_type, temporary_terms, {}); + 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 (simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - output << " % equation " << getBlockEquationID(block, i)+1 << " variable : " << sModel - << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl; - output << " "; - if (equ_type == EquationType::evaluate) + evaluation: + if (equ_type == EquationType::evaluate_s) { - output << tmp_output.str(); - output << " = "; - rhs->writeOutput(output, local_output_type, temporary_terms, {}); + e = getBlockEquationRenormalizedExpr(blk, eq); + lhs = e->arg1; + rhs = e->arg2; } - else if (equ_type == EquationType::evaluate_s) + else if (equ_type != EquationType::evaluate) { - output << "%" << tmp_output.str(); - output << " = "; - if (isBlockEquationRenormalized(block, i)) - { - rhs->writeOutput(output, local_output_type, temporary_terms, {}); - output << endl << " "; - tmp_output.str(""); - eq_node = getBlockEquationRenormalizedExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - lhs->writeOutput(output, local_output_type, temporary_terms, {}); - output << " = "; - rhs->writeOutput(output, local_output_type, temporary_terms, {}); - } - } - else - { - cerr << "Type mismatch for equation " << equation_ID+1 << endl; + cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl; exit(EXIT_FAILURE); } + output << indent_prefix; + lhs->writeOutput(output, local_output_type, temporary_terms, {}); + output << " = "; + rhs->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; break; case BlockSimulationType::solveBackwardSimple: case BlockSimulationType::solveForwardSimple: case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: - if (i < block_recursive) + if (eq < block_recursive_size) goto evaluation; - feedback_variables.push_back(variable_ID); - output << " % equation " << equation_ID+1 << " variable : " << sModel - << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl; - output << " " << "residual(" << i+1-block_recursive << ") = ("; + output << indent_prefix << "residual(" << eq+1-block_recursive_size << ") = ("; goto end; case BlockSimulationType::solveTwoBoundariesComplete: case BlockSimulationType::solveTwoBoundariesSimple: - if (i < block_recursive) + if (eq < block_recursive_size) goto evaluation; - feedback_variables.push_back(variable_ID); - output << " % equation " << equation_ID+1 << " variable : " << sModel - << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl; - Ufoss << " b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)"; - Uf[equation_ID] += Ufoss.str(); - Ufoss.str(""); - output << " residual(" << i+1-block_recursive << ", it_) = ("; + Uf[eq] << "b(" << eq+1-block_recursive_size << "+Per_J_) = -residual(" << eq+1-block_recursive_size << ", it_)"; + output << indent_prefix << "residual(" << eq+1-block_recursive_size << ", it_) = ("; goto end; default: end: - output << tmp_output.str(); + lhs->writeOutput(output, local_output_type, temporary_terms, {}); output << ") - ("; rhs->writeOutput(output, local_output_type, temporary_terms, {}); output << ");" << endl; -#ifdef CONDITION - if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - output << " condition(" << i+1 << ")=0;" << endl; -#endif } } + // The Jacobian if we have to solve the block // Write temporary terms for derivatives - write_eq_tt(blocks[block].size); + write_eq_tt(blocks[blk].size); - if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple - || simulation_type == BlockSimulationType::solveTwoBoundariesComplete) - output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl; - else - if (simulation_type == BlockSimulationType::solveBackwardSimple - || simulation_type == BlockSimulationType::solveForwardSimple - || simulation_type == BlockSimulationType::solveBackwardComplete - || simulation_type == BlockSimulationType::solveForwardComplete) - output << " % Jacobian " << endl << " if jacobian_eval" << endl; - else - output << " % Jacobian " << endl << " if jacobian_eval" << endl; - for (const auto &[indices, d] : blocks_derivatives[block]) + output << indent_prefix << "% Jacobian" << endl + << indent_prefix << "if jacobian_eval" << endl; + indent_prefix += " "; + for (const auto &[indices, d] : blocks_derivatives[blk]) { auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_endo[block].at({ var, lag }); - output << " g1(" << eq+1 << ", " << jacob_col+1 << ") = "; + int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag }); + output << indent_prefix << "g1(" << eq+1 << ", " << jacob_col+1 << ") = "; d->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; } - for (const auto &[indices, d] : blocks_derivatives_exo[block]) + for (const auto &[indices, d] : blocks_derivatives_exo[blk]) { auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_exo[block].at({ var, lag }); - output << " g1_x(" << eq+1 << ", " << jacob_col+1 << ") = "; + int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag }); + output << indent_prefix << "g1_x(" << eq+1 << ", " << jacob_col+1 << ") = "; d->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; } - for (const auto &[indices, d] : blocks_derivatives_exo_det[block]) + for (const auto &[indices, d] : blocks_derivatives_exo_det[blk]) { auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_exo_det[block].at({ var, lag }); - output << " g1_xd(" << eq+1 << ", " << jacob_col+1 << ") = "; + int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag }); + output << indent_prefix << "g1_xd(" << eq+1 << ", " << jacob_col+1 << ") = "; d->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; } - for (const auto &[indices, d] : blocks_derivatives_other_endo[block]) + for (const auto &[indices, d] : blocks_derivatives_other_endo[blk]) { auto [eq, var, lag] = indices; - int jacob_col = blocks_jacob_cols_other_endo[block].at({ var, lag }); - output << " g1_o(" << eq+1 << ", " << jacob_col+1 << ") = "; + int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag }); + output << indent_prefix << "g1_o(" << eq+1 << ", " << jacob_col+1 << ") = "; d->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; } - output << " varargout{1}=g1_x;" << endl - << " varargout{2}=g1_xd;" << endl - << " varargout{3}=g1_o;" << endl; + output << indent_prefix << "varargout{1}=g1_x;" << endl + << indent_prefix << "varargout{2}=g1_xd;" << endl + << indent_prefix << "varargout{3}=g1_o;" << endl; switch (simulation_type) { case BlockSimulationType::evaluateForward: case BlockSimulationType::evaluateBackward: - output << " end;" << endl - << " end;" << endl; + output << " end" << endl + << " end" << endl; break; case BlockSimulationType::solveBackwardSimple: case BlockSimulationType::solveForwardSimple: case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: output << " else" << endl; - for (const auto &[indices, id] : blocks_derivatives[block]) + for (const auto &[indices, d] : blocks_derivatives[blk]) { auto [eq, var, lag] = indices; - int eqr = getBlockEquationID(block, eq); - int varr = getBlockVariableID(block, var); if (lag == 0) { - output << " g1(" << eq+1 << ", " << var+1-block_recursive << ") = "; - id->writeOutput(output, local_output_type, temporary_terms, {}); - output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr)) - << "(" << lag - << ") " << varr+1 - << ", equation=" << eqr+1 << endl; + output << " g1(" << eq+1 << ", " << var+1-block_recursive_size << ") = "; + d->writeOutput(output, local_output_type, temporary_terms, {}); + output << ";" << endl; } - } - output << " end;" << endl; + output << " end" << endl; break; case BlockSimulationType::solveTwoBoundariesSimple: case BlockSimulationType::solveTwoBoundariesComplete: output << " else" << endl; - for (const auto &[indices, id] : blocks_derivatives[block]) + for (const auto &[indices, d] : blocks_derivatives[blk]) { auto [eq, var, lag] = indices; - int eqr = getBlockEquationID(block, eq); - int varr = getBlockVariableID(block, var); - ostringstream tmp_output; - if (eq >= block_recursive && var >= block_recursive) + int varr = getBlockVariableID(blk, var); + if (eq >= block_recursive_size && var >= block_recursive_size) { if (lag == 0) - Ufoss << "+g1(" << eq+1-block_recursive - << "+Per_J_, " << var+1-block_recursive - << "+Per_K_)*y(it_, " << varr+1 << ")"; + Uf[eq] << "+g1(" << eq+1-block_recursive_size + << "+Per_J_, " << var+1-block_recursive_size + << "+Per_K_)*y(it_, " << varr+1 << ")"; else if (lag == 1) - Ufoss << "+g1(" << eq+1-block_recursive - << "+Per_J_, " << var+1-block_recursive - << "+Per_y_)*y(it_+1, " << varr+1 << ")"; + Uf[eq] << "+g1(" << eq+1-block_recursive_size + << "+Per_J_, " << var+1-block_recursive_size + << "+Per_y_)*y(it_+1, " << varr+1 << ")"; else if (lag > 0) - Ufoss << "+g1(" << eq+1-block_recursive - << "+Per_J_, " << var+1-block_recursive - << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")"; + Uf[eq] << "+g1(" << eq+1-block_recursive_size + << "+Per_J_, " << var+1-block_recursive_size + << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")"; else - Ufoss << "+g1(" << eq+1-block_recursive - << "+Per_J_, " << var+1-block_recursive - << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")"; - Uf[eqr] += Ufoss.str(); - Ufoss.str(""); + Uf[eq] << "+g1(" << eq+1-block_recursive_size + << "+Per_J_, " << var+1-block_recursive_size + << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")"; + output << " "; if (lag == 0) - tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, " - << var+1-block_recursive << "+Per_K_) = "; + output << " g1(" << eq+1-block_recursive_size << "+Per_J_, " + << var+1-block_recursive_size << "+Per_K_) = "; else if (lag == 1) - tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, " - << var+1-block_recursive << "+Per_y_) = "; + output << " g1(" << eq+1-block_recursive_size << "+Per_J_, " + << var+1-block_recursive_size << "+Per_y_) = "; else if (lag > 0) - tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, " - << var+1-block_recursive << "+y_size*(it_+" << lag-1 << ")) = "; + output << " g1(" << eq+1-block_recursive_size << "+Per_J_, " + << var+1-block_recursive_size << "+y_size*(it_+" << lag-1 << ")) = "; else if (lag < 0) - tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, " - << var+1-block_recursive << "+y_size*(it_" << lag-1 << ")) = "; - output << " " << tmp_output.str(); - id->writeOutput(output, local_output_type, temporary_terms, {}); - output << ";"; - output << " %2 variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr)) - << "(" << lag << ") " << varr+1 - << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl; + output << " g1(" << eq+1-block_recursive_size << "+Per_J_, " + << var+1-block_recursive_size << "+y_size*(it_" << lag-1 << ")) = "; + d->writeOutput(output, local_output_type, temporary_terms, {}); + output << ";" << endl; } + } + for (int eq = 0; eq < block_size; eq++) + if (eq >= block_recursive_size) + output << " " << Uf[eq].str() << ";" << endl; -#ifdef CONDITION - output << " if (fabs(condition[" << eqr << "])= block_recursive) - output << " " << Uf[getBlockEquationID(block, i)] << ";" << endl; -#ifdef CONDITION - output << " if (fabs(condition(" << i+1 << "))Block_List[block].Max_Lead+ModelBlock->Block_List[block].Max_Lag; m++) - { - k = m-ModelBlock->Block_List[block].Max_Lag; - for (i = 0; i < ModelBlock->Block_List[block].IM_lead_lag[m].size; i++) - { - int eq = ModelBlock->Block_List[block].IM_lead_lag[m].Equ_Index[i]; - int var = ModelBlock->Block_List[block].IM_lead_lag[m].Var_Index[i]; - int u = ModelBlock->Block_List[block].IM_lead_lag[m].u[i]; - int eqr = ModelBlock->Block_List[block].IM_lead_lag[m].Equ[i]; - output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");" << endl; - } - } - for (i = 0; i < ModelBlock->Block_List[block].Size; i++) - output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");" << endl; -#endif - output << " end;" << endl - << " end;" << endl; + output << " end" << endl + << " end" << endl; break; default: break; @@ -1537,7 +1436,7 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, int num, } void -DynamicModel::writeSparseDynamicMFile(const string &basename) const +DynamicModel::writeDynamicBlockMFile(const string &basename) const { string sp; ofstream mDynamicModelFile; @@ -1857,8 +1756,6 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const << "end" << endl; mDynamicModelFile.close(); - - writeModelEquationsOrdered_M(basename); } void @@ -4748,7 +4645,10 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_d writeModelEquationsCode(basename); } else if (block && !bytecode) - writeSparseDynamicMFile(basename); + { + writeDynamicPerBlockMFiles(basename); + writeDynamicBlockMFile(basename); + } else if (use_dll) { writeDynamicCFile(basename); diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 698b22e0..cca44853 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -126,15 +126,15 @@ private: //! Writes dynamic model file (C version) /*! \todo add third derivatives handling */ void writeDynamicCFile(const string &basename) const; - //! Writes dynamic model file when SparseDLL option is on - void writeSparseDynamicMFile(const string &basename) const; //! Writes the dynamic model equations and its derivatives /*! \todo add third derivatives handling in C output */ void writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const; void writeDynamicModel(const string &basename, bool use_dll, bool julia) const; void writeDynamicModel(const string &basename, ostream &DynamicOutput, bool use_dll, bool julia) const; - //! Writes the Block reordred structure of the model in M output - void writeModelEquationsOrdered_M(const string &basename) const; + //! Writes the main dynamic function of block decomposed model (MATLAB version) + void writeDynamicBlockMFile(const string &basename) 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 reordred structure of the model in virtual machine bytecode void writeModelEquationsCode_Block(const string &basename, bool linear_decomposition) const; //! Writes the code of the model in virtual machine bytecode diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 5ac211b6..ba7c9e2d 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -128,51 +128,47 @@ StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instr } void -StaticModel::writeModelEquationsOrdered_M(const string &basename) const +StaticModel::writeStaticPerBlockMFiles(const string &basename) const { temporary_terms_t temporary_terms; // Temp terms written so far constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabStaticModelSparse; - //---------------------------------------------------------------------- - //For each block - for (int block = 0; block < static_cast(blocks.size()); block++) + 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[block].simulation_type; - int block_size = blocks[block].size; - int block_mfs = blocks[block].mfs_size; - int block_recursive = blocks[block].getRecursiveSize(); + BlockSimulationType simulation_type = blocks[blk].simulation_type; + int block_recursive_size = blocks[blk].getRecursiveSize(); - string filename = packageDir(basename + ".block") + "/static_" + to_string(block+1) + ".m"; + string filename = packageDir(basename + ".block") + "/static_" + to_string(blk+1) + ".m"; ofstream output; output.open(filename, ios::out | ios::binary); output << "%" << endl - << "% " << filename << " : Computes static model for Dynare" << endl + << "% " << filename << " : Computes static version of one block" << endl << "%" << endl << "% Warning : this file is generated automatically by Dynare" << endl << "% from model file (.mod)" << endl << endl - << "%/" << endl; + << "%" << endl; if (simulation_type == BlockSimulationType::evaluateBackward || simulation_type == BlockSimulationType::evaluateForward) - output << "function y = static_" << block+1 << "(y, x, params)" << endl; + output << "function y = static_" << blk+1 << "(y, x, params)" << endl; else - output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)" << endl; + output << "function [residual, y, g1] = static_" << blk+1 << "(y, x, params)" << endl; output << " % ////////////////////////////////////////////////////////////////////////" << endl - << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 + << " % //" << string(" Block ").substr(static_cast(log10(blk + 1))) << blk+1 << " //" << endl << " % // Simulation type " << BlockSim(simulation_type) << " //" << endl - << " % ////////////////////////////////////////////////////////////////////////" << endl - << " global options_;" << endl; - // The Temporary terms + << " % ////////////////////////////////////////////////////////////////////////" << endl; + if (simulation_type != BlockSimulationType::evaluateBackward && simulation_type != BlockSimulationType::evaluateForward) - output << " g1 = spalloc(" << block_mfs << ", " << block_mfs << ", " << blocks_derivatives[block].size() << ");" << endl; + output << " g1=spalloc(" << blocks[blk].mfs_size << "," << blocks[blk].mfs_size + << "," << blocks_derivatives[blk].size() << ");" << endl; // Declare global temp terms from this block and the previous ones bool global_keyword_written = false; - for (int blk2 = 0; blk2 <= block; blk2++) + for (int blk2 = 0; blk2 <= blk; blk2++) for (auto &eq_tt : blocks_temporary_terms[blk2]) for (auto tt : eq_tt) { @@ -189,20 +185,20 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const if (simulation_type != BlockSimulationType::evaluateBackward && simulation_type != BlockSimulationType::evaluateForward) - output << " residual=zeros(" << block_mfs << ",1);" << endl; + output << " residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl; // The equations deriv_node_temp_terms_t tef_terms; auto write_eq_tt = [&](int eq) { - for (auto it : blocks_temporary_terms[block][eq]) + for (auto it : blocks_temporary_terms[blk][eq]) { if (dynamic_cast(it)) it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, {}, tef_terms); output << " "; - it->writeOutput(output, local_output_type, blocks_temporary_terms[block][eq], {}, tef_terms); + it->writeOutput(output, local_output_type, blocks_temporary_terms[blk][eq], {}, tef_terms); output << " = "; it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms); temporary_terms.insert(it); @@ -210,58 +206,49 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const } }; - for (int i = 0; i < block_size; i++) + for (int eq = 0; eq < blocks[blk].size; eq++) { - write_eq_tt(i); + write_eq_tt(eq); - int variable_ID = getBlockVariableID(block, i); - int equation_ID = getBlockEquationID(block, i); - EquationType equ_type = getBlockEquationType(block, i); - string sModel = symbol_table.getName(symbol_table.getID(SymbolType::endogenous, variable_ID)); - BinaryOpNode *eq_node = getBlockEquationExpr(block, i); - expr_t lhs = eq_node->arg1, rhs = eq_node->arg2; - ostringstream tmp_output; - lhs->writeOutput(tmp_output, local_output_type, temporary_terms, {}); + 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: - output << " "; - if (equ_type == EquationType::evaluate) + if (equ_type == EquationType::evaluate_s) { - output << tmp_output.str() << " = "; - rhs->writeOutput(output, local_output_type, temporary_terms, {}); + e = getBlockEquationRenormalizedExpr(blk, eq); + lhs = e->arg1; + rhs = e->arg2; } - else if (equ_type == EquationType::evaluate_s) + else if (equ_type != EquationType::evaluate) { - eq_node = getBlockEquationRenormalizedExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - lhs->writeOutput(output, local_output_type, temporary_terms, {}); - output << " = "; - rhs->writeOutput(output, local_output_type, temporary_terms, {}); - } - else - { - cerr << "Type mismatch for equation " << equation_ID+1 << endl; + cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl; exit(EXIT_FAILURE); } + output << " "; + lhs->writeOutput(output, local_output_type, temporary_terms, {}); + output << " = "; + rhs->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; break; case BlockSimulationType::solveBackwardSimple: case BlockSimulationType::solveForwardSimple: case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: - if (i < block_recursive) + if (eq < block_recursive_size) goto evaluation; - output << " " << "residual(" << i+1-block_recursive << ") = (" - << tmp_output.str() << ") - ("; + output << " residual(" << eq+1-block_recursive_size << ") = ("; + lhs->writeOutput(output, local_output_type, temporary_terms, {}); + output << ") - ("; rhs->writeOutput(output, local_output_type, temporary_terms, {}); output << ");" << endl; break; default: - cerr << "Incorrect type for block " << block+1 << endl; + cerr << "Incorrect type for block " << blk+1 << endl; exit(EXIT_FAILURE); } } @@ -273,13 +260,13 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: // Write temporary terms for derivatives - write_eq_tt(blocks[block].size); + write_eq_tt(blocks[blk].size); - for (const auto &[indices, id] : blocks_derivatives[block]) + for (const auto &[indices, d] : blocks_derivatives[blk]) { auto [eq, var, ignore] = indices; - output << " g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = "; - id->writeOutput(output, local_output_type, temporary_terms, {}); + output << " g1(" << eq+1-block_recursive_size << ", " << var+1-block_recursive_size << ") = "; + d->writeOutput(output, local_output_type, temporary_terms, {}); output << ";" << endl; } break; @@ -1736,8 +1723,8 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, writeModelEquationsCode(basename); else if (block && !bytecode) { - writeModelEquationsOrdered_M(basename); - writeStaticBlockMFSFile(basename); + writeStaticPerBlockMFiles(basename); + writeStaticBlockMFile(basename); } else if (use_dll) { @@ -1761,7 +1748,7 @@ StaticModel::exoPresentInEqs() const } void -StaticModel::writeStaticBlockMFSFile(const string &basename) const +StaticModel::writeStaticBlockMFile(const string &basename) const { string filename = packageDir(basename) + "/static.m"; @@ -1779,29 +1766,27 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const << " var_index = [];" << endl << endl << " switch nblock" << endl; - int nb_blocks = blocks.size(); - - for (int b = 0; b < nb_blocks; b++) + for (int blk = 0; blk < static_cast(blocks.size()); blk++) { set local_var; - output << " case " << b+1 << endl; + output << " case " << blk+1 << endl; - BlockSimulationType simulation_type = blocks[b].simulation_type; + BlockSimulationType simulation_type = blocks[blk].simulation_type; if (simulation_type == BlockSimulationType::evaluateBackward || simulation_type == BlockSimulationType::evaluateForward) { - output << " y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl; - ostringstream tmp; - for (int i = 0; i < blocks[b].size; i++) - tmp << " " << getBlockVariableID(b, i)+1; - output << " var_index = [" << tmp.str() << "];" << endl + output << " y_tmp = " << basename << ".block.static_" << blk+1 << "(y, x, params);" << endl + << " var_index = ["; + for (int var = 0; var < blocks[blk].size; var++) + output << " " << getBlockVariableID(blk, var)+1; + output << "];" << endl << " residual = y(var_index) - y_tmp(var_index);" << endl << " y = y_tmp;" << endl; } else - output << " [residual, y, g1] = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl; + output << " [residual, y, g1] = " << basename << ".block.static_" << blk+1 << "(y, x, params);" << endl; } output << " end" << endl diff --git a/src/StaticModel.hh b/src/StaticModel.hh index c05fef59..a6f7ffd5 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -45,11 +45,11 @@ private: //! Writes the static model equations and its derivatives void writeStaticModel(const string &basename, ostream &StaticOutput, bool use_dll, bool julia) const; - //! Writes the static function calling the block to solve (Matlab version) - void writeStaticBlockMFSFile(const string &basename) const; + //! Writes the main static function of block decomposed model (MATLAB version) + void writeStaticBlockMFile(const string &basename) const; - //! Writes the Block reordred structure of the model in M output - void writeModelEquationsOrdered_M(const string &basename) const; + //! Writes the per-block static files of block decomposed model (MATLAB version) + void writeStaticPerBlockMFiles(const string &basename) const; //! Writes the code of the Block reordred structure of the model in virtual machine bytecode void writeModelEquationsCode_Block(const string &basename) const;