From 2e136cae56c33b4ab734a079f5bfebae16ea7b80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Fri, 24 Apr 2020 12:29:02 +0200 Subject: [PATCH] Block decomposition: simplify computation of block derivatives MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit — use a std::map for storing block derivatives — remove redundant ModelTree::first_chain_rule_derivatives structure — remove unused codepaths in StaticModel — DynamicModel: simplify code that determines the type of derivatives in a block. We now use a slightly different categorization. — by the way, fix the max lead/lag information for blocks that are obtained via merging. A workaround was previously implemented in DynamicModel::get_Derivative(), but it is no longer needed with this fix. --- src/DynamicModel.cc | 203 +++++++++++++++++++------------------------- src/DynamicModel.hh | 22 +++-- src/ModelTree.cc | 23 ++--- src/ModelTree.hh | 16 +--- src/StaticModel.cc | 176 ++++++++++---------------------------- src/StaticModel.hh | 4 +- 6 files changed, 162 insertions(+), 282 deletions(-) diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index a877ce0b..8d4eb129 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -165,9 +165,10 @@ DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_n } void -DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const +DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const map_idx_t &map_idx) const { - if (auto it = first_chain_rule_derivatives.find({ eqr, varr, lag }); it != first_chain_rule_derivatives.end()) + if (auto it = blocks_derivatives[blk].find({ eq, var, lag }); + it != blocks_derivatives[blk].end()) it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false); else { @@ -182,7 +183,6 @@ DynamicModel::computeTemporaryTermsOrdered() map> first_occurence; map reference_count; BinaryOpNode *eq_node; - first_chain_rule_derivatives_t::const_iterator it_chr; ostringstream tmp_s; v_temporary_terms.clear(); map_idx.clear(); @@ -211,11 +211,8 @@ DynamicModel::computeTemporaryTermsOrdered() eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i); } } - for (const auto &it : blocks_derivatives[block]) - { - expr_t id = get<3>(it); - id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); - } + for (const auto &[ignore, id] : blocks_derivatives[block]) + id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); for (const auto &it : derivative_endo[block]) it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); for (const auto &it : derivative_other_endo[block]) @@ -241,11 +238,8 @@ DynamicModel::computeTemporaryTermsOrdered() eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i); } } - for (const auto &it : blocks_derivatives[block]) - { - expr_t id = get<3>(it); - id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); - } + for (const auto &[ignore, id] : blocks_derivatives[block]) + id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); for (const auto &it : derivative_endo[block]) it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); for (const auto &it : derivative_other_endo[block]) @@ -267,11 +261,8 @@ DynamicModel::computeTemporaryTermsOrdered() eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); } } - for (const auto &it : blocks_derivatives[block]) - { - expr_t id = get<3>(it); - id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); - } + for (const auto &[ignore, id] : blocks_derivatives[block]) + id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); for (const auto &it : derivative_endo[block]) it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); for (const auto &it : derivative_other_endo[block]) @@ -339,8 +330,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const int prev_lag; int prev_var, count_col, count_col_endo, count_col_exo, count_col_exo_det, count_col_other_endo; map, expr_t> tmp_block_endo_derivative; - for (const auto &it : blocks_derivatives[block]) - tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it); + for (const auto &[idx, d] : blocks_derivatives[block]) + tmp_block_endo_derivative[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d; prev_var = 999999999; prev_lag = -9999999; count_col_endo = 0; @@ -762,8 +753,9 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: output << " else" << endl; - for (const auto &[eq, var, lag, id] : blocks_derivatives[block]) + for (const auto &[indices, id] : blocks_derivatives[block]) { + auto [eq, var, lag] = indices; int eqr = getBlockEquationID(block, eq); int varr = getBlockVariableID(block, var); if (lag == 0) @@ -782,8 +774,9 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const case BlockSimulationType::solveTwoBoundariesSimple: case BlockSimulationType::solveTwoBoundariesComplete: output << " else" << endl; - for (const auto &[eq, var, lag, id] : blocks_derivatives[block]) + for (const auto &[indices, id] : blocks_derivatives[block]) { + auto [eq, var, lag] = indices; int eqr = getBlockEquationID(block, eq); int varr = getBlockVariableID(block, var); ostringstream tmp_output; @@ -1185,8 +1178,8 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id file_open = true; } map, expr_t> tmp_block_endo_derivative; - for (const auto &it : blocks_derivatives[block]) - tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it); + for (const auto &[idx, d] : blocks_derivatives[block]) + tmp_block_endo_derivative[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d; map, expr_t> tmp_exo_derivative; for (const auto &it : derivative_exo[block]) tmp_exo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second; @@ -1393,8 +1386,9 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id case BlockSimulationType::solveTwoBoundariesComplete: case BlockSimulationType::solveTwoBoundariesSimple: count_u = feedback_variables.size(); - for (const auto &[eq, var, lag, ignore] : blocks_derivatives[block]) + for (const auto &[indices, ignore] : blocks_derivatives[block]) { + auto [eq, var, lag] = indices; int eqr = getBlockEquationID(block, eq); int varr = getBlockVariableID(block, var); if (eq >= block_recursive and var >= block_recursive) @@ -1419,7 +1413,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id Uf[eqr].Ufl->lag = lag; FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag); fnumexpr.write(code_file, instruction_number); - compileChainRuleDerivative(code_file, instruction_number, eqr, varr, lag, map_idx); + compileChainRuleDerivative(code_file, instruction_number, block, eq, var, lag, map_idx); FSTPU_ fstpu(count_u); fstpu.write(code_file, instruction_number); count_u++; @@ -1808,8 +1802,9 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, int num, int block_size = blocks[num].size; int block_mfs = blocks[num].mfs_size; int block_recursive = blocks[num].getRecursiveSize(); - for (const auto &[eq, var, lag, ignore] : blocks_derivatives[num]) + for (const auto &[indices, ignore] : blocks_derivatives[num]) { + auto [eq, var, lag] = indices; if (lag != 0 && !is_two_boundaries) continue; if (eq >= block_recursive && var >= block_recursive) @@ -3240,8 +3235,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de tmp_s.str(""); count_lead_lag_incidence = 0; dynamic_jacob_map_t reordered_dynamic_jacobian; - for (const auto &it : blocks_derivatives[block]) - reordered_dynamic_jacobian[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it); + for (const auto &[idx, d] : blocks_derivatives[block]) + reordered_dynamic_jacobian[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d; output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];" << endl; int last_var = -1; vector local_state_var; @@ -3415,10 +3410,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de if (block == 0) { set> row_state_var_incidence; - for (const auto &it : blocks_derivatives[block]) - if (auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(block, get<1>(it))+1); + for (const auto &[idx, ignore] : blocks_derivatives[block]) + if (auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(block, get<1>(idx))+1); it_state_var != state_var.end()) - if (auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(block, get<0>(it))+1); + if (auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(block, get<0>(idx))+1); 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(); @@ -4973,108 +4968,80 @@ DynamicModel::writeRevXrefs(ostream &output, const map, set> } } -map, int> -DynamicModel::get_Derivatives(int block) -{ - int max_lag, max_lead; - map, int> Derivatives; - BlockSimulationType simulation_type = blocks[block].simulation_type; - if (simulation_type == BlockSimulationType::evaluateBackward - || simulation_type == BlockSimulationType::evaluateForward) - { - max_lag = 1; - max_lead = 1; - blocks[block].max_lag = max_lag; - blocks[block].max_lead = max_lead; - } - else - { - max_lag = blocks[block].max_lag; - max_lead = blocks[block].max_lead; - } - int block_size = blocks[block].size; - int block_nb_recursive = blocks[block].getRecursiveSize(); - for (int lag = -max_lag; lag <= max_lead; lag++) - { - for (int eq = 0; eq < block_size; eq++) - { - int eqr = getBlockEquationID(block, eq); - for (int var = 0; var < block_size; var++) - { - int varr = getBlockVariableID(block, var); - if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end()) - { - bool OK = true; - if (auto its = Derivatives.find({ lag, eq, var, eqr, varr }); - its != Derivatives.end() && its->second == 2) - OK = false; - if (OK) - { - if (getBlockEquationType(block, eq) == EquationType::evaluate_s - && eq < block_nb_recursive) - //It's a normalized equation, we have to recompute the derivative using chain rule derivative function - Derivatives[{ lag, eq, var, eqr, varr }] = 1; - else - //It's a feedback equation we can use the derivatives - Derivatives[{ lag, eq, var, eqr, varr }] = 0; - } - if (var < block_nb_recursive) - { - int eqs = getBlockEquationID(block, var); - for (int vars = block_nb_recursive; vars < block_size; vars++) - { - int varrs = getBlockVariableID(block, vars); - //A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation) - if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end()) - Derivatives[{ lag, eq, vars, eqr, varrs }] = 2; - } - } - } - } - } - } - return Derivatives; +map, DynamicModel::BlockDerivativeType> +DynamicModel::determineBlockDerivativesType(int blk) +{ + map, BlockDerivativeType> derivType; + int size = blocks[blk].size; + int nb_recursive = blocks[blk].getRecursiveSize(); + for (int lag = -blocks[blk].max_lag; lag <= blocks[blk].max_lead; lag++) + for (int eq = 0; eq < size; eq++) + for (int var = 0; var < size; var++) + if (int eq_orig = getBlockEquationID(blk, eq), var_orig = getBlockVariableID(blk, var); + dynamic_jacobian.find({ lag, eq_orig, var_orig }) != dynamic_jacobian.end()) + { + if (getBlockEquationType(blk, eq) == EquationType::evaluate_s + && eq < nb_recursive) + /* It’s a normalized recursive equation, we have to recompute + the derivative using the chain rule */ + derivType[{ lag, eq, var }] = BlockDerivativeType::normalizedChainRule; + else if (derivType.find({ lag, eq, var }) == derivType.end()) + derivType[{ lag, eq, var }] = BlockDerivativeType::standard; + + if (var < nb_recursive) + for (int feedback_var = nb_recursive; feedback_var < size; feedback_var++) + if (derivType.find({ lag, var, feedback_var }) != derivType.end()) + /* A new derivative needs to be computed using the chain rule + (a feedback variable appears in the recursive equation + defining the current variable) */ + derivType[{ lag, eq, feedback_var }] = BlockDerivativeType::chainRule; + } + return derivType; } void DynamicModel::computeChainRuleJacobian() { - map recursive_variables; int nb_blocks = blocks.size(); blocks_derivatives.resize(nb_blocks); - for (int block = 0; block < nb_blocks; block++) + for (int blk = 0; blk < nb_blocks; blk++) { - block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives; - recursive_variables.clear(); - int block_nb_recursives = blocks[block].getRecursiveSize(); - for (int i = 0; i < block_nb_recursives; i++) + int nb_recursives = blocks[blk].getRecursiveSize(); + + // Create a map from recursive vars to their defining (normalized) equation + map recursive_vars; + for (int i = 0; i < nb_recursives; i++) { - if (getBlockEquationType(block, i) == EquationType::evaluate_s) - recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i); + int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(blk, i)), 0); + if (getBlockEquationType(blk, i) == EquationType::evaluate_s) + recursive_vars[deriv_id] = getBlockEquationRenormalizedExpr(blk, i); else - recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i); + recursive_vars[deriv_id] = getBlockEquationExpr(blk, i); } - auto Derivatives = get_Derivatives(block); - for (const auto &it : Derivatives) + + // Compute the block derivatives + for (const auto &[indices, derivType] : determineBlockDerivativesType(blk)) { - int Deriv_type = it.second; - auto [lag, eq, var, eqr, varr] = it.first; - if (Deriv_type == 0) - first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }]; - else if (Deriv_type == 1) - first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables); - else if (Deriv_type == 2) + auto [lag, eq, var] = indices; + int eq_orig = getBlockEquationID(blk, eq), var_orig = getBlockVariableID(blk, var); + int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, var_orig), lag); + expr_t d{nullptr}; + switch (derivType) { - if (getBlockEquationType(block, eq) == EquationType::evaluate_s - && eq < block_nb_recursives) - first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables); - else - first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables); + case BlockDerivativeType::standard: + d = derivatives[1][{ eq_orig, deriv_id }]; + break; + case BlockDerivativeType::chainRule: + d = equations[eq_orig]->getChainRuleDerivative(deriv_id, recursive_vars); + break; + case BlockDerivativeType::normalizedChainRule: + d = equation_type_and_normalized_equation[eq_orig].second->getChainRuleDerivative(deriv_id, recursive_vars); + break; } - tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]); + + blocks_derivatives[blk][{ eq, var, lag }] = d; } - blocks_derivatives[block] = tmp_derivatives; } } diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index af64ca2b..9cd29c13 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -128,15 +128,19 @@ private: void writeSetAuxiliaryVariables(const string &basename, bool julia) const; void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const; - //! Computes jacobian and prepares for equation normalization - /*! Using values from initval/endval blocks and parameter initializations: - - computes the jacobian for the model w.r. to contemporaneous variables - - removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff) - */ - //void evaluateJacobian(const eval_context_t &eval_context, jacob_map *j_m, bool dynamic); + // Used by determineBlockDerivativesType() + enum class BlockDerivativeType + { + standard, + chainRule, + normalizedChainRule + }; + + /* For each tuple (lag, eq, var) within the given block, determine the type + of the derivative to be computed. Indices are within the block (i.e. + between 0 and blocks[blk].size-1). */ + map, BlockDerivativeType> determineBlockDerivativesType(int blk); - //! return a map on the block jacobian - map, int> get_Derivatives(int block); //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables void computeChainRuleJacobian(); @@ -150,7 +154,7 @@ private: //! Write derivative code of an equation w.r. to a variable void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const; //! Write chain rule derivative code of an equation w.r. to a variable - void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, const map_idx_t &map_idx) const; + void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const map_idx_t &map_idx) const; //! Get the type corresponding to a derivation ID SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override; diff --git a/src/ModelTree.cc b/src/ModelTree.cc index a7bc507e..7eff6fef 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -103,17 +103,15 @@ ModelTree::copyHelper(const ModelTree &m) nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)}; for (const auto &it : m.dynamic_jacobian) dynamic_jacobian[it.first] = f(it.second); - for (const auto &it : m.first_chain_rule_derivatives) - first_chain_rule_derivatives[it.first] = f(it.second); for (const auto &it : m.equation_type_and_normalized_equation) equation_type_and_normalized_equation.emplace_back(it.first, f(it.second)); for (const auto &it : m.blocks_derivatives) { - block_derivatives_equation_variable_laglead_nodeid_t v; + map, expr_t> v; for (const auto &it2 : it) - v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2))); + v[it2.first] = f(it2.second); blocks_derivatives.push_back(v); } @@ -209,7 +207,6 @@ ModelTree::operator=(const ModelTree &m) endo_idx_block2orig = m.endo_idx_block2orig; eq_idx_orig2block = m.eq_idx_orig2block; endo_idx_orig2block = m.endo_idx_orig2block; - first_chain_rule_derivatives.clear(); map_idx = m.map_idx; equation_type_and_normalized_equation.clear(); blocks_derivatives.clear(); @@ -965,8 +962,12 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) // Merge the current block into the previous one blocks[i-1].size++; blocks[i-1].mfs_size = blocks[i-1].size; - blocks[i-1].max_lag = max(blocks[i-1].max_lag, max_lag); - blocks[i-1].max_lead = max(blocks[i-1].max_lead, max_lead); + /* For max lag/lead, the max of the two blocks is not enough. + We need to consider the case where a variable of the + previous block appears with a lag/lead in the current one + (the reverse case is excluded, by construction). */ + blocks[i-1].max_lag = is_lag ? 1 : max(blocks[i-1].max_lag, max_lag); + blocks[i-1].max_lead = is_lead ? 1 : max(blocks[i-1].max_lead, max_lead); blocks[i-1].n_static += blocks[i].n_static; blocks[i-1].n_forward += blocks[i].n_forward; blocks[i-1].n_backward += blocks[i].n_backward; @@ -1010,12 +1011,13 @@ ModelTree::determineLinearBlocks() { BlockSimulationType simulation_type = blocks[block].simulation_type; int block_size = blocks[block].size; - block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block]; + auto derivatives_block = blocks_derivatives[block]; int first_variable_position = blocks[block].first_equation; if (simulation_type == BlockSimulationType::solveBackwardComplete || simulation_type == BlockSimulationType::solveForwardComplete) - for (const auto &[ignore, ignore2, lag, d1] : derivatives_block) + for (const auto &[indices, d1] : derivatives_block) { + int lag = get<2>(indices); if (lag == 0) { set> endogenous; @@ -1031,8 +1033,9 @@ ModelTree::determineLinearBlocks() } else if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - for (const auto &[ignore, ignore2, lag, d1] : derivatives_block) + for (const auto &[indices, d1] : derivatives_block) { + int lag = get<2>(indices); set> endogenous; d1->collectEndogenous(endogenous); if (endogenous.size() > 0) diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 34862c3c..04875f98 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -55,12 +55,6 @@ using equation_type_and_normalized_equation_t = vector>; -//! for a block contains derivatives tuple -using block_derivatives_equation_variable_laglead_nodeid_t = vector>; - -//! for all blocks derivatives description -using blocks_derivatives_t = vector; - //! Shared code for static and dynamic models class ModelTree : public DataTree { @@ -178,17 +172,15 @@ protected: Set by updateReverseVariableEquationOrderings() */ vector eq_idx_orig2block, endo_idx_orig2block; - //! Store the derivatives or the chainrule derivatives:map, expr_t> - using first_chain_rule_derivatives_t = map, expr_t>; - first_chain_rule_derivatives_t first_chain_rule_derivatives; - map_idx_t map_idx; //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation equation_type_and_normalized_equation_t equation_type_and_normalized_equation; - //! for all blocks derivatives description - blocks_derivatives_t blocks_derivatives; + /* Stores derivatives of each block. + The tuple is: equation number (inside the block), variable number (inside + the block), lead/lag */ + vector, expr_t>> blocks_derivatives; //! Map the derivatives for a block tuple using derivative_t = map, expr_t>; diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 503ce8f3..e4244577 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -139,10 +139,10 @@ StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_nu } void -StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const +StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const { - if (auto it = first_chain_rule_derivatives.find({ eqr, varr, lag }); - it != first_chain_rule_derivatives.end()) + if (auto it = blocks_derivatives[blk].find({ eq, var, lag }); + it != blocks_derivatives[blk].end()) it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false); else { @@ -157,7 +157,6 @@ StaticModel::computeTemporaryTermsOrdered() map> first_occurence; map reference_count; BinaryOpNode *eq_node; - first_chain_rule_derivatives_t::const_iterator it_chr; ostringstream tmp_s; map_idx.clear(); @@ -192,11 +191,8 @@ StaticModel::computeTemporaryTermsOrdered() eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i); } } - for (const auto &it : blocks_derivatives[block]) - { - expr_t id = get<3>(it); - id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1); - } + for (const auto &[ignore, id] : blocks_derivatives[block]) + id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1); v_temporary_terms_inuse[block] = {}; computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]); } @@ -218,11 +214,8 @@ StaticModel::computeTemporaryTermsOrdered() eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i); } } - for (const auto &it : blocks_derivatives[block]) - { - expr_t id = get<3>(it); - id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); - } + for (const auto &[ignore, id] : blocks_derivatives[block]) + id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1); } for (int block = 0; block < nb_blocks; block++) @@ -241,11 +234,8 @@ StaticModel::computeTemporaryTermsOrdered() eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); } } - for (const auto &it : blocks_derivatives[block]) - { - expr_t id = get<3>(it); - id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); - } + for (const auto &[ignore, id] : blocks_derivatives[block]) + id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); for (int i = 0; i < blocks[block].size; i++) for (const auto &it : v_temporary_terms[block][i]) it->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); @@ -453,8 +443,9 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const case BlockSimulationType::solveForwardSimple: case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: - for (const auto &[eq, var, ignore, id] : blocks_derivatives[block]) + for (const auto &[indices, id] : blocks_derivatives[block]) { + auto &[eq, var, ignore] = indices; int eqr = getBlockEquationID(block, eq); int varr = getBlockVariableID(block, var); output << " g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = "; @@ -817,8 +808,9 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: count_u = feedback_variables.size(); - for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[block]) + for (const auto &[indices, ignore2] : blocks_derivatives[block]) { + auto [eq, var, ignore] = indices; int eqr = getBlockEquationID(block, eq); int varr = getBlockVariableID(block, var); if (eq >= block_recursive && var >= block_recursive) @@ -838,7 +830,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map Uf[eqr].Ufl->var = varr; FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr); fnumexpr.write(code_file, instruction_number); - compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx, temporary_terms); + compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, map_idx, temporary_terms); FSTPSU_ fstpsu(count_u); fstpsu.write(code_file, instruction_number); count_u++; @@ -1005,14 +997,15 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveForwardComplete: count_u = feedback_variables.size(); - for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[block]) + for (const auto &[indices, ignore2] : blocks_derivatives[block]) { + auto &[eq, var, ignore] = indices; int eqr = getBlockEquationID(block, eq); int varr = getBlockVariableID(block, var); FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0); fnumexpr.write(code_file, instruction_number); - compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2 /*temporary_terms*/); + compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, map_idx2[block], tt2 /*temporary_terms*/); FSTPG2_ fstpg2(eq, var); fstpg2.write(code_file, instruction_number); @@ -1055,8 +1048,9 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, int num, int block_size = blocks[num].size; int block_mfs = blocks[num].mfs_size; int block_recursive = blocks[num].getRecursiveSize(); - for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[num]) + for (const auto &[indices, ignore2] : blocks_derivatives[num]) { + auto [eq, var, ignore] = indices; int lag = 0; if (eq >= block_recursive && var >= block_recursive) { @@ -2104,121 +2098,43 @@ StaticModel::addAllParamDerivId(set &deriv_id_set) deriv_id_set.insert(i + symbol_table.endo_nbr()); } -map, int> -StaticModel::get_Derivatives(int block) -{ - map, int> Derivatives; - int block_size = blocks[block].size; - int block_nb_recursive = blocks[block].getRecursiveSize(); - int lag = 0; - for (int eq = 0; eq < block_size; eq++) - { - int eqr = getBlockEquationID(block, eq); - for (int var = 0; var < block_size; var++) - { - int varr = getBlockVariableID(block, var); - if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end()) - { - bool OK = true; - if (auto its = Derivatives.find({ lag, eq, var, eqr, varr }); - its != Derivatives.end() && its->second == 2) - OK = false; - - if (OK) - { - if (getBlockEquationType(block, eq) == EquationType::evaluate_s - && eq < block_nb_recursive) - //It's a normalized equation, we have to recompute the derivative using chain rule derivative function - Derivatives[{ lag, eq, var, eqr, varr }] = 1; - else - //It's a feedback equation we can use the derivatives - Derivatives[{ lag, eq, var, eqr, varr }] = 0; - } - if (var < block_nb_recursive) - { - int eqs = getBlockEquationID(block, var); - for (int vars = block_nb_recursive; vars < block_size; vars++) - { - int varrs = getBlockVariableID(block, vars); - //A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation) - if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end()) - Derivatives[{ lag, eq, vars, eqr, varrs }] = 2; - } - } - } - } - } - - return Derivatives; -} - void StaticModel::computeChainRuleJacobian() { - map recursive_variables; int nb_blocks = blocks.size(); blocks_derivatives.resize(nb_blocks); - for (int block = 0; block < nb_blocks; block++) + for (int blk = 0; blk < nb_blocks; blk++) { - block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives; - recursive_variables.clear(); - BlockSimulationType simulation_type = blocks[block].simulation_type; - int block_size = blocks[block].size; - int block_nb_recursives = blocks[block].getRecursiveSize(); - if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) + int nb_recursives = blocks[blk].getRecursiveSize(); + + map recursive_vars; + for (int i = 0; i < nb_recursives; i++) { - for (int i = 0; i < block_nb_recursives; i++) + int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(blk, i)), 0); + if (getBlockEquationType(blk, i) == EquationType::evaluate_s) + recursive_vars[deriv_id] = getBlockEquationRenormalizedExpr(blk, i); + else + recursive_vars[deriv_id] = getBlockEquationExpr(blk, i); + } + + assert(blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesSimple + && blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesComplete); + + int size = blocks[blk].size; + + for (int eq = nb_recursives; eq < size; eq++) + { + int eq_orig = getBlockEquationID(blk, eq); + for (int var = nb_recursives; var < size; var++) { - if (getBlockEquationType(block, i) == EquationType::evaluate_s) - recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i); - else - recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i); - } - auto Derivatives = get_Derivatives(block); - for (const auto &it : Derivatives) - { - int Deriv_type = it.second; - auto [lag, eq, var, eqr, varr] = it.first; - if (Deriv_type == 0) - first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }]; - else if (Deriv_type == 1) - first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables); - else if (Deriv_type == 2) - { - if (getBlockEquationType(block, eq) == EquationType::evaluate_s - && eq < block_nb_recursives) - first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables); - else - first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables); - } - tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]); + int var_orig = getBlockVariableID(blk, var); + expr_t d1 = equations[eq_orig]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, var_orig), 0), recursive_vars); + if (d1 == Zero) + continue; + + blocks_derivatives[blk][{ eq, var, 0 }] = d1; } } - else - { - for (int i = 0; i < block_nb_recursives; i++) - { - if (getBlockEquationType(block, i) == EquationType::evaluate_s) - recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i); - else - recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i); - } - for (int eq = block_nb_recursives; eq < block_size; eq++) - { - int eqr = getBlockEquationID(block, eq); - for (int var = block_nb_recursives; var < block_size; var++) - { - int varr = getBlockVariableID(block, var); - expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), 0), recursive_variables); - if (d1 == Zero) - continue; - first_chain_rule_derivatives[{ eqr, varr, 0 }] = d1; - tmp_derivatives.emplace_back(eq, var, 0, first_chain_rule_derivatives[{ eqr, varr, 0 }]); - } - } - } - blocks_derivatives[block] = tmp_derivatives; } } diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 8414e68d..1a93f313 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -77,7 +77,7 @@ private: //! Write derivative code of an equation w.r. to a variable void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const; //! Write chain rule derivative code of an equation w.r. to a variable - void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const; + void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const; //! Get the type corresponding to a derivation ID SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override; @@ -87,8 +87,6 @@ private: int getSymbIDByDerivID(int deriv_id) const noexcept(false) override; //! Compute the column indices of the static Jacobian void computeStatJacobianCols(); - //! return a map on the block jacobian - map, int> get_Derivatives(int block); //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables void computeChainRuleJacobian();