diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 2cd7fdd5976b39da14c2adca56598adeb53d18d0..a34e356cdec38a35bd49533e2deb1f6bf9816ad3 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3293,34 +3293,16 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
if (paramsDerivsOrder > 0 && !no_tmp_terms)
computeParamsDerivativesTemporaryTerms();
- if (!computingPassBlock(eval_context, no_tmp_terms) && block)
+ computingPassBlock(eval_context, no_tmp_terms);
+ if (block_decomposed)
+ computeBlockDynJacobianCols();
+ if (!block_decomposed && block)
{
cerr << "ERROR: Block decomposition requested but failed." << endl;
exit(EXIT_FAILURE);
}
}
-bool
-DynamicModel::computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms)
-{
- auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
- if (!computeNonSingularNormalization(contemporaneous_jacobian))
- return false;
- auto [prologue, epilogue] = computePrologueAndEpilogue();
- auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
- equationTypeDetermination(first_order_endo_derivatives, mfs);
- cout << "Finding the optimal block decomposition of the " << modelClassName() << "..." << endl;
- computeBlockDecomposition(prologue, epilogue);
- reduceBlockDecomposition();
- printBlockDecomposition();
- computeChainRuleJacobian();
- determineLinearBlocks();
- computeBlockDynJacobianCols();
- if (!no_tmp_terms)
- computeBlockTemporaryTerms();
- return true;
-}
-
void
DynamicModel::computeXrefs()
{
@@ -3462,9 +3444,11 @@ DynamicModel::determineBlockDerivativesType(int blk)
void
DynamicModel::computeChainRuleJacobian()
{
- int nb_blocks = blocks.size();
+ size_t nb_blocks { blocks.size() };
+
blocks_derivatives.resize(nb_blocks);
- for (int blk = 0; blk < nb_blocks; blk++)
+
+ for (int blk {0}; blk < static_cast<int>(nb_blocks); blk++)
{
int nb_recursives = blocks[blk].getRecursiveSize();
@@ -3507,52 +3491,71 @@ DynamicModel::computeChainRuleJacobian()
blocks_derivatives[blk][{ eq, var, lag }] = d;
}
}
-}
-void
-DynamicModel::computeBlockDynJacobianCols()
-{
- int nb_blocks = blocks.size();
+ /* Also store information and derivatives w.r.t. other types of variables
+ (for the stochastic mode) */
blocks_derivatives_other_endo.resize(nb_blocks);
blocks_derivatives_exo.resize(nb_blocks);
blocks_derivatives_exo_det.resize(nb_blocks);
blocks_other_endo.resize(nb_blocks);
blocks_exo.resize(nb_blocks);
blocks_exo_det.resize(nb_blocks);
+ for (auto &[indices, d1] : derivatives[1])
+ {
+ auto [eq_orig, deriv_id] { vectorToTuple<2>(indices) };
+ int block_eq { eq2block[eq_orig] };
+ int eq { getBlockInitialEquationID(block_eq, eq_orig) };
+ int var { getTypeSpecificIDByDerivID(deriv_id) };
+ int lag { getLagByDerivID(deriv_id) };
+ switch (getTypeByDerivID(indices[1]))
+ {
+ case SymbolType::endogenous:
+ if (block_eq != endo2block[var])
+ {
+ blocks_derivatives_other_endo[block_eq][{ eq, var, lag }] = d1;
+ blocks_other_endo[block_eq].insert(var);
+ }
+ break;
+ case SymbolType::exogenous:
+ blocks_derivatives_exo[block_eq][{ eq, var, lag }] = d1;
+ blocks_exo[block_eq].insert(var);
+ break;
+ case SymbolType::exogenousDet:
+ blocks_derivatives_exo_det[block_eq][{ eq, var, lag }] = d1;
+ blocks_exo_det[block_eq].insert(var);
+ break;
+ default:
+ break;
+ }
+ }
+}
+
+void
+DynamicModel::computeBlockDynJacobianCols()
+{
+ size_t nb_blocks { blocks.size() };
// Structures used for lexicographic ordering over (lag, var ID)
vector<set<pair<int, int>>> dynamic_endo(nb_blocks), dynamic_other_endo(nb_blocks),
dynamic_exo(nb_blocks), dynamic_exo_det(nb_blocks);
for (auto & [indices, d1] : derivatives[1])
{
- int eq_orig = indices[0];
- int block_eq = eq2block[eq_orig];
- int eq = getBlockInitialEquationID(block_eq, eq_orig);
- int var { getTypeSpecificIDByDerivID(indices[1]) };
- int lag = getLagByDerivID(indices[1]);
- switch (getTypeByDerivID(indices[1]))
+ auto [eq_orig, deriv_id] { vectorToTuple<2>(indices) };
+ int block_eq { eq2block[eq_orig] };
+ int var { getTypeSpecificIDByDerivID(deriv_id) };
+ int lag { getLagByDerivID(deriv_id) };
+ switch (getTypeByDerivID(deriv_id))
{
case SymbolType::endogenous:
if (block_eq == endo2block[var])
- {
- int var_in_block = getBlockInitialVariableID(block_eq, var);
- dynamic_endo[block_eq].emplace(lag, var_in_block);
- }
+ dynamic_endo[block_eq].emplace(lag, getBlockInitialVariableID(block_eq, var));
else
- {
- blocks_derivatives_other_endo[block_eq][{ eq, var, lag }] = derivatives[1][{ eq_orig, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
- blocks_other_endo[block_eq].insert(var);
- dynamic_other_endo[block_eq].emplace(lag, var);
- }
+ dynamic_other_endo[block_eq].emplace(lag, var);
break;
case SymbolType::exogenous:
- blocks_derivatives_exo[block_eq][{ eq, var, lag }] = derivatives[1][{ eq_orig, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
- blocks_exo[block_eq].insert(var);
dynamic_exo[block_eq].emplace(lag, var);
break;
case SymbolType::exogenousDet:
- blocks_derivatives_exo_det[block_eq][{ eq, var, lag }] = derivatives[1][{ eq_orig, getDerivID(symbol_table.getID(SymbolType::exogenousDet, var), lag) }];
- blocks_exo_det[block_eq].insert(var);
dynamic_exo_det[block_eq].emplace(lag, var);
break;
default:
@@ -3565,7 +3568,7 @@ DynamicModel::computeBlockDynJacobianCols()
blocks_jacob_cols_other_endo.resize(nb_blocks);
blocks_jacob_cols_exo.resize(nb_blocks);
blocks_jacob_cols_exo_det.resize(nb_blocks);
- for (int blk = 0; blk < nb_blocks; blk++)
+ for (size_t blk {0}; blk < nb_blocks; blk++)
{
for (int index{0};
auto [lag, var] : dynamic_endo[blk])
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index b6a3fe780966cc14f6fe26ac02cd67c70b5c5e46..9e0c1780e0ce4315fa975c79a69361de41136895 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -169,8 +169,7 @@ private:
between 0 and blocks[blk].size-1). */
map<tuple<int, int, int>, BlockDerivativeType> determineBlockDerivativesType(int blk);
- //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
- void computeChainRuleJacobian();
+ void computeChainRuleJacobian() override;
string reform(const string &name) const;
@@ -299,12 +298,6 @@ private:
return blocks_jacob_cols_endo[blk].at({ var, lag });
}
- /* Compute block decomposition, its derivatives and temporary terms.
- Meant to be overriden in derived classes which don’t support block
- decomposition (currently Epilogue). Returns a boolean indicating success
- (failure can happen in normalization). */
- virtual bool computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms);
-
protected:
string
modelClassName() const override
diff --git a/src/ModelEquationBlock.cc b/src/ModelEquationBlock.cc
index 58db1c8a1aacdc205eae895f47c016a8b44529ab..dfe4ef54d9cf1dd6a4aa4ada3a2469b35c7d93b5 100644
--- a/src/ModelEquationBlock.cc
+++ b/src/ModelEquationBlock.cc
@@ -30,12 +30,11 @@ PlannerObjective::PlannerObjective(SymbolTable &symbol_table_arg,
{
}
-bool
+void
PlannerObjective::computingPassBlock([[maybe_unused]] const eval_context_t &eval_context,
[[maybe_unused]] bool no_tmp_terms)
{
// Disable block decomposition on planner objective
- return false;
}
OrigRamseyDynamicModel::OrigRamseyDynamicModel(SymbolTable &symbol_table_arg,
@@ -521,10 +520,9 @@ Epilogue::writeOutput(ostream &output) const
SymbolList{move(symbol_list)}.writeOutput("M_.epilogue_var_list_", output);
}
-bool
+void
Epilogue::computingPassBlock([[maybe_unused]] const eval_context_t &eval_context,
[[maybe_unused]] bool no_tmp_terms)
{
// Disable block decomposition on epilogue blocks
- return false;
}
diff --git a/src/ModelEquationBlock.hh b/src/ModelEquationBlock.hh
index a34fb6dcfcf887dab547747823c225de00f210ec..8efe6aaf7aa24f98132b7e7afe7416bd04ece960 100644
--- a/src/ModelEquationBlock.hh
+++ b/src/ModelEquationBlock.hh
@@ -40,7 +40,7 @@ protected:
}
private:
- bool computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms) override;
+ void computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms) override;
};
class OrigRamseyDynamicModel : public DynamicModel
@@ -141,7 +141,7 @@ private:
//! Helper for public writeEpilogueFile
void writeStaticEpilogueFile(const string &basename) const;
void writeDynamicEpilogueFile(const string &basename) const;
- bool computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms) override;
+ void computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms) override;
};
#endif
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 364a1777cabab273b98c1d2e7ce6212449a9c5f3..e375c68468d0a5548536769a1de0cfbbfd3fe1fd 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1874,3 +1874,23 @@ ModelTree::joinMEXCompilationThreads()
for (auto &it : mex_compilation_threads)
it.join();
}
+
+void
+ModelTree::computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms)
+{
+ auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
+ if (!computeNonSingularNormalization(contemporaneous_jacobian))
+ return;
+ auto [prologue, epilogue] = computePrologueAndEpilogue();
+ auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
+ equationTypeDetermination(first_order_endo_derivatives, mfs);
+ cout << "Finding the optimal block decomposition of the " << modelClassName() << "..." << endl;
+ computeBlockDecomposition(prologue, epilogue);
+ reduceBlockDecomposition();
+ printBlockDecomposition();
+ computeChainRuleJacobian();
+ determineLinearBlocks();
+ if (!no_tmp_terms)
+ computeBlockTemporaryTerms();
+ block_decomposed = true;
+}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 87ac681cdd727b7efba8d52a8997db5beb2a018b..68887054021f26a37a00eff6bc707190d1fa1baa 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -189,6 +189,9 @@ protected:
};
};
+ // Whether block decomposition has been successfully computed
+ bool block_decomposed {false};
+
// Stores various informations on the blocks
vector<BlockInfo> blocks;
@@ -321,6 +324,7 @@ protected:
//! Writes LaTeX model file
void writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const;
+private:
//! 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>;
@@ -394,6 +398,7 @@ protected:
//! Determine for each block if it is linear or not
void determineLinearBlocks();
+protected:
//! Return the type of equation belonging to the block
EquationType
getBlockEquationType(int blk, int eq) const
@@ -444,11 +449,23 @@ protected:
};
//! Initialize equation_reordered & variable_reordered
void initializeVariablesAndEquations();
+
+private:
//! Returns the 1st derivatives w.r.t. endogenous in a different format
/*! Returns a map (equation, type-specific ID, lag) → derivative.
Assumes that derivatives have already been computed. */
map<tuple<int, int, int>, expr_t> collectFirstOrderDerivativesEndogenous();
+protected:
+ //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
+ virtual void computeChainRuleJacobian() = 0;
+
+ /* Compute block decomposition, its derivatives and temporary terms. Meant to
+ be overriden in derived classes which don’t support block decomposition
+ (currently Epilogue and PlannerObjective). Sets “block_decomposed” to true
+ in case of success. */
+ virtual void computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms);
+
/* Get column number within Jacobian of a given block.
“var” is the block-specific endogenous variable index. */
virtual int getBlockJacobianEndoCol(int blk, int var, int lag) const = 0;
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 33ed5c7ea850be04de1920a5d96ed830f46cb004..cc02959535ce7455b730e8491032b5c1ae1fa8f1 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -375,33 +375,14 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
if (paramsDerivsOrder > 0 && !no_tmp_terms)
computeParamsDerivativesTemporaryTerms();
- if (!computingPassBlock(eval_context, no_tmp_terms) && block)
+ computingPassBlock(eval_context, no_tmp_terms);
+ if (!block_decomposed && block)
{
cerr << "ERROR: Block decomposition requested but failed. If your model does not have a steady state, you may want to try the 'no_static' option of the 'model' block." << endl;
exit(EXIT_FAILURE);
}
}
-bool
-StaticModel::computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms)
-{
- auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
- if (!computeNonSingularNormalization(contemporaneous_jacobian))
- return false;
- auto [prologue, epilogue] = computePrologueAndEpilogue();
- auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
- equationTypeDetermination(first_order_endo_derivatives, mfs);
- cout << "Finding the optimal block decomposition of the " << modelClassName() << "..." << endl;
- computeBlockDecomposition(prologue, epilogue);
- reduceBlockDecomposition();
- printBlockDecomposition();
- computeChainRuleJacobian();
- determineLinearBlocks();
- if (!no_tmp_terms)
- computeBlockTemporaryTerms();
- return true;
-}
-
void
StaticModel::writeStaticMFile(const string &basename) const
{
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 471a309e834000a3511cc7697f0453baee690b1c..d76202d10e700d4e964528971b8fccf9d8b9be37 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -91,8 +91,7 @@ private:
return symbol_table.endo_nbr();
}
- //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
- void computeChainRuleJacobian();
+ void computeChainRuleJacobian() override;
// Helper for writing MATLAB/Octave functions for residuals/derivatives and their temporary terms
void writeStaticMFileHelper(const string &basename,
@@ -113,12 +112,6 @@ private:
return var;
}
- /* Compute block decomposition, its derivatives and temporary terms.
- Meant to be overriden in derived classes which don’t support block
- decomposition (currently PlannerObjective). Returns a boolean indicating success
- (failure can happen in normalization). */
- virtual bool computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms);
-
// Write the block structure of the model in the driver file
void writeBlockDriverOutput(ostream &output) const;