From 325e6487d4307e452529ea8bc374af07eecd879e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 24 Feb 2023 15:12:07 -0500
Subject: [PATCH] Bytecode: move computation routines out of the Evaluate class

---
 mex/sources/bytecode/Evaluate.cc     | 217 ---------------------------
 mex/sources/bytecode/Evaluate.hh     |  16 +-
 mex/sources/bytecode/Interpreter.cc  | 135 +++++++++++++++++
 mex/sources/bytecode/Interpreter.hh  |   4 +
 mex/sources/bytecode/SparseMatrix.cc |  82 ++++++++++
 mex/sources/bytecode/SparseMatrix.hh |  11 ++
 6 files changed, 233 insertions(+), 232 deletions(-)

diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index 222fd4cbd..c474347f2 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/mex/sources/bytecode/Evaluate.cc
@@ -2255,107 +2255,6 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
 #endif
 }
 
-void
-Evaluate::evaluate_over_periods(bool forward)
-{
-  if (steady_state)
-    compute_block_time(0, false, false);
-  else
-    {
-      if (forward)
-        {
-          for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
-            compute_block_time(0, false, false);
-          it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
-        }
-      else
-        {
-          for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
-            compute_block_time(0, false, false);
-          it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
-        }
-    }
-}
-
-void
-Evaluate::solve_simple_one_periods()
-{
-  bool cvg = false;
-  int iter = 0;
-  double ya;
-  double slowc = 1;
-  res1 = 0;
-  while (!(cvg || iter > maxit_))
-    {
-      Per_y_ = it_*y_size;
-      ya = y[Block_Contain[0].Variable + Per_y_];
-      compute_block_time(0, false, false);
-      if (!isfinite(res1))
-        {
-          res1 = std::numeric_limits<double>::quiet_NaN();
-          while ((isinf(res1) || isnan(res1)) && (slowc > 1e-9))
-            {
-              compute_block_time(0, false, false);
-              if (!isfinite(res1))
-                {
-                  slowc /= 1.5;
-                  mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
-                  y[Block_Contain[0].Variable + Per_y_] = ya - slowc * divide(r[0], g1[0]);
-                }
-            }
-        }
-      double rr;
-      rr = r[0];
-      cvg = (fabs(rr) < solve_tolf);
-      if (cvg)
-        continue;
-      try
-        {
-          y[Block_Contain[0].Variable + Per_y_] += -slowc *divide(rr, g1[0]);
-        }
-      catch (FloatingPointException &fpeh)
-        {
-          mexPrintf("%s\n      \n", fpeh.message.c_str());
-          mexPrintf("      Singularity in block %d", block_num+1);
-        }
-      iter++;
-    }
-  if (!cvg)
-    throw FatalException{"In Solve Forward simple, convergence not achieved in block "
-                         + to_string(block_num+1) + ", after " + to_string(iter) + " iterations"};
-}
-
-void
-Evaluate::solve_simple_over_periods(bool forward)
-{
-  g1 = static_cast<double *>(mxMalloc(sizeof(double)));
-  test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double));
-  r = static_cast<double *>(mxMalloc(sizeof(double)));
-  test_mxMalloc(r, __LINE__, __FILE__, __func__, sizeof(double));
-  if (steady_state)
-    {
-      it_ = 0;
-      solve_simple_one_periods();
-    }
-  else
-    {
-      if (forward)
-        {
-          for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
-            solve_simple_one_periods();
-          it_= periods+y_kmin-1; // Do not leave it_ in inconsistent state
-        }
-      else
-        {
-          for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
-            solve_simple_one_periods();
-          it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
-        }
-    }
-  mxFree(g1);
-  mxFree(r);
-}
-
 void
 Evaluate::gotoBlock(int block)
 {
@@ -2373,122 +2272,6 @@ Evaluate::gotoBlock(int block)
   u_count_int = fb->get_u_count_int();
 }
 
-void
-Evaluate::compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx)
-{
-  res1 = 0;
-  *_res1 = 0;
-  *_res2 = 0;
-  *_max_res = 0;
-  for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
-    {
-      Per_u_ = (it_-y_kmin)*u_count_int;
-      Per_y_ = it_*y_size;
-      int shift = (it_-y_kmin) * size;
-      compute_block_time(Per_u_, false, no_derivatives);
-      if (!(isnan(res1) || isinf(res1)))
-        for (int i = 0; i < size; i++)
-          {
-            double rr;
-            rr = r[i];
-            res[i+shift] = rr;
-            if (max_res < fabs(rr))
-              {
-                *_max_res = fabs(rr);
-                *_max_res_idx = i;
-              }
-            *_res2 += rr*rr;
-            *_res1 += fabs(rr);
-          }
-      else
-        return;
-    }
-  it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
-  return;
-}
-
-bool
-Evaluate::compute_complete(bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
-{
-  bool result;
-  res1 = 0;
-  compute_block_time(0, false, no_derivatives);
-  if (!(isnan(res1) || isinf(res1)))
-    {
-      _res1 = 0;
-      _res2 = 0;
-      _max_res = 0;
-      for (int i = 0; i < size; i++)
-        {
-          double rr;
-          rr = r[i];
-          if (max_res < fabs(rr))
-            {
-              _max_res = fabs(rr);
-              _max_res_idx = i;
-            }
-          _res2 += rr*rr;
-          _res1 += fabs(rr);
-        }
-      result = true;
-    }
-  else
-    result = false;
-  return result;
-}
-
-bool
-Evaluate::compute_complete(double lambda, double *crit)
-{
-  double res1_ = 0, res2_ = 0, max_res_ = 0;
-  int max_res_idx_ = 0;
-  if (steady_state)
-    {
-      it_ = 0;
-      for (int i = 0; i < size; i++)
-        {
-          int eq = index_vara[i];
-          y[eq] = ya[eq] + lambda * direction[eq];
-        }
-      Per_u_ = 0;
-      Per_y_ = 0;
-      if (compute_complete(true, res1, res2, max_res, max_res_idx))
-        res2_ = res2;
-      else
-        return false;
-    }
-  else
-    {
-      for (int it = y_kmin; it < periods+y_kmin; it++)
-        for (int i = 0; i < size; i++)
-          {
-            int eq = index_vara[i];
-            y[eq+it*y_size] = ya[eq+it*y_size] + lambda * direction[eq+it*y_size];
-          }
-      for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
-        {
-          Per_u_ = (it_-y_kmin)*u_count_int;
-          Per_y_ = it_*y_size;
-          if (compute_complete(true, res1, res2, max_res, max_res_idx))
-            {
-              res2_ += res2;
-              res1_ += res1;
-              if (max_res > max_res_)
-                {
-                  max_res = max_res_;
-                  max_res_idx = max_res_idx_;
-                }
-            }
-          else
-            return false;
-        }
-      it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
-    }
-  mexPrintf("  lambda=%e, res2=%e\n", lambda, res2_);
-  *crit = res2_/2;
-  return true;
-}
-
 void
 Evaluate::printCurrentBlock()
 {
diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh
index 086e9a495..16b22aa32 100644
--- a/mex/sources/bytecode/Evaluate.hh
+++ b/mex/sources/bytecode/Evaluate.hh
@@ -93,20 +93,10 @@ protected:
   double divide(double a, double b);
   double log1(double a);
   double log10_1(double a);
-  void evaluate_over_periods(bool forward);
-  void solve_simple_one_periods();
-  void solve_simple_over_periods(bool forward);
   void compute_block_time(int Per_u_, bool evaluate, bool no_derivatives);
-  int Per_u_, Per_y_;
   int it_;
-  int maxit_;
-  double *direction;
-  double solve_tolf;
   bool print_error;
-  double res1, res2, max_res;
-  int max_res_idx;
-
-  int *index_vara;
+  double res1;
 
   int block_num; // Index of the current block
   int size; // Size of the current block
@@ -161,10 +151,6 @@ public:
   Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool steady_state_arg, int periods_arg, BasicSymbolTable &symbol_table_arg);
   // TODO: integrate into the constructor
   void loadCodeFile(const filesystem::path &codfile);
-  bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
-  void compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
-
-  bool compute_complete(double lambda, double *crit);
 
   int
   get_block_number() const
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 9515de231..8f4e06971 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -65,6 +65,141 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
   print_it = print_it_arg;
 }
 
+void
+Interpreter::evaluate_over_periods(bool forward)
+{
+  if (steady_state)
+    compute_block_time(0, false, false);
+  else
+    {
+      if (forward)
+        {
+          for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
+            compute_block_time(0, false, false);
+          it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
+        }
+      else
+        {
+          for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
+            compute_block_time(0, false, false);
+          it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
+        }
+    }
+}
+
+void
+Interpreter::solve_simple_one_periods()
+{
+  bool cvg = false;
+  int iter = 0;
+  double ya;
+  double slowc = 1;
+  res1 = 0;
+  while (!(cvg || iter > maxit_))
+    {
+      Per_y_ = it_*y_size;
+      ya = y[Block_Contain[0].Variable + Per_y_];
+      compute_block_time(0, false, false);
+      if (!isfinite(res1))
+        {
+          res1 = std::numeric_limits<double>::quiet_NaN();
+          while ((isinf(res1) || isnan(res1)) && (slowc > 1e-9))
+            {
+              compute_block_time(0, false, false);
+              if (!isfinite(res1))
+                {
+                  slowc /= 1.5;
+                  mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
+                  y[Block_Contain[0].Variable + Per_y_] = ya - slowc * divide(r[0], g1[0]);
+                }
+            }
+        }
+      double rr;
+      rr = r[0];
+      cvg = (fabs(rr) < solve_tolf);
+      if (cvg)
+        continue;
+      try
+        {
+          y[Block_Contain[0].Variable + Per_y_] += -slowc *divide(rr, g1[0]);
+        }
+      catch (FloatingPointException &fpeh)
+        {
+          mexPrintf("%s\n      \n", fpeh.message.c_str());
+          mexPrintf("      Singularity in block %d", block_num+1);
+        }
+      iter++;
+    }
+  if (!cvg)
+    throw FatalException{"In Solve Forward simple, convergence not achieved in block "
+                         + to_string(block_num+1) + ", after " + to_string(iter) + " iterations"};
+}
+
+void
+Interpreter::solve_simple_over_periods(bool forward)
+{
+  g1 = static_cast<double *>(mxMalloc(sizeof(double)));
+  test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double));
+  r = static_cast<double *>(mxMalloc(sizeof(double)));
+  test_mxMalloc(r, __LINE__, __FILE__, __func__, sizeof(double));
+  if (steady_state)
+    {
+      it_ = 0;
+      solve_simple_one_periods();
+    }
+  else
+    {
+      if (forward)
+        {
+          for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
+            solve_simple_one_periods();
+          it_= periods+y_kmin-1; // Do not leave it_ in inconsistent state
+        }
+      else
+        {
+          for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
+            solve_simple_one_periods();
+          it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
+        }
+    }
+  mxFree(g1);
+  mxFree(r);
+}
+
+void
+Interpreter::compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx)
+{
+  res1 = 0;
+  *_res1 = 0;
+  *_res2 = 0;
+  *_max_res = 0;
+  for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
+    {
+      Per_u_ = (it_-y_kmin)*u_count_int;
+      Per_y_ = it_*y_size;
+      int shift = (it_-y_kmin) * size;
+      compute_block_time(Per_u_, false, no_derivatives);
+      if (!(isnan(res1) || isinf(res1)))
+        for (int i = 0; i < size; i++)
+          {
+            double rr;
+            rr = r[i];
+            res[i+shift] = rr;
+            if (max_res < fabs(rr))
+              {
+                *_max_res = fabs(rr);
+                *_max_res_idx = i;
+              }
+            *_res2 += rr*rr;
+            *_res1 += fabs(rr);
+          }
+      else
+        return;
+    }
+  it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
+  return;
+}
+
 void
 Interpreter::evaluate_a_block(bool initialization, bool single_block, const string &bin_base_name)
 {
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 99f41ab66..f38386813 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -41,6 +41,10 @@ private:
   bool global_temporary_terms;
   bool print;
   int col_x, col_y;
+  void evaluate_over_periods(bool forward);
+  void solve_simple_one_periods();
+  void solve_simple_over_periods(bool forward);
+  void compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
 protected:
   void evaluate_a_block(bool initialization, bool single_block, const string &bin_base_name);
   int simulate_a_block(const vector_table_conditional_local_type &vector_table_conditional_local, bool single_block, const string &bin_base_name);
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 12f04eafe..5cc4f27f5 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -1824,6 +1824,88 @@ dynSparseMatrix::Sparse_transpose(const mxArray *A_m)
   return C_m;
 }
 
+bool
+dynSparseMatrix::compute_complete(bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
+{
+  bool result;
+  res1 = 0;
+  compute_block_time(0, false, no_derivatives);
+  if (!(isnan(res1) || isinf(res1)))
+    {
+      _res1 = 0;
+      _res2 = 0;
+      _max_res = 0;
+      for (int i = 0; i < size; i++)
+        {
+          double rr;
+          rr = r[i];
+          if (max_res < fabs(rr))
+            {
+              _max_res = fabs(rr);
+              _max_res_idx = i;
+            }
+          _res2 += rr*rr;
+          _res1 += fabs(rr);
+        }
+      result = true;
+    }
+  else
+    result = false;
+  return result;
+}
+
+bool
+dynSparseMatrix::compute_complete(double lambda, double *crit)
+{
+  double res1_ = 0, res2_ = 0, max_res_ = 0;
+  int max_res_idx_ = 0;
+  if (steady_state)
+    {
+      it_ = 0;
+      for (int i = 0; i < size; i++)
+        {
+          int eq = index_vara[i];
+          y[eq] = ya[eq] + lambda * direction[eq];
+        }
+      Per_u_ = 0;
+      Per_y_ = 0;
+      if (compute_complete(true, res1, res2, max_res, max_res_idx))
+        res2_ = res2;
+      else
+        return false;
+    }
+  else
+    {
+      for (int it = y_kmin; it < periods+y_kmin; it++)
+        for (int i = 0; i < size; i++)
+          {
+            int eq = index_vara[i];
+            y[eq+it*y_size] = ya[eq+it*y_size] + lambda * direction[eq+it*y_size];
+          }
+      for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
+        {
+          Per_u_ = (it_-y_kmin)*u_count_int;
+          Per_y_ = it_*y_size;
+          if (compute_complete(true, res1, res2, max_res, max_res_idx))
+            {
+              res2_ += res2;
+              res1_ += res1;
+              if (max_res > max_res_)
+                {
+                  max_res = max_res_;
+                  max_res_idx = max_res_idx_;
+                }
+            }
+          else
+            return false;
+        }
+      it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
+    }
+  mexPrintf("  lambda=%e, res2=%e\n", lambda, res2_);
+  *crit = res2_/2;
+  return true;
+}
+
 bool
 dynSparseMatrix::mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc)
 {
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index 32a4733bd..e8eb2013f 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -195,6 +195,17 @@ protected:
   int stack_solve_algo, solve_algo;
   int minimal_solving_periods;
   bool print_it;
+  int Per_u_, Per_y_;
+  int maxit_;
+  double *direction;
+  double solve_tolf;
+  double res2, max_res;
+  int max_res_idx;
+  int *index_vara;
+
+  bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
+
+  bool compute_complete(double lambda, double *crit);
 };
 
 #endif // _SPARSEMATRIX_HH