From c48248fc0d65ed70fb13dd508ed0254f113cbcb6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 30 Nov 2022 14:43:44 +0100
Subject: [PATCH] Implement time-recursive block decomposition tuned for purely
backward/forward/static models
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Such a decomposition has to be simulated with periods as the outer loop and
blocks as the inner loop.
It is enabled by default for purely backward/forward/static models, as long as
the “block” option is not given. In that case, “mfs” is also set to 3 by
default (until that value becomes the global default).
M_.time_recursive_block_decomposition is set to “true” when that decomposition
has been performed, “false” otherwise for the traditional decomposition (the
latter has to be simulated with blocks as the outer loop and periods as the
inner loop).
---
src/DynamicModel.cc | 11 ++++++++++-
src/ModelTree.cc | 16 ++++++++++------
src/ModelTree.hh | 19 +++++++++++++++++++
3 files changed, 39 insertions(+), 7 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 1a2def69..0c50e7db 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1350,6 +1350,8 @@ void
DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
const vector<int> &state_var, bool estimation_present) const
{
+ output << "M_.block_structure.time_recursive = " << boolalpha << time_recursive_block_decomposition << ";" << endl;
+
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
{
int block_size = blocks[blk].size;
@@ -2959,6 +2961,11 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_c
if (paramsDerivsOrder > 0 && !no_tmp_terms)
computeParamsDerivativesTemporaryTerms();
+ if (!block && (max_endo_lag == 0 || max_endo_lead == 0))
+ {
+ time_recursive_block_decomposition = true;
+ mfs = 3; // FIXME: remove this line when mfs=3 becomes the global default
+ }
computingPassBlock(eval_context, no_tmp_terms);
if (block_decomposed)
computeBlockDynJacobianCols();
@@ -3077,7 +3084,9 @@ DynamicModel::determineBlockDerivativesType(int blk)
map<tuple<int, int, int>, BlockDerivativeType> derivType;
int size = blocks[blk].size;
int nb_recursive = blocks[blk].getRecursiveSize();
- for (int lag = -blocks[blk].max_lag; lag <= blocks[blk].max_lead; lag++)
+ for (int lag {time_recursive_block_decomposition ? 0 : -blocks[blk].max_lag};
+ lag <= (time_recursive_block_decomposition ? 0 : blocks[blk].max_lead);
+ lag++)
for (int eq = 0; eq < size; eq++)
{
set<pair<int, int>> endos_and_lags;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 372ca706..51060f45 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -166,6 +166,7 @@ ModelTree::ModelTree(const ModelTree &m) :
eq_idx_orig2block{m.eq_idx_orig2block},
endo_idx_orig2block{m.endo_idx_orig2block},
block_decomposed{m.block_decomposed},
+ time_recursive_block_decomposition{m.time_recursive_block_decomposition},
blocks{m.blocks},
endo2block{m.endo2block},
eq2block{m.eq2block},
@@ -210,6 +211,7 @@ ModelTree::operator=(const ModelTree &m)
equation_type_and_normalized_equation.clear();
blocks_derivatives.clear();
block_decomposed = m.block_decomposed;
+ time_recursive_block_decomposition = m.time_recursive_block_decomposition;
blocks = m.blocks;
endo2block = m.endo2block;
eq2block = m.eq2block;
@@ -408,7 +410,8 @@ ModelTree::computePrologueAndEpilogue()
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;
+ if (!time_recursive_block_decomposition || lag == 0)
+ IM[i * n + endo2eq[endo]] = true;
}
bool something_has_been_done;
@@ -662,7 +665,7 @@ ModelTree::computeBlockDecomposition(int prologue, int epilogue)
For detecting dependencies between variables, use the symbolic adjacency
matrix */
VariableDependencyGraph G(nb_simvars);
- for (const auto &[key, value] : computeSymbolicJacobian(false))
+ for (const auto &[key, value] : computeSymbolicJacobian(time_recursive_block_decomposition))
{
auto [eq, endo] = key;
if (eq_idx_orig2block[eq] >= prologue
@@ -722,10 +725,11 @@ ModelTree::computeBlockDecomposition(int prologue, int epilogue)
feedback set */
for (int i = 0; i < nb_simvars; i++)
if (equation_type_and_normalized_equation[eq_idx_block2orig[i+prologue]].first == EquationType::solve
- || variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0
- || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
- || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
- || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0
+ || (!time_recursive_block_decomposition &&
+ (variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0
+ || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
+ || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
+ || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0))
|| mfs == 0)
add_edge(vertex(i, G), vertex(i, G), G);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index c6312782..de892f53 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -212,6 +212,25 @@ protected:
// Whether block decomposition has been successfully computed
bool block_decomposed {false};
+ /* Whether the block decomposition to compute is time-recursive (i.e. the
+ model can be simulated as a whole period-by-period).
+
+ If true, only contemporaneous occurrences of variables are considered when
+ computing the block structure; leads and lags are essentially treated as
+ exogenous, i.e. they are ignored. Such a decomposition only makes sense
+ for models that are purely backward/forward/static. When using the
+ resulting block decomposition to simulate the model, periods must be the
+ outer loop, and blocks the inner loop.
+
+ If false, then the full lead/lag structure is taken into account when
+ computing the block structure. This is the only option if there are both
+ leads and lags. When using the
+ resulting block decomposition to simulate the model, blocks must be the
+ outer loop, and periods the inner loop.
+
+ Of course, this setting does not make any difference on StaticModel. */
+ bool time_recursive_block_decomposition {false};
+
// Stores various informations on the blocks
vector<BlockInfo> blocks;
--
GitLab