From 2d3e3eff6f938fb58a893ee1153e45e61bfd26d4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Mon, 24 Apr 2023 17:49:54 +0200
Subject: [PATCH] Block decomposition: add specialized normalization algorithm
 for purely backward models

If the model is purely backward, determine whether all original equations have
a single contemporaneous endogenous on the LHS. If this is the case, then first
try a normalization by enforcing that each original equation is matched with
the endogenous on the LHS.

This helps with the simulation of purely backward models, where equations are
renormalized with mfs=3, since it produces a simpler system to be recursively
evaluated/solved.
---
 src/ModelTree.cc | 48 ++++++++++++++++++++++++++++++++++++++++++++++++
 src/ModelTree.hh | 11 +++++++++++
 2 files changed, 59 insertions(+)

diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 7eef0981..72320d52 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -283,6 +283,25 @@ ModelTree::computeNonSingularNormalization(const eval_context_t &eval_context)
 
   cout << "Normalizing the " << modelClassName() << "..." << endl;
 
+  /* If the model is purely backward, determine whether all original equations
+     have a single contemporaneous endogenous on the LHS. If this is the case,
+     then first try a normalization by enforcing that each original equation is
+     matched with the endogenous on the LHS. */
+  if (time_recursive_block_decomposition)
+    {
+      auto [normalize_by_lhs, lhs_symbolic_jacobian] { computeLeftHandSideSymbolicJacobian() };
+      if (normalize_by_lhs)
+        try
+          {
+            computeNormalization(lhs_symbolic_jacobian);
+            return true;
+          }
+        catch (ModelNormalizationFailed &e)
+          {
+            cerr << "WARNING: All equations are written so that a single contemporaneous endogenous variable appears on the left-hand side. This suggests a natural normalization of the model. However, variable " << e.unmatched_endo << " could not be matched with an equation. Check whether this is desired." << endl;
+          }
+    }
+
   auto contemporaneous_jacobian {evaluateAndReduceJacobian(eval_context)};
 
   // Compute the maximum value of each row of the contemporaneous Jacobian matrix
@@ -1826,6 +1845,35 @@ ModelTree::computeSymbolicJacobian(bool contemporaneous_only) const
   return symbolic_jacobian;
 }
 
+pair<bool, ModelTree::jacob_map_t>
+ModelTree::computeLeftHandSideSymbolicJacobian() const
+{
+  jacob_map_t lhs_symbolic_jacobian;
+  auto not_contemporaneous = [](const pair<int, int> &p) { return p.second != 0; };
+
+  for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
+    if (equations_lineno[eq]) // Hand-written equation: test whether LHS has single contemporaneous endo
+      {
+        set<pair<int, int>> endos_and_lags;
+        equations[eq]->arg1->collectEndogenous(endos_and_lags);
+        erase_if(endos_and_lags, not_contemporaneous);
+        if (endos_and_lags.size() == 1)
+          lhs_symbolic_jacobian.try_emplace({ eq, endos_and_lags.begin()->first }, 1);
+        else
+          return { false, {} };
+      }
+    else // Generated equation: keep all endos on both LHS and RHS
+      {
+        set<pair<int, int>> endos_and_lags;
+        equations[eq]->collectEndogenous(endos_and_lags);
+        erase_if(endos_and_lags, not_contemporaneous);
+        for (const auto &[endo, lag] : endos_and_lags)
+          lhs_symbolic_jacobian.try_emplace({ eq, endo }, 1);
+      }
+
+  return { true, lhs_symbolic_jacobian };
+}
+
 void
 ModelTree::updateReverseVariableEquationOrderings()
 {
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 7ba465a5..2cc2be0c 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -433,6 +433,17 @@ private:
      variables; otherwise also considers leads and lags. */
   jacob_map_t computeSymbolicJacobian(bool contemporaneous_only) const;
 
+  /* Compute a pseudo-Jacobian whose all elements are either zero or one.
+     For the equations that were originally written by the user (identified as
+     those having an associated line number), checks whether there is a single
+     contemporaneous endogenous on the left-hand side; if yes, only this
+     endogenous is associated with a one on the line of the corresponding
+     equation; otherwise, returns false as the first output argument and
+     aborts the computation.
+     For the other equations, fills the corresponding lines as is done
+     by computeSymbolicJacobian(true). */
+  pair<bool, jacob_map_t> computeLeftHandSideSymbolicJacobian() const;
+
   // Compute {var,eq}_idx_orig2block from {var,eq}_idx_block2orig
   void updateReverseVariableEquationOrderings();