From 53add0b2fe94b8358fd48ecaf939f9ca0de00c3f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Mon, 28 Nov 2022 13:40:49 +0100
Subject: [PATCH] Block decomposition: when falling back to symbolic
 normalization, use contemporaneous symbolic Jacobian

Previously, it would use a symbolic Jacobian where leads and lags were also
taken into account.
---
 src/ModelTree.cc | 12 ++++++------
 src/ModelTree.hh |  6 ++++--
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index ef94feaa..b8fb83fd 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -331,9 +331,8 @@ ModelTree::computeNonSingularNormalization(const jacob_map_t &contemporaneous_ja
     {
       cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
       /* If no non-singular normalization can be found, try to find a
-         normalization even with a potential singularity.
-         TODO: Explain why symbolic_jacobian is not contemporaneous. */
-      auto symbolic_jacobian = computeSymbolicJacobian();
+         normalization even with a potential singularity. */
+      auto symbolic_jacobian = computeSymbolicJacobian(true);
       found_normalization = computeNormalization(symbolic_jacobian, true);
     }
 
@@ -661,7 +660,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())
+  for (const auto &[key, value] : computeSymbolicJacobian(false))
     {
       auto [eq, endo] = key;
       if (eq_idx_orig2block[eq] >= prologue
@@ -1878,7 +1877,7 @@ ModelTree::collectFirstOrderDerivativesEndogenous()
 }
 
 ModelTree::jacob_map_t
-ModelTree::computeSymbolicJacobian() const
+ModelTree::computeSymbolicJacobian(bool contemporaneous_only) const
 {
   jacob_map_t symbolic_jacobian;
   for (int i = 0; i < static_cast<int>(equations.size()); i++)
@@ -1886,7 +1885,8 @@ ModelTree::computeSymbolicJacobian() const
       set<pair<int, int>> endos_and_lags;
       equations[i]->collectEndogenous(endos_and_lags);
       for (const auto &[endo, lag] : endos_and_lags)
-        symbolic_jacobian[{ i, endo }] = 1;
+        if (!contemporaneous_only || lag == 0)
+          symbolic_jacobian[{ i, endo }] = 1;
     }
   return symbolic_jacobian;
 }
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 8b080513..c6312782 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -424,8 +424,10 @@ private:
   static set<filesystem::path> mex_compilation_done;
 
   /* Compute a pseudo-Jacobian whose all elements are either zero or one,
-     depending on whether the variable symbolically appears in the equation */
-  jacob_map_t computeSymbolicJacobian() const;
+     depending on whether the variable symbolically appears in the equation. If
+     contemporaneous_only=true, only considers contemporaneous occurences of
+     variables; otherwise also considers leads and lags. */
+  jacob_map_t computeSymbolicJacobian(bool contemporaneous_only) const;
 
   // Compute {var,eq}_idx_orig2block from {var,eq}_idx_block2orig
   void updateReverseVariableEquationOrderings();
-- 
GitLab