diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index c1cd8d9e7cefad8bbfa80a3a62e2c2d9e2dcab40..a106ae19f3135b277452c5004eeacb643f709772 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -2453,7 +2453,7 @@ ModelComparisonStatement::writeJsonOutput(ostream &output) const
output << "}";
}
-PlannerObjectiveStatement::PlannerObjectiveStatement(const StaticModel &model_tree_arg) :
+PlannerObjectiveStatement::PlannerObjectiveStatement(const PlannerObjective &model_tree_arg) :
model_tree{model_tree_arg}
{
}
@@ -2473,7 +2473,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct,
mod_file_struct.planner_objective_present = true;
}
-const StaticModel &
+const PlannerObjective &
PlannerObjectiveStatement::getPlannerObjective() const
{
return model_tree;
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 9888fbc114f0e62d1dcbb52af4063e48b5c3e41e..1cb72657dabe589a52489f65509b0e2a74b07d18 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -580,10 +580,10 @@ public:
class PlannerObjectiveStatement : public Statement
{
private:
- StaticModel model_tree;
+ PlannerObjective model_tree;
bool computing_pass_called{false};
public:
- explicit PlannerObjectiveStatement(const StaticModel &model_tree_arg);
+ explicit PlannerObjectiveStatement(const PlannerObjective &model_tree_arg);
/*! \todo check there are only endogenous variables at the current period in the objective
(no exogenous, no lead/lag) */
void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
@@ -592,7 +592,7 @@ public:
void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
void writeJsonOutput(ostream &output) const override;
//! Return a reference the Planner Objective model tree
- const StaticModel &getPlannerObjective() const;
+ const PlannerObjective &getPlannerObjective() const;
};
class BVARDensityStatement : public Statement
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 2d6a75358afa9bfc6fee031caef41601e4a4537f..85037d4dcce09c630db097418a98f47e8582b8f5 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -575,13 +575,17 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
};
int jacobian_ncols_exo {symbol_table.exo_nbr()};
int jacobian_ncols_exo_det {symbol_table.exo_det_nbr()};
+ vector<int> eq_idx(equations.size());
+ iota(eq_idx.begin(), eq_idx.end(), 0);
+ vector<int> endo_idx(symbol_table.endo_nbr());
+ iota(endo_idx.begin(), endo_idx.end(), 0);
code_file << FBEGINBLOCK_{symbol_table.endo_nbr(),
simulation_type,
0,
symbol_table.endo_nbr(),
- endo_idx_block2orig,
- eq_idx_block2orig,
+ endo_idx,
+ eq_idx,
false,
symbol_table.endo_nbr(),
max_endo_lag,
@@ -3285,47 +3289,41 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeParamsDerivatives(paramsDerivsOrder);
}
+ computeTemporaryTerms(!use_dll, no_tmp_terms);
- if (block)
- {
- auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
-
- computeNonSingularNormalization(contemporaneous_jacobian, true);
-
- auto [prologue, epilogue] = computePrologueAndEpilogue();
-
- auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
-
- equationTypeDetermination(first_order_endo_derivatives, mfs);
-
- cout << "Finding the optimal block decomposition of the model ..." << endl;
-
- computeBlockDecomposition(prologue, epilogue);
-
- reduceBlockDecomposition();
-
- printBlockDecomposition();
-
- computeChainRuleJacobian();
+ /* Must be called after computeTemporaryTerms(), because it depends on
+ temporary_terms_mlv to be filled */
+ if (paramsDerivsOrder > 0 && !no_tmp_terms)
+ computeParamsDerivativesTemporaryTerms();
- determineLinearBlocks();
-
- computeBlockDynJacobianCols();
-
- if (!no_tmp_terms)
- computeBlockTemporaryTerms();
- }
- else
+ if (!computingPassBlock(eval_context, no_tmp_terms) && block)
{
- computeTemporaryTerms(!use_dll, no_tmp_terms);
-
- /* Must be called after computeTemporaryTerms(), because it depends on
- temporary_terms_mlv to be filled */
- if (paramsDerivsOrder > 0 && !no_tmp_terms)
- computeParamsDerivativesTemporaryTerms();
+ 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, true))
+ 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 dynamic model..." << endl;
+ computeBlockDecomposition(prologue, epilogue);
+ reduceBlockDecomposition();
+ printBlockDecomposition();
+ computeChainRuleJacobian();
+ determineLinearBlocks();
+ computeBlockDynJacobianCols();
+ if (!no_tmp_terms)
+ computeBlockTemporaryTerms();
+ return true;
+}
+
void
DynamicModel::computeXrefs()
{
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 8232d22a5a69253f51c05ab77a19a2cf00c50f02..37051525e6b6b90bcd2854b076476a8b7070ddc8 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -297,6 +297,12 @@ 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);
+
public:
DynamicModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
@@ -313,7 +319,7 @@ public:
//! Write cross references
void writeXrefs(ostream &output) const;
- //! Execute computations (variable sorting + derivation)
+ //! Execute computations (variable sorting + derivation + block decomposition)
/*!
\param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
\param derivsOrder order of derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true when order >= 2)
diff --git a/src/ModFile.cc b/src/ModFile.cc
index a0ddb4d8d75321bafa42d4a3743a2180e57c571d..bb0b8ca8c82de104a167dd93f762c0982fea3dee 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -496,7 +496,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
else
pos = pos2;
assert(pos);
- const StaticModel &planner_objective = pos->getPlannerObjective();
+ const PlannerObjective &planner_objective = pos->getPlannerObjective();
/*
clone the model then clone the new equations back to the original because
diff --git a/src/ModelEquationBlock.cc b/src/ModelEquationBlock.cc
index 9b0551f46b7677308ca8f0789f6fbc3da28230c7..fe0d649146333b3842a4c4d7ef06b82ada493b87 100644
--- a/src/ModelEquationBlock.cc
+++ b/src/ModelEquationBlock.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2010-2021 Dynare Team
+ * Copyright © 2010-2022 Dynare Team
*
* This file is part of Dynare.
*
@@ -23,6 +23,21 @@
#include "ModelEquationBlock.hh"
+PlannerObjective::PlannerObjective(SymbolTable &symbol_table_arg,
+ NumericalConstants &num_constants_arg,
+ ExternalFunctionsTable &external_functions_table_arg) :
+ StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg}
+{
+}
+
+bool
+PlannerObjective::computingPassBlock([[maybe_unused]] const eval_context_t &eval_context,
+ [[maybe_unused]] bool no_tmp_terms)
+{
+ // Disable block decomposition on planner objective
+ return false;
+}
+
SteadyStateModel::SteadyStateModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
ExternalFunctionsTable &external_functions_table_arg,
@@ -488,3 +503,11 @@ Epilogue::writeOutput(ostream &output) const
symbol_list.push_back(symbol_table.getName(symb_id));
SymbolList{move(symbol_list)}.writeOutput("M_.epilogue_var_list_", output);
}
+
+bool
+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 d2382cdb715ae985007042c7146149019eb8b47c..eb255e03ef61e6e53152fcb2fb9babb5e4cea6db 100644
--- a/src/ModelEquationBlock.hh
+++ b/src/ModelEquationBlock.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2010-2019 Dynare Team
+ * Copyright © 2010-2022 Dynare Team
*
* This file is part of Dynare.
*
@@ -26,6 +26,16 @@
#include "DynamicModel.hh"
#include "WarningConsolidation.hh"
+class PlannerObjective : public StaticModel
+{
+public:
+ PlannerObjective(SymbolTable &symbol_table_arg,
+ NumericalConstants &num_constants_arg,
+ ExternalFunctionsTable &external_functions_table_arg);
+private:
+ bool computingPassBlock(const eval_context_t &eval_context, bool no_tmp_terms) override;
+};
+
class SteadyStateModel : public DataTree
{
private:
@@ -98,6 +108,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;
};
#endif
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index d69b605ba525996dc88b9dbb10968908a4c127ef..8da7b1ad8a2ebd67ca8e79e9b3dc299aadc2acf8 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -242,7 +242,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
it != mate_map.begin() + n)
{
if (verbose)
- cerr << "ERROR: Could not normalize the " << (dynamic ? "dynamic" : "static") << " model. Variable "
+ cerr << "Could not normalize the " << (dynamic ? "dynamic" : "static") << " model. Variable "
<< symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
<< " is not in the maximum cardinality matching." << endl;
check = false;
@@ -250,13 +250,21 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
return check;
}
-void
+bool
ModelTree::computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, bool dynamic)
{
- cout << "Normalizing the " << (dynamic ? "dynamic" : "static") << " model..." << endl;
-
int n = equations.size();
+ /* Optimal policy models (discretionary, or Ramsey before computing FOCs) do
+ not have as many equations as variables. */
+ if (n != symbol_table.endo_nbr())
+ {
+ cout << "The " << (dynamic ? "dynamic" : "static") << " model cannot be normalized, since it does not have as many equations as variables." << endl;
+ return false;
+ }
+
+ cout << "Normalizing the " << (dynamic ? "dynamic" : "static") << " model..." << endl;
+
// Compute the maximum value of each row of the contemporaneous Jacobian matrix
vector max_val(n, 0.0);
for (const auto &[eq_and_endo, val] : contemporaneous_jacobian)
@@ -306,18 +314,9 @@ ModelTree::computeNonSingularNormalization(const jacob_map_t &contemporaneous_ja
found_normalization = computeNormalization(symbolic_jacobian, dynamic, true);
}
- if (!found_normalization)
- {
- /* Some models don’t have a steady state, and this can cause the
- normalization to fail (e.g. if some variable only appears in a diff(),
- it will disappear from the static model). Suggest the “no_static”
- option as a possible solution. */
- if (!dynamic)
- cerr << "If your model does not have a steady state, you may want to try the 'no_static' option of the 'model' block." << endl;
- /* The last call to computeNormalization(), which was verbose, already
- printed an error message, so we can immediately exit. */
- exit(EXIT_FAILURE);
- }
+ /* NB: If normalization failed, an explanatory message has been printed by the last call
+ to computeNormalization(), which has verbose=true */
+ return found_normalization;
}
ModelTree::jacob_map_t
@@ -332,25 +331,20 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context) const
int eq = indices[0];
int var { getTypeSpecificIDByDerivID(deriv_id) };
int lag = getLagByDerivID(deriv_id);
- double val = 0;
- try
- {
- val = d1->eval(eval_context);
- }
- catch (ExprNode::EvalExternalFunctionException &e)
- {
- val = 1;
- }
- catch (ExprNode::EvalException &e)
- {
- cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1;
- if (equations_lineno[eq])
- cerr << " (line " << *equations_lineno[eq] << ")";
- cerr << " and variable " << getNameByDerivID(deriv_id) << "(" << lag << ") !" << endl;
- d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModel, {}, {});
- cerr << endl;
- exit(EXIT_FAILURE);
- }
+ double val { [&]
+ {
+ try
+ {
+ return d1->eval(eval_context);
+ }
+ catch (ExprNode::EvalExternalFunctionException &e)
+ {
+ return 1.0;
+ }
+ /* Other types of EvalException should not happen (all symbols should
+ have a value; we don’t evaluate an equal sign) */
+ }() };
+
if ((isnan(val) || fabs(val) >= cutoff) && lag == 0)
contemporaneous_jacobian[{ eq, var }] = val;
}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index a03e96597f35b7d0dbdb7ad8b789d5591d423a6f..ecb49614c135e0dcb464a81039c0e50759b443b7 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -352,8 +352,9 @@ protected:
If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e. use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to reflect all the edges in the symbolic incidence matrix.
If no matching is found with a zero cutoff, an error message is printed.
The resulting normalization is stored in endo2eq.
+ Returns a boolean indicating success.
*/
- void computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, bool dynamic);
+ bool computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, bool dynamic);
//! Evaluate the jacobian (w.r.t. endogenous) and suppress all the elements below the cutoff
/*! Returns the contemporaneous_jacobian.
Elements below the cutoff are discarded. External functions are evaluated to 1. */
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 9d6e27da116a48c23772c1126d410cc007f6786a..39110bfb6f9cc1f51576e8d03866318c2145526e 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2154,9 +2154,9 @@ ParsingDriver::run_model_comparison()
void
ParsingDriver::begin_planner_objective()
{
- planner_objective = make_unique<StaticModel>(mod_file->symbol_table,
- mod_file->num_constants,
- mod_file->external_functions_table);
+ planner_objective = make_unique<PlannerObjective>(mod_file->symbol_table,
+ mod_file->num_constants,
+ mod_file->external_functions_table);
set_current_data_tree(planner_objective.get());
}
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index 6cafa8daf321bc46f81524b00e570fb606ee8238..b91aeabf5c58d8f77264ba428b0fec022f8d1bd3 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -111,7 +111,7 @@ private:
int declare_symbol(const string &name, SymbolType type, const string &tex_name, const vector<pair<string, string>> &partition_value);
//! Temporary store for the planner objective
- unique_ptr<StaticModel> planner_objective;
+ unique_ptr<PlannerObjective> planner_objective;
//! Temporary store for the occbin_constraints statement
unique_ptr<DataTree> occbin_constraints_tree;
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 45d8723cb0eaf73a5ecbd2f6c59b3c78443a16a6..8636bd76407a9f908c53d98ddd44c6c84b205a85 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -23,6 +23,7 @@
#include <cassert>
#include <algorithm>
#include <sstream>
+#include <numeric>
#include "StaticModel.hh"
#include "DynamicModel.hh"
@@ -245,6 +246,10 @@ StaticModel::writeStaticBytecode(const string &basename) const
int u_count_int { writeBytecodeBinFile(basename + "/model/bytecode/static.bin", false) };
BytecodeWriter code_file {basename + "/model/bytecode/static.cod"};
+ vector<int> eq_idx(equations.size());
+ iota(eq_idx.begin(), eq_idx.end(), 0);
+ vector<int> endo_idx(symbol_table.endo_nbr());
+ iota(endo_idx.begin(), endo_idx.end(), 0);
// Declare temporary terms and the (single) block
code_file << FDIMST_{static_cast<int>(temporary_terms_derivatives[0].size()
@@ -253,8 +258,8 @@ StaticModel::writeStaticBytecode(const string &basename) const
BlockSimulationType::solveForwardComplete,
0,
symbol_table.endo_nbr(),
- endo_idx_block2orig,
- eq_idx_block2orig,
+ endo_idx,
+ eq_idx,
false,
symbol_table.endo_nbr(),
0,
@@ -362,44 +367,40 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
computeParamsDerivatives(paramsDerivsOrder);
}
- if (block)
- {
- auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
-
- computeNonSingularNormalization(contemporaneous_jacobian, 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 model ..." << endl;
-
- computeBlockDecomposition(prologue, epilogue);
-
- reduceBlockDecomposition();
-
- printBlockDecomposition();
+ computeTemporaryTerms(true, no_tmp_terms);
- computeChainRuleJacobian();
+ /* Must be called after computeTemporaryTerms(), because it depends on
+ temporary_terms_mlv to be filled */
+ if (paramsDerivsOrder > 0 && !no_tmp_terms)
+ computeParamsDerivativesTemporaryTerms();
- determineLinearBlocks();
-
- if (!no_tmp_terms)
- computeBlockTemporaryTerms();
- }
- else
+ if (!computingPassBlock(eval_context, no_tmp_terms) && block)
{
- computeTemporaryTerms(true, no_tmp_terms);
-
- /* Must be called after computeTemporaryTerms(), because it depends on
- temporary_terms_mlv to be filled */
- if (paramsDerivsOrder > 0 && !no_tmp_terms)
- computeParamsDerivativesTemporaryTerms();
+ 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, false))
+ 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 static model..." << 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 092ebdc046af0fb1ef5feebe99504963d6f088e6..da65283905cc5a4069e114d64bbf491744138a5f 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -110,6 +110,12 @@ 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);
+
public:
StaticModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants,
@@ -124,7 +130,7 @@ public:
//! Writes information about the static model to the driver file
void writeDriverOutput(ostream &output, bool block) const;
- //! Execute computations (variable sorting + derivation)
+ //! Execute computations (variable sorting + derivation + block decomposition)
/*!
\param eval_context evaluation context for normalization
\param no_tmp_terms if true, no temporary terms will be computed in the static files