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

Factorize computing pass for block decomposition 

Also add “block_decomposed” data member for tracking whether the block
decomposition has been successful.
master
Sébastien Villemot 2022-09-28 16:31:51 +02:00
parent 6e69eea3cf
commit 4c58451d83
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
8 changed files with 97 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
src/DynamicModel.cc
View File

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

9
src/DynamicModel.hh
View File

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

6
src/ModelEquationBlock.cc
View File

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

4
src/ModelEquationBlock.hh
View File

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

20
src/ModelTree.cc
View File

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

17
src/ModelTree.hh
View File

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

23
src/StaticModel.cc
View File

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

9
src/StaticModel.hh
View File

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