stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Always compute block decomposition, even if “block” option was not passed 

If block decomposition fails, error out if “block” option was passed, but not
otherwise.

This commit does not modify the generated files.

This is a preliminary step for dynare#1859.
master
Sébastien Villemot 2022-09-13 15:43:13 +02:00
parent 1aead92cd5
commit d67f569035
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
13 changed files with 162 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/ComputingTasks.cc
View File

@ -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;

6
src/ComputingTasks.hh
View File

@ -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

72
src/DynamicModel.cc
View File

@ -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,45 +3289,39 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeParamsDerivatives(paramsDerivsOrder);
}
computeTemporaryTerms(!use_dll, no_tmp_terms);
if (block)
/* Must be called after computeTemporaryTerms(), because it depends on
temporary_terms_mlv to be filled */
if (paramsDerivsOrder > 0 && !no_tmp_terms)
computeParamsDerivativesTemporaryTerms();
if (!computingPassBlock(eval_context, no_tmp_terms) && 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();
determineLinearBlocks();
computeBlockDynJacobianCols();
if (!no_tmp_terms)
computeBlockTemporaryTerms();
cerr << "ERROR: Block decomposition requested but failed." << endl;
exit(EXIT_FAILURE);
}
else
{
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();
}
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

8
src/DynamicModel.hh
View File

@ -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)

2
src/ModFile.cc
View File

@ -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

25
src/ModelEquationBlock.cc
View File

@ -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;
}

13
src/ModelEquationBlock.hh
View File

@ -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

64
src/ModelTree.cc
View File

@ -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;
}

3
src/ModelTree.hh
View File

@ -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. */

6
src/ParsingDriver.cc
View File

@ -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());
}

2
src/ParsingDriver.hh
View File

@ -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;

71
src/StaticModel.cc
View File

@ -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,42 +367,38 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
computeParamsDerivatives(paramsDerivsOrder);
}
if (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();
if (!computingPassBlock(eval_context, no_tmp_terms) && 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();
computeChainRuleJacobian();
determineLinearBlocks();
if (!no_tmp_terms)
computeBlockTemporaryTerms();
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);
}
else
{
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();
}
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

8
src/StaticModel.hh
View File

@ -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