Implement time-recursive block decomposition tuned for purely backward/forward/static models

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

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

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

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;