Refactor methods for writing static and dynamic files

– factorize common code between the static and the dynamic version
– reorganise language-specific code into dedicated functions
– use a function template in the main helper to do some computations
  at compile-time (using constexpr features)

src/DataTree.cc

@@ -832,12 +832,6 @@ DataTree::addAllParamDerivId([[maybe_unused]] set<int> &deriv_id_set)
 {
 }
 
-int
-DataTree::getDynJacobianCol([[maybe_unused]] int deriv_id) const noexcept(false)
-{
-  throw UnknownDerivIDException();
-}
-
 bool
 DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
 {

src/DataTree.hh

@@ -305,8 +305,21 @@ public:
   virtual SymbolType getTypeByDerivID(int deriv_id) const noexcept(false);
   virtual int getLagByDerivID(int deriv_id) const noexcept(false);
   virtual int getSymbIDByDerivID(int deriv_id) const noexcept(false);
-  //! Returns the column of the dynamic Jacobian associated to a derivation ID
+
-  virtual int getDynJacobianCol(int deriv_id) const noexcept(false);
+  //! Returns the column of the Jacobian associated to a derivation ID
+  virtual int
+  getJacobianCol([[maybe_unused]] int deriv_id) const
+  {
+    throw UnknownDerivIDException();
+  }
+
+  //! Returns the number of columns of the Jacobian
+  virtual int
+  getJacobianColsNbr() const
+  {
+    throw UnknownDerivIDException();
+  }
+
   //! Adds to the set all the deriv IDs corresponding to parameters
   virtual void addAllParamDerivId(set<int> &deriv_id_set);
 

src/DynamicModel.hh

@@ -96,10 +96,6 @@ private:
   //! Nonzero equations in the Hessian
   set<int> nonzero_hessian_eqs;
 
-  //! Number of columns of dynamic jacobian
-  /*! Set by computeDerivID()s and computeDynJacobianCols() */
-  int dynJacobianColsNbr{0};
-
   //! Creates mapping for variables and equations they are present in
   map<int, set<int>> variableMapping;
 
@@ -124,15 +120,10 @@ private:
   //! Writes dynamic model file (Matlab version)
   void writeDynamicMFile(const string &basename) const;
   //! Writes dynamic model file (Julia version)
-  void writeDynamicJuliaFile(const string &dynamic_basename) const;
+  void writeDynamicJuliaFile(const string &basename) const;
   //! Writes dynamic model file (C version)
   /*! \todo add third derivatives handling */
   void writeDynamicCFile(const string &basename) const;
-  //! Writes the dynamic model equations and its derivatives
-  /*! \todo add third derivatives handling in C output */
-  void writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
-  void writeDynamicModel(const string &basename, bool use_dll, bool julia) const;
-  void writeDynamicModel(const string &basename, ostream &DynamicOutput, bool use_dll, bool julia) const;
   //! Writes the main dynamic function of block decomposed model (MATLAB version)
   void writeDynamicBlockMFile(const string &basename) const;
   //! Writes the main dynamic function of block decomposed model (C version)
@@ -463,7 +454,22 @@ public:
   void writeLatexOriginalFile(const string &basename, bool write_equation_tags) const;
 
   int getDerivID(int symb_id, int lag) const noexcept(false) override;
-  int getDynJacobianCol(int deriv_id) const noexcept(false) override;
+
+  int
+  getJacobianCol(int deriv_id) const override
+  {
+    if (auto it = dyn_jacobian_cols_table.find(deriv_id);
+        it == dyn_jacobian_cols_table.end())
+      throw UnknownDerivIDException();
+    else
+      return it->second;
+  }
+  int
+  getJacobianColsNbr() const override
+  {
+    return dyn_jacobian_cols_table.size();
+  }
+
   void addAllParamDerivId(set<int> &deriv_id_set) override;
 
   //! Returns true indicating that this is a dynamic model

src/ExprNode.cc

@@ -1074,7 +1074,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case ExprNodeOutputType::juliaDynamicModel:
         case ExprNodeOutputType::matlabDynamicModel:
         case ExprNodeOutputType::CDynamicModel:
-          i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
+          i = datatree.getJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
           output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case ExprNodeOutputType::CStaticModel:

src/ModelTree.cc

@@ -1326,14 +1326,6 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, d
   output << "]";
 }
 
-void
-ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
-{
-  temporary_terms_t tt;
-  temporary_terms_idxs_t ttidxs;
-  writeModelEquations(output, output_type, tt);
-}
-
 void
 ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
                                const temporary_terms_t &temporary_terms) const
@@ -1599,19 +1591,6 @@ ModelTree::set_cutoff_to_zero()
   cutoff = 0;
 }
 
-void
-ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
-{
-  if (isJuliaOutput(output_type))
-    output << "    @inbounds ";
-  output << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (isMatlabOutput(output_type) || isJuliaOutput(output_type))
-    output << eq_nb + 1 << "," << col_nb + 1;
-  else
-    output << eq_nb + col_nb *equations.size();
-  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
-}
-
 void
 ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
 {

src/ModelTree.hh

@@ -28,6 +28,7 @@
 #include <array>
 #include <filesystem>
 #include <optional>
+#include <cassert>
 
 #include "DataTree.hh"
 #include "EquationTags.hh"
@@ -254,9 +255,13 @@ protected:
                                              ostream &output, ExprNodeOutputType output_type,
                                              deriv_node_temp_terms_t &tef_terms) const;
   //! Writes model equations
-  void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const;
   void writeModelEquations(ostream &output, ExprNodeOutputType output_type,
                            const temporary_terms_t &temporary_terms) const;
+
+  // Returns outputs for derivatives and temporary terms at each derivation order
+  template<ExprNodeOutputType output_type>
+  pair<vector<ostringstream>, vector<ostringstream>> writeModelFileHelper() const;
+
   //! Writes JSON model equations
   //! if residuals = true, we are writing the dynamic/static model.
   //! Otherwise, just the model equations (with line numbers, no tmp terms)
@@ -448,9 +453,6 @@ public:
   void reorderAuxiliaryEquations();
   //! Find equations of the form “variable=constant”, excluding equations with “mcp” tag (see dynare#1697)
   void findConstantEquationsWithoutMcpTag(map<VariableNode *, NumConstNode *> &subst_table) const;
-  //! Helper for writing the Jacobian elements in MATLAB and C
-  /*! Writes either (i+1,j+1) or [i+j*no_eq] */
-  void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
   /* Given an expression, searches for the first equation that has exactly this
      expression on the LHS, and returns the RHS of that equation.
      If no such equation can be found, throws an ExprNode::MatchFailureExpression */
@@ -511,4 +513,139 @@ public:
   }
 };
 
+template<ExprNodeOutputType output_type>
+pair<vector<ostringstream>, vector<ostringstream>>
+ModelTree::writeModelFileHelper() const
+{
+  vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
+  vector<ostringstream> tt_output(derivatives.size()); // Temp terms output (at all orders)
+
+  deriv_node_temp_terms_t tef_terms;
+  temporary_terms_t temp_term_union;
+
+  writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs,
+                                        tt_output[0], output_type, tef_terms);
+
+  writeTemporaryTerms(temporary_terms_derivatives[0],
+                      temp_term_union,
+                      temporary_terms_idxs,
+                      tt_output[0], output_type, tef_terms);
+
+  writeModelEquations(d_output[0], output_type, temp_term_union);
+
+  // Writing Jacobian
+  if (!derivatives[1].empty())
+    {
+      writeTemporaryTerms(temporary_terms_derivatives[1],
+                          temp_term_union,
+                          temporary_terms_idxs,
+                          tt_output[1], output_type, tef_terms);
+
+      for (const auto &[indices, d1] : derivatives[1])
+        {
+          auto [eq, var] = vectorToTuple<2>(indices);
+
+          if constexpr(isJuliaOutput(output_type))
+            d_output[1] << "    @inbounds ";
+          d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
+          if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type))
+            d_output[1] << eq + 1 << "," << getJacobianCol(var) + 1;
+          else
+            d_output[1] << eq + getJacobianCol(var)*equations.size();
+          d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+          d1->writeOutput(d_output[1], output_type,
+                          temp_term_union, temporary_terms_idxs, tef_terms);
+          d_output[1] << ";" << endl;
+        }
+    }
+
+  // Write derivatives for order ≥ 2
+  for (size_t i = 2; i < derivatives.size(); i++)
+    if (!derivatives[i].empty())
+      {
+        writeTemporaryTerms(temporary_terms_derivatives[i],
+                            temp_term_union,
+                            temporary_terms_idxs,
+                            tt_output[i], output_type, tef_terms);
+
+        /* When creating the sparse matrix (in MATLAB or C mode), since storage
+           is in column-major order, output the first column, then the second,
+           then the third. This gives a significant performance boost in use_dll
+           mode (at both compilation and runtime), because it facilitates memory
+           accesses and expression reusage. */
+        ostringstream i_output, j_output, v_output;
+
+        for (int k{0}; // Current line index in the 3-column matrix
+             const auto &[vidx, d] : derivatives[i])
+          {
+            int eq{vidx[0]};
+
+            int col_idx{0};
+            for (size_t j = 1; j < vidx.size(); j++)
+              {
+                col_idx *= getJacobianColsNbr();
+                col_idx += getJacobianCol(vidx[j]);
+              }
+
+            if constexpr(isJuliaOutput(output_type))
+              {
+                d_output[i] << "    @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = ";
+                d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms);
+                d_output[i] << endl;
+              }
+            else
+              {
+                i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                         << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                         << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                         << "=" << eq + 1 << ";" << endl;
+                j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                         << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                         << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                         << "=" << col_idx + 1 << ";" << endl;
+                v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                         << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+                d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
+                v_output << ";" << endl;
+
+                k++;
+              }
+
+            // Output symetric elements at order 2
+            if (i == 2 && vidx[1] != vidx[2])
+              {
+                int col_idx_sym{getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1])};
+
+                if constexpr(isJuliaOutput(output_type))
+                  d_output[2] << "    @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
+                              << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
+                else
+                  {
+                    i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                             << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                             << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                             << "=" << eq + 1 << ";" << endl;
+                    j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                             << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                             << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                             << "=" << col_idx_sym + 1 << ";" << endl;
+                    v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                             << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                             << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="
+                             << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                             << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                             << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
+
+                    k++;
+                  }
+              }
+          }
+        if constexpr(!isJuliaOutput(output_type))
+          d_output[i] << i_output.str() << j_output.str() << v_output.str();
+      }
+
+  return { move(d_output), move(tt_output) };
+}
+
 #endif

src/StaticModel.cc

@@ -939,7 +939,73 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
 void
 StaticModel::writeStaticMFile(const string &basename) const
 {
-  writeStaticModel(basename, false, false);
+  auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::matlabStaticModel>();
+
+  // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
+  map<string, string> tmp_paren_vars;
+  bool message_printed{false};
+  for (auto &it : tt_output)
+    fixNestedParenthesis(it, tmp_paren_vars, message_printed);
+  for (auto &it : d_output)
+    fixNestedParenthesis(it, tmp_paren_vars, message_printed);
+
+  ostringstream init_output, end_output;
+  init_output << "residual = zeros(" << equations.size() << ", 1);";
+  end_output << "if ~isreal(residual)" << endl
+             << "  residual = real(residual)+imag(residual).^2;" << endl
+             << "end";
+  writeStaticModelHelper(basename, "static_resid", "residual", "static_resid_tt",
+                         temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(),
+                         "", init_output, end_output,
+                         d_output[0], tt_output[0]);
+
+  init_output.str("");
+  end_output.str("");
+  init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");";
+  end_output << "if ~isreal(g1)" << endl
+             << "    g1 = real(g1)+2*imag(g1);" << endl
+             << "end";
+  writeStaticModelHelper(basename, "static_g1", "g1", "static_g1_tt",
+                         temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
+                         "static_resid_tt",
+                         init_output, end_output,
+                         d_output[1], tt_output[1]);
+  writeWrapperFunctions(basename, "g1");
+
+  // For order ≥ 2
+  int ncols{symbol_table.endo_nbr()};
+  int ntt{static_cast<int>(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())};
+  for (size_t i{2}; i < derivatives.size(); i++)
+    {
+      ncols *= symbol_table.endo_nbr();
+      ntt += temporary_terms_derivatives[i].size();
+      string gname{"g" + to_string(i)};
+      string gprevname{"g" + to_string(i-1)};
+
+      init_output.str("");
+      end_output.str("");
+      if (derivatives[i].size())
+        {
+          init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
+                      << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
+                      << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
+          end_output << gname << " = sparse("
+                     << gname << "_i," << gname << "_j," << gname << "_v,"
+                     << equations.size() << "," << ncols << ");";
+        }
+      else
+        init_output << gname << " = sparse([],[],[]," << equations.size() << "," << ncols << ");";
+      writeStaticModelHelper(basename, "static_" + gname, gname,
+                             "static_" + gname + "_tt",
+                             ntt,
+                             "static_" + gprevname + "_tt",
+                             init_output, end_output,
+                             d_output[i], tt_output[i]);
+      if (i <= 3)
+        writeWrapperFunctions(basename, gname);
+    }
+
+  writeStaticMatlabCompatLayer(basename);
 }
 
 void
@@ -1099,249 +1165,121 @@ StaticModel::writeStaticMatlabCompatLayer(const string &basename) const
 }
 
 void
-StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const
+StaticModel::writeStaticCFile(const string &basename) const
 {
-  writeStaticModel("", StaticOutput, use_dll, julia);
+  // Writing comments and function definition command
-}
+  string filename{basename + "/model/src/static.c"};
 
-void
+  int ntt{static_cast<int>(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size())};
-StaticModel::writeStaticModel(const string &basename, bool use_dll, bool julia) const
-{
-  ofstream StaticOutput;
-  writeStaticModel(basename, StaticOutput, use_dll, julia);
-}
 
-void
+  ofstream output{filename, ios::out | ios::binary};
-StaticModel::writeStaticModel(const string &basename,
+  if (!output.is_open())
-                              ostream &StaticOutput, bool use_dll, bool julia) const
-{
-  vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
-  vector<ostringstream> tt_output(derivatives.size()); // Temp terms output (at all orders)
-
-  ExprNodeOutputType output_type = (use_dll ? ExprNodeOutputType::CStaticModel :
-                                    julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel);
-
-  deriv_node_temp_terms_t tef_terms;
-  temporary_terms_t temp_term_union;
-
-  writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs,
-                                        tt_output[0], output_type, tef_terms);
-
-  writeTemporaryTerms(temporary_terms_derivatives[0],
-                      temp_term_union,
-                      temporary_terms_idxs,
-                      tt_output[0], output_type, tef_terms);
-
-  writeModelEquations(d_output[0], output_type, temp_term_union);
-
-  int nrows = equations.size();
-  int JacobianColsNbr = symbol_table.endo_nbr();
-  int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
-
-  auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); };
-
-  // Write Jacobian w.r. to endogenous only
-  if (!derivatives[1].empty())
     {
-      writeTemporaryTerms(temporary_terms_derivatives[1],
+      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
-                          temp_term_union,
+      exit(EXIT_FAILURE);
-                          temporary_terms_idxs,
-                          tt_output[1], output_type, tef_terms);
-
-      for (const auto & [indices, d1] : derivatives[1])
-        {
-          auto [eq, var] = vectorToTuple<2>(indices);
-
-          jacobianHelper(d_output[1], eq, getJacobCol(var), output_type);
-          d_output[1] << "=";
-          d1->writeOutput(d_output[1], output_type,
-                          temp_term_union, temporary_terms_idxs, tef_terms);
-          d_output[1] << ";" << endl;
-        }
     }
 
-  // Write derivatives for order ≥ 2
+  output << "/*" << endl
-  for (size_t i = 2; i < derivatives.size(); i++)
+         << " * " << filename << " : Computes static model for Dynare" << endl
-    if (!derivatives[i].empty())
+         << " *" << endl
-      {
+         << " * Warning : this file is generated automatically by Dynare" << endl
-        writeTemporaryTerms(temporary_terms_derivatives[i],
+         << " *           from model file (.mod)" << endl << endl
-                            temp_term_union,
+         << " */" << endl
-                            temporary_terms_idxs,
+         << endl
-                            tt_output[i], output_type, tef_terms);
+         << "#include <math.h>" << endl
+         << "#include <stdlib.h>" << endl
+         << R"(#include "mex.h")" << endl
+         << endl;
 
-        /* When creating the sparse matrix (in MATLAB or C mode), since storage
+  // Write function definition if BinaryOpcode::powerDeriv is used
-           is in column-major order, output the first column, then the second,
+  writePowerDeriv(output);
-           then the third. This gives a significant performance boost in use_dll
-           mode (at both compilation and runtime), because it facilitates memory
-           accesses and expression reusage. */
-        ostringstream i_output, j_output, v_output;
 
-        for (int k{0}; // Current line index in the 3-column matrix
+  output << endl;
-             const auto &[vidx, d] : derivatives[i])
-          {
-            int eq = vidx[0];
 
-            int col_idx = 0;
+  auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::CStaticModel>();
-            for (size_t j = 1; j < vidx.size(); j++)
-              {
-                col_idx *= JacobianColsNbr;
-                col_idx += getJacobCol(vidx[j]);
-              }
 
-            if (output_type == ExprNodeOutputType::juliaStaticModel)
-              {
-                d_output[i] << "    @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = ";
-                d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms);
-                d_output[i] << endl;
-              }
-            else
-              {
-                i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                         << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                         << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                         << "=" << eq + 1 << ";" << endl;
-                j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                         << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                         << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                         << "=" << col_idx + 1 << ";" << endl;
-                v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                         << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-                d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
-                v_output << ";" << endl;
-
-                k++;
-              }
-
-            // Output symetric elements at order 2
-            if (i == 2 && vidx[1] != vidx[2])
-              {
-                int col_idx_sym = getJacobCol(vidx[2]) * JacobianColsNbr + getJacobCol(vidx[1]);
-
-                if (output_type == ExprNodeOutputType::juliaStaticModel)
-                  d_output[2] << "    @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
-                              << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
-                else
-                  {
-                    i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                             << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                             << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                             << "=" << eq + 1 << ";" << endl;
-                    j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                             << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                             << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                             << "=" << col_idx_sym + 1 << ";" << endl;
-                    v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                             << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                             << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="
-                             << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                             << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                             << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
-
-                    k++;
-                  }
-              }
-          }
-        if (output_type != ExprNodeOutputType::juliaStaticModel)
-          d_output[i] << i_output.str() << j_output.str() << v_output.str();
-      }
-
-  if (output_type == ExprNodeOutputType::matlabStaticModel)
-    {
-      // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
-      map<string, string> tmp_paren_vars;
-      bool message_printed = false;
-      for (auto &it : tt_output)
-        fixNestedParenthesis(it, tmp_paren_vars, message_printed);
-      for (auto &it : d_output)
-        fixNestedParenthesis(it, tmp_paren_vars, message_printed);
-
-      ostringstream init_output, end_output;
-      init_output << "residual = zeros(" << equations.size() << ", 1);";
-      end_output << "if ~isreal(residual)" << endl
-                 << "  residual = real(residual)+imag(residual).^2;" << endl
-                 << "end";
-      writeStaticModelHelper(basename, "static_resid", "residual", "static_resid_tt",
-                             temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(),
-                             "", init_output, end_output,
-                             d_output[0], tt_output[0]);
-
-      init_output.str("");
-      end_output.str("");
-      init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");";
-      end_output << "if ~isreal(g1)" << endl
-                 << "    g1 = real(g1)+2*imag(g1);" << endl
-                 << "end";
-      writeStaticModelHelper(basename, "static_g1", "g1", "static_g1_tt",
-                             temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
-                             "static_resid_tt",
-                             init_output, end_output,
-                             d_output[1], tt_output[1]);
-      writeWrapperFunctions(basename, "g1");
-
-      // For order ≥ 2
-      int ncols = JacobianColsNbr;
-      int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size();
-      for (size_t i = 2; i < derivatives.size(); i++)
-        {
-          ncols *= JacobianColsNbr;
-          ntt += temporary_terms_derivatives[i].size();
-          string gname = "g" + to_string(i);
-          string gprevname = "g" + to_string(i-1);
-
-          init_output.str("");
-          end_output.str("");
-          if (derivatives[i].size())
-            {
-              init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
-                          << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
-                          << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
-              end_output << gname << " = sparse("
-                         << gname << "_i," << gname << "_j," << gname << "_v,"
-                         << nrows << "," << ncols << ");";
-            }
-          else
-            init_output << gname << " = sparse([],[],[]," << nrows << "," << ncols << ");";
-          writeStaticModelHelper(basename, "static_" + gname, gname,
-                                 "static_" + gname + "_tt",
-                                 ntt,
-                                 "static_" + gprevname + "_tt",
-                                 init_output, end_output,
-                                 d_output[i], tt_output[i]);
-          if (i <= 3)
-            writeWrapperFunctions(basename, gname);
-        }
-
-      writeStaticMatlabCompatLayer(basename);
-    }
-  else if (output_type == ExprNodeOutputType::CStaticModel)
-    {
   for (size_t i = 0; i < d_output.size(); i++)
     {
-          string funcname = i == 0 ? "resid" : "g" + to_string(i);
+      string funcname{i == 0 ? "resid" : "g" + to_string(i)};
-          StaticOutput << "void static_" << funcname << "_tt(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl
+      output << "void static_" << funcname << "_tt(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl
              << "{" << endl
              << tt_output[i].str()
              << "}" << endl
              << endl
              << "void static_" << funcname << "(const double *restrict y, const double *restrict x, const double *restrict params, const double *restrict T, ";
       if (i == 0)
-            StaticOutput << "double *restrict residual";
+        output << "double *restrict residual";
       else if (i == 1)
-            StaticOutput << "double *restrict g1";
+        output << "double *restrict g1";
       else
-            StaticOutput << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v";
+        output << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v";
-          StaticOutput << ")" << endl
+      output << ")" << endl
                    << "{" << endl;
       if (i == 0)
-            StaticOutput << "  double lhs, rhs;" << endl;
+        output << "  double lhs, rhs;" << endl;
-          StaticOutput << d_output[i].str()
+      output << d_output[i].str()
              << "}" << endl
              << endl;
     }
-    }
+
-  else
+  output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
-    {
+         << "{" << endl
+         << "  if (nrhs > 3)" << endl
+         << R"(    mexErrMsgTxt("Accepts at most 3 output arguments");)" << endl
+         << "  if (nrhs != 3)" << endl
+         << R"(    mexErrMsgTxt("Requires exactly 3 input arguments");)" << endl
+         << "  double *y = mxGetPr(prhs[0]);" << endl
+         << "  double *x = mxGetPr(prhs[1]);" << endl
+         << "  double *params = mxGetPr(prhs[2]);" << endl
+         << endl
+         << "  double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
+         << endl
+         << "  if (nlhs >= 1)" << endl
+         << "    {" << endl
+         << "      plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
+         << "      double *residual = mxGetPr(plhs[0]);" << endl
+         << "      static_resid_tt(y, x, params, T);" << endl
+         << "      static_resid(y, x, params, T, residual);" << endl
+         << "    }" << endl
+         << endl
+         << "  if (nlhs >= 2)" << endl
+         << "    {" << endl
+         << "      plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
+         << "      double *g1 = mxGetPr(plhs[1]);" << endl
+         << "      static_g1_tt(y, x, params, T);" << endl
+         << "      static_g1(y, x, params, T, g1);" << endl
+         << "    }" << endl
+         << endl
+         << "  if (nlhs >= 3)" << endl
+         << "    {" << endl
+         << "      mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+         << "      mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+         << "      mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+         << "      static_g2_tt(y, x, params, T);" << endl
+         << "      static_g2(y, x, params, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
+         << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
+         << "      mxArray *n = mxCreateDoubleScalar(" << symbol_table.endo_nbr()*symbol_table.endo_nbr() << ");" << endl
+         << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
+         << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
+         << "      plhs[2] = plhs_sparse[0];" << endl
+         << "      mxDestroyArray(g2_i);" << endl
+         << "      mxDestroyArray(g2_j);" << endl
+         << "      mxDestroyArray(g2_v);" << endl
+         << "      mxDestroyArray(m);" << endl
+         << "      mxDestroyArray(n);" << endl
+         << "    }" << endl
+         << endl
+         << "  free(T);" << endl
+         << "}" << endl;
+
+  output.close();
+}
+
+void
+StaticModel::writeStaticJuliaFile(const string &basename) const
+{
+  auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::juliaStaticModel>();
+
   stringstream output;
   output << "module " << basename << "Static" << endl
          << "#" << endl
@@ -1467,6 +1405,7 @@ StaticModel::writeStaticModel(const string &basename,
          << "end" << endl << endl;
 
   // staticG2!
+  int hessianColsNbr{symbol_table.endo_nbr() * symbol_table.endo_nbr()};
   output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
          << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    @assert length(T) >= "
@@ -1494,12 +1433,12 @@ StaticModel::writeStaticModel(const string &basename,
          << "end" << endl << endl;
 
   // staticG3!
-      int ncols = hessianColsNbr * JacobianColsNbr;
+  int ncols{hessianColsNbr * symbol_table.endo_nbr()};
   output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
          << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    @assert length(T) >= "
          << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl
-             << "    @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+         << "    @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl
          << "    @assert length(y) == " << symbol_table.endo_nbr() << endl
          << "    @assert length(x) == " << symbol_table.exo_nbr() << endl
          << "    @assert length(params) == " << symbol_table.param_nbr() << endl
@@ -1546,100 +1485,6 @@ StaticModel::writeStaticModel(const string &basename,
   output << "end" << endl;
 
   writeToFileIfModified(output, basename + "Static.jl");
-    }
-}
-
-void
-StaticModel::writeStaticCFile(const string &basename) const
-{
-  // Writing comments and function definition command
-  string filename = basename + "/model/src/static.c";
-
-  int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
-
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "/*" << endl
-         << " * " << filename << " : Computes static model for Dynare" << endl
-         << " *" << endl
-         << " * Warning : this file is generated automatically by Dynare" << endl
-         << " *           from model file (.mod)" << endl << endl
-         << " */" << endl
-         << endl
-         << "#include <math.h>" << endl
-         << "#include <stdlib.h>" << endl
-         << R"(#include "mex.h")" << endl
-         << endl;
-
-  // Write function definition if BinaryOpcode::powerDeriv is used
-  writePowerDeriv(output);
-
-  output << endl;
-
-  writeStaticModel(output, true, false);
-
-  output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
-         << "{" << endl
-         << "  if (nrhs > 3)" << endl
-         << R"(    mexErrMsgTxt("Accepts at most 3 output arguments");)" << endl
-         << "  if (nrhs != 3)" << endl
-         << R"(    mexErrMsgTxt("Requires exactly 3 input arguments");)" << endl
-         << "  double *y = mxGetPr(prhs[0]);" << endl
-         << "  double *x = mxGetPr(prhs[1]);" << endl
-         << "  double *params = mxGetPr(prhs[2]);" << endl
-         << endl
-         << "  double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
-         << endl
-         << "  if (nlhs >= 1)" << endl
-         << "    {" << endl
-         << "      plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
-         << "      double *residual = mxGetPr(plhs[0]);" << endl
-         << "      static_resid_tt(y, x, params, T);" << endl
-         << "      static_resid(y, x, params, T, residual);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 2)" << endl
-         << "    {" << endl
-         << "      plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
-         << "      double *g1 = mxGetPr(plhs[1]);" << endl
-         << "      static_g1_tt(y, x, params, T);" << endl
-         << "      static_g1(y, x, params, T, g1);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 3)" << endl
-         << "    {" << endl
-         << "      mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      static_g2_tt(y, x, params, T);" << endl
-         << "      static_g2(y, x, params, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
-         << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
-         << "      mxArray *n = mxCreateDoubleScalar(" << symbol_table.endo_nbr()*symbol_table.endo_nbr() << ");" << endl
-         << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
-         << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
-         << "      plhs[2] = plhs_sparse[0];" << endl
-         << "      mxDestroyArray(g2_i);" << endl
-         << "      mxDestroyArray(g2_j);" << endl
-         << "      mxDestroyArray(g2_v);" << endl
-         << "      mxDestroyArray(m);" << endl
-         << "      mxDestroyArray(n);" << endl
-         << "    }" << endl
-         << endl
-         << "  free(T);" << endl
-         << "}" << endl;
-
-  output.close();
-}
-
-void
-StaticModel::writeStaticJuliaFile(const string &basename) const
-{
-  writeStaticModel(basename, false, true);
 }
 
 void

src/StaticModel.hh

@@ -43,9 +43,6 @@ private:
   //! Writes static model file (Julia version)
   void writeStaticJuliaFile(const string &basename) const;
 
-  //! Writes the static model equations and its derivatives
-  void writeStaticModel(const string &basename, ostream &StaticOutput, bool use_dll, bool julia) const;
-
   //! Writes the main static function of block decomposed model (MATLAB version)
   void writeStaticBlockMFile(const string &basename) const;
 
@@ -89,8 +86,18 @@ private:
   int getLagByDerivID(int deriv_id) const noexcept(false) override;
   //! Get the symbol ID corresponding to a derivation ID
   int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
-  //! Compute the column indices of the static Jacobian
+
-  void computeStatJacobianCols();
+  int
+  getJacobianCol(int deriv_id) const override
+  {
+    return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id));
+  }
+  int
+  getJacobianColsNbr() const override
+  {
+    return symbol_table.endo_nbr();
+  }
+
   //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
   void computeChainRuleJacobian();
 
@@ -106,9 +113,6 @@ private:
   //! Create a legacy *_static.m file for Matlab/Octave not yet using the temporary terms array interface
   void writeStaticMatlabCompatLayer(const string &name) const;
 
-  void writeStaticModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
-  void writeStaticModel(const string &dynamic_basename, bool use_dll, bool julia) const;
-
 public:
   StaticModel(SymbolTable &symbol_table_arg,
               NumericalConstants &num_constants,