diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 90451ed2ca606cc9c7cf284a7c869c261b7ed9a6..1525b30ab301ae83b2f97e1a0168f06b0b2b32fd 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<tuple<int, int, int>, 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;
-
- 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;
- }
+ {
+ set<pair<int, int>> 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;
+ }
+ }
return derivType;
}
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 7c69901e69cc111185259012f01bf11ad45fb97b..f92e4db9ca6153caad1b696deae32df94b37adad 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::jacob_map_t, ModelTree::jacob_map_t>
-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<bool> IM(n*n, false);
- if (cutoff == 0)
- for (int i = 0; i < n; i++)
- {
- set<pair<int, int>> 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<pair<int, int>> 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<int> &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<pair<int, int>> 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<pair<int, int>> 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<pair<int, int>> 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 7166e3cf35d65e5d159a4e9a6ad15fbbf272545d..2bffcf811748304040b9e27f1c97040c6aa80234 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<tuple<int, int, int>, 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<int> 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<pair<int, int>, 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<jacob_map_t, jacob_map_t> 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<jacob_map_t, jacob_map_t> 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<tuple<int, int, int>, 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 3e76cddf14992d974356960121f4687854e6d559..e7c84cad8c3c8dde43f37de6aafcb7ae053d9b5b 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);