From 39407083bedeafa45a0a86633f743488d335f457 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 28 Apr 2020 18:00:26 +0200 Subject: [PATCH] Block decomposition: ensure that the cutoff option only impacts normalization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously, the cutoff option would also impact the block decomposition itself, since it would had an influence on the incidence matrix used for computing the blocks and their derivatives. The problem is that, in the general case, it’s quite possible that an element of the numerical Jacobian be zero at the evaluation point, while being quite different from zero along the simulation path. A typical example is an expression of the type x*y, where y is an endogenous and x is an exogenous not present in the initval block (and hence initialized to zero). --- src/DynamicModel.cc | 49 +++++++++-------- src/ModelTree.cc | 129 +++++++++++++------------------------------- src/ModelTree.hh | 21 +++----- src/StaticModel.cc | 4 +- 4 files changed, 72 insertions(+), 131 deletions(-) diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 90451ed2..1525b30a 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -3228,7 +3228,7 @@ 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; + map, expr_t> reordered_dynamic_jacobian; 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; @@ -4812,7 +4812,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian); - computePrologueAndEpilogue(static_jacobian); + computePrologueAndEpilogue(); auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous(); @@ -4820,7 +4820,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO cout << "Finding the optimal block decomposition of the model ..." << endl; - auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false); + auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(); reduceBlocksAndTypeDetermination(linear_decomposition); @@ -4971,26 +4971,31 @@ DynamicModel::determineBlockDerivativesType(int blk) 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; + { + set> endos_and_lags; + int eq_orig = getBlockEquationID(blk, eq); + equations[eq_orig]->collectEndogenous(endos_and_lags); + for (int var = 0; var < size; var++) + if (int var_orig = getBlockVariableID(blk, var); + endos_and_lags.find({ var_orig, lag }) != endos_and_lags.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; - } + 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; } diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 7c69901e..f92e4db9 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -101,8 +101,6 @@ ModelTree::copyHelper(const ModelTree &m) trend_symbols_map[it.first] = f(it.second); for (const auto &it : m.nonstationary_symbols_map) 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.equation_type_and_normalized_equation) equation_type_and_normalized_equation.emplace_back(it.first, f(it.second)); @@ -200,7 +198,6 @@ ModelTree::operator=(const ModelTree &m) trend_symbols_map.clear(); nonstationary_symbols_map.clear(); - dynamic_jacobian.clear(); eq_idx_block2orig = m.eq_idx_block2orig; endo_idx_block2orig = m.endo_idx_block2orig; eq_idx_orig2block = m.eq_idx_orig2block; @@ -285,7 +282,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo } void -ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian) +ModelTree::computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, double cutoff, const jacob_map_t &static_jacobian) { cout << "Normalizing the model..." << endl; @@ -338,33 +335,6 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian TODO: Explain why symbolic_jacobian is not contemporaneous. */ auto symbolic_jacobian = computeSymbolicJacobian(); found_normalization = computeNormalization(symbolic_jacobian, true); - if (found_normalization) - { - /* Update the Jacobian matrices by ensuring that they have a - numerical value associated to each symbolic occurrence. - TODO: Explain why this is needed. Maybe in order to have the - right incidence matrix in computePrologueAndEpilogue? */ - for (const auto &[eq_and_endo, ignore] : symbolic_jacobian) - { - if (static_jacobian.find(eq_and_endo) == static_jacobian.end()) - static_jacobian[eq_and_endo] = 0; - if (dynamic_jacobian.find({ 0, eq_and_endo.first, eq_and_endo.second }) == dynamic_jacobian.end()) - dynamic_jacobian[{ 0, eq_and_endo.first, eq_and_endo.second }] = Zero; - if (contemporaneous_jacobian.find(eq_and_endo) == contemporaneous_jacobian.end()) - contemporaneous_jacobian[eq_and_endo] = 0; - try - { - if (derivatives[1].find({ eq_and_endo.first, getDerivID(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second), 0) }) == derivatives[1].end()) - derivatives[1][{ eq_and_endo.first, getDerivID(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second), 0) }] = Zero; - } - catch (DataTree::UnknownDerivIDException &e) - { - cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second)) - << " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl; - exit(EXIT_FAILURE); - } - } - } } if (!found_normalization) @@ -375,7 +345,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian } pair -ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) +ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) const { jacob_map_t contemporaneous_jacobian, static_jacobian; int nb_elements_contemporaneous_jacobian = 0; @@ -423,17 +393,10 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double static_jacobian[{ eq, var }] += val; else static_jacobian[{ eq, var }] = val; - - dynamic_jacobian[{ lag, eq, var }] = d1; } } } - // Get rid of the elements of the Jacobian matrix below the cutoff - // TODO: Why? - for (const auto &it : jacobian_elements_to_delete) - derivatives[1].erase(it); - if (jacobian_elements_to_delete.size() > 0) { cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded." << endl @@ -514,7 +477,7 @@ ModelTree::computeNaturalNormalization() } void -ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian) +ModelTree::computePrologueAndEpilogue() { const int n = equations.size(); @@ -533,17 +496,13 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian) “eq_idx_block2orig” (resp. “endo_idx_block2orig”). Stored in row-major order. */ vector IM(n*n, false); - if (cutoff == 0) - for (int i = 0; i < n; i++) - { - set> endos_and_lags; - equations[i]->collectEndogenous(endos_and_lags); - for (const auto &[endo, lag] : endos_and_lags) - IM[i * n + endo2eq[endo]] = true; - } - else - for (const auto &[eq_and_endo, val] : static_jacobian) - IM[eq_and_endo.first * n + endo2eq[eq_and_endo.second]] = true; + for (int i = 0; i < n; i++) + { + set> endos_and_lags; + equations[i]->collectEndogenous(endos_and_lags); + for (auto [endo, lag] : endos_and_lags) + IM[i * n + endo2eq[endo]] = true; + } bool something_has_been_done; // Find the prologue equations and place first the AR(1) shock equations first @@ -677,22 +636,24 @@ ModelTree::getVariableLeadLagByBlock(const vector &endo2simblock) const }; lag_lead_vector_t variable_lag_lead(nb_endo, { 0, 0 }), equation_lag_lead(nb_endo, { 0, 0 }); - for (const auto &[key, value] : dynamic_jacobian) + for (int eq = 0; eq < nb_endo; eq++) { - auto [lag, eq, endo] = key; - if (belong_to_same_block(endo, eq)) - { - variable_lag_lead[endo].first = max(variable_lag_lead[endo].first, -lag); - variable_lag_lead[endo].second = max(variable_lag_lead[endo].second, lag); - equation_lag_lead[eq].first = max(equation_lag_lead[eq].first, -lag); - equation_lag_lead[eq].second = max(equation_lag_lead[eq].second, lag); - } + set> endos_and_lags; + equations[eq]->collectEndogenous(endos_and_lags); + for (auto [endo, lag] : endos_and_lags) + if (belong_to_same_block(endo, eq)) + { + variable_lag_lead[endo].first = max(variable_lag_lead[endo].first, -lag); + variable_lag_lead[endo].second = max(variable_lag_lead[endo].second, lag); + equation_lag_lead[eq].first = max(equation_lag_lead[eq].first, -lag); + equation_lag_lead[eq].second = max(equation_lag_lead[eq].second, lag); + } } return { equation_lag_lead, variable_lag_lead }; } lag_lead_vector_t -ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_) +ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock() { int nb_var = symbol_table.endo_nbr(); int nb_simvars = nb_var - prologue - epilogue; @@ -700,11 +661,10 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob /* Construct the graph representing the dependencies between all variables that do not belong to the prologue or the epilogue. - For detecting dependencies between variables, use the static jacobian, - except when the cutoff is zero, in which case use the symbolic adjacency + For detecting dependencies between variables, use the symbolic adjacency matrix */ VariableDependencyGraph G(nb_simvars); - for (const auto &[key, value] : cutoff == 0 ? computeSymbolicJacobian() : static_jacobian) + for (const auto &[key, value] : computeSymbolicJacobian()) { auto [eq, endo] = key; if (eq_idx_orig2block[eq] >= prologue @@ -869,29 +829,14 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) set> endos_and_lags; equations[eq_idx_block2orig[eq]]->collectEndogenous(endos_and_lags); for (const auto &[endo, lag] : endos_and_lags) - { - if (linear_decomposition) - { - if (dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo }) - != dynamic_jacobian.end()) - { - max_lead = max(lag, max_lead); - max_lag = max(-lag, max_lag); - } - } - else - { - if (find(endo_idx_block2orig.begin()+first_eq, - endo_idx_block2orig.begin()+first_eq+blocks[i].size, endo) - != endo_idx_block2orig.begin()+first_eq+blocks[i].size - && dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo }) - != dynamic_jacobian.end()) - { - max_lead = max(lag, max_lead); - max_lag = max(-lag, max_lag); - } - } - } + if (linear_decomposition || + find(endo_idx_block2orig.begin()+first_eq, + endo_idx_block2orig.begin()+first_eq+blocks[i].size, endo) + != endo_idx_block2orig.begin()+first_eq+blocks[i].size) + { + max_lead = max(lag, max_lead); + max_lag = max(-lag, max_lag); + } } // Determine the block type @@ -940,12 +885,10 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) for (int j = blocks[i-1].first_equation; j < blocks[i-1].first_equation+blocks[i-1].size; j++) { - if (dynamic_jacobian.find({ -1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] }) - != dynamic_jacobian.end()) - is_lag = true; - if (dynamic_jacobian.find({ +1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] }) - != dynamic_jacobian.end()) - is_lead = true; + set> endos_and_lags; + equations[eq_idx_block2orig[first_eq]]->collectEndogenous(endos_and_lags); + is_lag = endos_and_lags.find({ endo_idx_block2orig[j], -1 }) != endos_and_lags.end(); + is_lead = endos_and_lags.find({ endo_idx_block2orig[j], 1 }) != endos_and_lags.end(); } BlockSimulationType prev_Type = blocks[i-1].simulation_type; diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 7166e3cf..2bffcf81 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -155,13 +155,6 @@ protected: //! Nonstationary variables and their deflators nonstationary_symbols_map_t nonstationary_symbols_map; - //! Sparse matrix of double to store the values of the Jacobian - /*! First index is lag, second index is equation number, third index is endogenous type specific ID */ - using dynamic_jacob_map_t = map, expr_t>; - - //! The jacobian without the elements below the cutoff - dynamic_jacob_map_t dynamic_jacobian; - /* Maps indices of equations in the block-decomposition order into original equation IDs */ vector eq_idx_block2orig; @@ -255,7 +248,7 @@ protected: //! Writes LaTeX model file void writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const; - //! Sparse matrix of double to store the values of the Jacobian + //! Sparse matrix of double to store the values of the static Jacobian /*! First index is equation number, second index is endogenous type specific ID */ using jacob_map_t = map, double>; @@ -287,24 +280,24 @@ protected: If no matching is found with a zero cutoff, an error message is printed. The resulting normalization is stored in endo2eq. */ - void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian); + void computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, double cutoff, const jacob_map_t &static_jacobian); //! Try to find a natural normalization if all equations are matched to an endogenous variable on the LHS bool computeNaturalNormalization(); //! Evaluate the jacobian (w.r.t. endogenous) and suppress all the elements below the cutoff - /*! Returns a pair (contemporaneous_jacobian, static_jacobian). Also fills - dynamic_jacobian. Elements below the cutoff are discarded. External functions are evaluated to 1. */ - pair evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose); + /*! Returns a pair (contemporaneous_jacobian, static_jacobian). + Elements below the cutoff are discarded. External functions are evaluated to 1. */ + pair evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) const; //! Select and reorder the non linear equations of the model void select_non_linear_equations_and_variables(); //! Search the equations and variables belonging to the prologue and the epilogue of the model - void computePrologueAndEpilogue(const jacob_map_t &static_jacobian); + void computePrologueAndEpilogue(); //! Determine the type of each equation of model and try to normalize the unnormalized equation void equationTypeDetermination(const map, expr_t> &first_order_endo_derivatives, int mfs); /* Compute the block decomposition and for a non-recusive block find the minimum feedback set Initializes the “blocks” structure, and fills the following fields: size, mfs_size, n_static, n_forward, n_backward, n_mixed. */ - lag_lead_vector_t computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_); + lag_lead_vector_t computeBlockDecompositionAndFeedbackVariablesForEachBlock(); /* Reduce the number of block by merging the same type of equations in the prologue and the epilogue, and determine the type of each block. diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 3e76cddf..e7c84cad 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -1125,7 +1125,7 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian); - computePrologueAndEpilogue(static_jacobian); + computePrologueAndEpilogue(); auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous(); @@ -1133,7 +1133,7 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co cout << "Finding the optimal block decomposition of the model ..." << endl; - auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false); + auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(); reduceBlocksAndTypeDetermination(false);