From 85a0208e64a3c553987fb81748fe8d7d9bade384 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 21 Jul 2022 17:22:08 +0200
Subject: [PATCH] Bugfix with temporary terms in block+bytecode
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

There were actually two distinct bugs, leading to incorrect results in some
corner cases:

– in the “evaluate” mode of the bytecode MEX, the temporary terms of the
  derivatives “evaluate” blocks were not evaluated at runtime; but these
  temporary terms may be needed for residuals of subsequent blocks;

– when the bytecode MEX was only computing residuals of the model (and not 1st
  order derivatives), the temporary terms of the derivatives were not evaluated
  at runtime; but these temporary terms may be needed for residuals of subsequent
  blocks.

(manually cherry picked from commit a97a41f6c03c60c0b628d2a9caceec8e182b6222)
---
 src/DynamicModel.cc | 26 +++++++++++---------------
 src/StaticModel.cc  | 22 ++++++++++++++++------
 2 files changed, 27 insertions(+), 21 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 939fe3ed..e92589b0 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1214,17 +1214,20 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
               fstpr.write(code_file, instruction_number);
             }
         }
+
+      /* Write temporary terms for derivatives. This is done before FENDEQU,
+         because residuals of a subsequent block may depend on temporary terms for
+         the derivatives of the present block.
+
+         Also note that in the case of “evaluate” blocks, derivatives are not
+         computed in the “evaluate” mode; still their temporary terms must be
+         computed even in that mode, because for the same reason as above they may
+         be needed in subsequent blocks. */
+      write_eq_tt(blocks[block].size);
+
       FENDEQU_ fendequ;
       fendequ.write(code_file, instruction_number);
 
-      /* If the block is not of type “evaluate backward/forward”, then we write
-         the temporary terms for derivatives at this point, i.e. before the
-         JMPIFEVAL, because they will be needed in both “simulate” and
-         “evaluate” modes. */
-      if (simulation_type != BlockSimulationType::evaluateBackward
-          && simulation_type != BlockSimulationType::evaluateForward)
-        write_eq_tt(blocks[block].size);
-
       // Get the current code_file position and jump if eval = true
       streampos pos1 = code_file.tellp();
       FJMPIFEVAL_ fjmp_if_eval(0);
@@ -1345,13 +1348,6 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
       code_file.seekp(pos3);
       prev_instruction_number = instruction_number;
 
-      /* If the block is of type “evaluate backward/forward”, then write the
-         temporary terms for derivatives at this point, because they have not
-         been written before the JMPIFEVAL. */
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        write_eq_tt(blocks[block].size);
-
       // Write the derivatives for the “evaluate” mode
       for (const auto &[indices, d] : blocks_derivatives[block])
         {
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index ffea4fb8..7adf0a45 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -700,6 +700,17 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
               fstpr.write(code_file, instruction_number);
             }
         }
+
+      /* Write temporary terms for derivatives. This is done before FENDEQU,
+         because residuals of a subsequent block may depend on temporary terms for
+         the derivatives of the present block.
+
+         Also note that in the case of “evaluate” blocks, derivatives are not
+         computed in the “evaluate” mode; still their temporary terms must be
+         computed even in that mode, because for the same reason as above they may
+         be needed in subsequent blocks. */
+      write_eq_tt(blocks[block].size);
+
       FENDEQU_ fendequ;
       fendequ.write(code_file, instruction_number);
 
@@ -707,9 +718,6 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
       if (simulation_type != BlockSimulationType::evaluateBackward
           && simulation_type != BlockSimulationType::evaluateForward)
         {
-          // Write temporary terms for derivatives
-          write_eq_tt(blocks[block].size);
-
           switch (simulation_type)
             {
             case BlockSimulationType::solveBackwardSimple:
@@ -876,14 +884,16 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
               fstpr.write(code_file, instruction_number);
             }
         }
+
+      // Write temporary terms for derivatives
+      // Same remark as above
+      write_eq_tt(blocks[block].size);
+
       FENDEQU_ fendequ_l;
       fendequ_l.write(code_file, instruction_number);
 
       // The Jacobian if we have to solve the block determinsitic bloc
 
-      // Write temporary terms for derivatives
-      write_eq_tt(blocks[block].size);
-
       switch (simulation_type)
         {
         case BlockSimulationType::solveBackwardSimple:
-- 
GitLab