Derivation engine w.r.t. endogenous generalized to any order

src/ComputingTasks.cc

@@ -2167,7 +2167,7 @@ PlannerObjectiveStatement::getPlannerObjective() const
 void
 PlannerObjectiveStatement::computingPass()
 {
-  model_tree.computingPass({}, false, true, true, 0, false, false, false);
+  model_tree.computingPass(3, 0, {}, false, false, false, false);
   computing_pass_called = true;
 }
 

src/DynamicModel.cc

@@ -4367,11 +4367,11 @@ DynamicModel::substitutePacExpectation()
 }
 
 void
-DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
+DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
-                            bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition)
+                            bool bytecode, bool nopreprocessoroutput, bool linear_decomposition)
 {
-  assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder));
+  assert(jacobianExo || (derivsOrder < 2 && paramsDerivsOrder == 0));
 
   initializeVariablesAndEquations();
 
@@ -4393,38 +4393,18 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
   // Launch computations
   if (!nopreprocessoroutput)
-    {
-      if (linear_decomposition)
-        cout << "Computing nonlinear dynamic model derivatives:" << endl
-             << " - order 1" << endl;
-      else
-        cout << "Computing dynamic model derivatives:" << endl
-             << " - order 1" << endl;
-    }
-
-  computeJacobian(vars);
+    cout << "Computing " << (linear_decomposition ? "nonlinear " : "")
+         << "dynamic model derivatives (order " << derivsOrder << ")." << endl;
 
-  if (hessian)
-    {
-      if (!nopreprocessoroutput)
-        cout << " - order 2" << endl;
-      computeHessian(vars);
-    }
+  computeDerivatives(derivsOrder, vars);
 
   if (paramsDerivsOrder > 0)
     {
       if (!nopreprocessoroutput)
-        cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
+        cout << "Computing dynamic model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl;
       computeParamsDerivatives(paramsDerivsOrder);
     }
 
-  if (thirdDerivatives)
-    {
-      if (!nopreprocessoroutput)
-        cout << " - order 3" << endl;
-      computeThirdDerivatives(vars);
-    }
-
   jacob_map_t contemporaneous_jacobian, static_jacobian;
   map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives;
   // for each block contains pair<Size, Feddback_variable>

src/DynamicModel.hh

@@ -281,14 +281,13 @@ public:
   //! Execute computations (variable sorting + derivation)
   /*!
     \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
-    \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
-    \param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
+    \param derivsOrder order of derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true when order >= 2)
     \param paramsDerivsOrder order of derivatives w.r. to a pair (endo/exo/exo_det, parameter) to be computed (>0 implies jacobianExo = true)
     \param eval_context evaluation context for normalization
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
-  void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
-                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition);
+  void computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
+                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool nopreprocessoroutput, bool linear_decomposition);
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
 

src/ModFile.cc

@@ -701,7 +701,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
         {
           non_linear_equations_dynamic_model = dynamic_model;
           non_linear_equations_dynamic_model.set_cutoff_to_zero();
-          non_linear_equations_dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
+          non_linear_equations_dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
         }
       if (!no_static)
         {
@@ -711,13 +711,14 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
               || mod_file_struct.calib_smoother_present)
             static_model.set_cutoff_to_zero();
 
-          const bool static_hessian = mod_file_struct.identification_present
-            || mod_file_struct.estimation_analytic_derivation;
+          int derivsOrder = 1;
           int paramsDerivsOrder = 0;
           if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
-            paramsDerivsOrder = params_derivs_order;
-          static_model.computingPass(global_eval_context, no_tmp_terms, static_hessian,
-                                     false, paramsDerivsOrder, block, byte_code, nopreprocessoroutput);
+            {
+              derivsOrder = 2;
+              paramsDerivsOrder = params_derivs_order;
+            }
+          static_model.computingPass(derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, byte_code, nopreprocessoroutput);
         }
       // Set things to compute for dynamic model
       if (mod_file_struct.perfect_foresight_solver_present || mod_file_struct.check_present
@@ -727,7 +728,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
           || mod_file_struct.calib_smoother_present)
         {
           if (mod_file_struct.perfect_foresight_solver_present)
-            dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
+            dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
           else
             {
               if (mod_file_struct.stoch_simul_present
@@ -740,25 +741,21 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
                   cerr << "ERROR: Incorrect order option..." << endl;
                   exit(EXIT_FAILURE);
                 }
-              bool hessian = mod_file_struct.order_option >= 2
-                || mod_file_struct.identification_present
-                || mod_file_struct.estimation_analytic_derivation
-                || linear
-                || output == FileOutputType::second
-                || output == FileOutputType::third;
-              bool thirdDerivatives = mod_file_struct.order_option == 3
-                || mod_file_struct.estimation_analytic_derivation
-                || output == FileOutputType::third;
+              int derivsOrder = mod_file_struct.order_option;
+              if (mod_file_struct.identification_present || linear || output == FileOutputType::second)
+                derivsOrder = max(derivsOrder, 2);
+              if (mod_file_struct.estimation_analytic_derivation || output == FileOutputType::third)
+                derivsOrder = max(derivsOrder, 3);
               int paramsDerivsOrder = 0;
               if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
                 paramsDerivsOrder = params_derivs_order;
-              dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
+              dynamic_model.computingPass(true, derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
               if (linear && mod_file_struct.ramsey_model_present)
-                orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
+                orig_ramsey_dynamic_model.computingPass(true, 2, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
             }
         }
       else // No computing task requested, compute derivatives up to 2nd order by default
-        dynamic_model.computingPass(true, true, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
+        dynamic_model.computingPass(true, 2, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
 
       map<int, string> eqs;
       if (mod_file_struct.ramsey_model_present)
@@ -784,7 +781,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
   for (auto & statement : statements)
     statement->computingPass();
 
-  epilogue.computingPass(true, true, false, 0, global_eval_context, true, false, false, false, true, false);
+  epilogue.computingPass(true, 2, 0, global_eval_context, true, false, false, false, true, false);
 }
 
 void

src/ModelTree.cc

@@ -1276,77 +1276,45 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
 }
 
 void
-ModelTree::computeJacobian(const set<int> &vars)
+ModelTree::computeDerivatives(int order, const set<int> &vars)
 {
-  for (int var : vars)
-    {
-      for (int eq = 0; eq < (int) equations.size(); eq++)
-        {
-          expr_t d1 = equations[eq]->getDerivative(var);
-          if (d1 == Zero)
-            continue;
-          derivatives[1][{ eq, var }] = d1;
-          ++NNZDerivatives[1];
-        }
-    }
-}
-
-void
-ModelTree::computeHessian(const set<int> &vars)
-{
-  for (const auto &it : derivatives[1])
-    {
-      int eq, var1;
-      tie(eq, var1) = vectorToTuple<2>(it.first);
-      expr_t d1 = it.second;
-
-      // Store only second derivatives with var2 <= var1
-      for (int var2 : vars)
-        {
-          if (var2 > var1)
-            continue;
+  assert (order >= 1);
 
-          expr_t d2 = d1->getDerivative(var2);
-          if (d2 == Zero)
-            continue;
-          derivatives[2][{ eq, var1, var2 }] = d2;
-          if (var2 == var1)
-            ++NNZDerivatives[2];
-          else
-            NNZDerivatives[2] += 2;
-        }
-    }
-}
-
-void
-ModelTree::computeThirdDerivatives(const set<int> &vars)
-{
-  for (const auto &it : derivatives[2])
-    {
-      int eq, var1, var2;
-      tie(eq, var1, var2) = vectorToTuple<3>(it.first);
-      // By construction, var2 <= var1
+  // Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
+  derivatives.resize(max(static_cast<size_t>(order), derivatives.size()));
+  NNZDerivatives.resize(max(static_cast<size_t>(order), NNZDerivatives.size()), 0);
 
-      expr_t d2 = it.second;
+  // First-order derivatives
+  for (int var : vars)
+    for (int eq = 0; eq < (int) equations.size(); eq++)
+      {
+        expr_t d1 = equations[eq]->getDerivative(var);
+        if (d1 == Zero)
+          continue;
+        derivatives[1][{ eq, var }] = d1;
+        ++NNZDerivatives[1];
+      }
 
-      // Store only third derivatives such that var3 <= var2 <= var1
-      for (int var3 : vars)
+  // Higher-order derivatives
+  for (int o = 2; o <= order; o++)
+    for (const auto &it : derivatives[o-1])
+      for (int var : vars)
         {
-          if (var3 > var2)
+          if (it.first.back() > var)
             continue;
 
-          expr_t d3 = d2->getDerivative(var3);
-          if (d3 == Zero)
+          expr_t d = it.second->getDerivative(var);
+          if (d == Zero)
             continue;
-          derivatives[3][{ eq, var1, var2, var3 }] = d3;
-          if (var3 == var2 && var2 == var1)
-            ++NNZDerivatives[3];
-          else if (var3 == var2 || var2 == var1)
-            NNZDerivatives[3] += 3;
-          else
-            NNZDerivatives[3] += 6;
+
+          vector<int> indices{it.first};
+          indices.push_back(var);
+          // At this point, indices of endogenous variables are sorted in non-decreasing order
+          derivatives[o][indices] = d;
+          do
+            NNZDerivatives[o]++;
+          while (next_permutation(next(indices.begin()), indices.end()));
         }
-    }
 }
 
 void

src/ModelTree.hh

@@ -91,7 +91,8 @@ protected:
   /*! Index 0 is not used, index 1 contains first derivatives, ...
      For each derivation order, stores a map whose key is a vector of integer: the
      first integer is the equation index, the remaining ones are the derivation
-     IDs of variables */
+     IDs of variables (in non-decreasing order, to avoid storing symmetric
+     elements several times) */
   vector<map<vector<int>, expr_t>> derivatives;
 
   //! Number of non-zero derivatives
@@ -104,7 +105,7 @@ protected:
   derivation order w.r.t. parameters). For e.g., { 1, 2 } corresponds to the jacobian
   differentiated twice w.r.t. to parameters.
   In inner maps, the vector of integers consists of: the equation index, then
-  the derivation IDs of endogenous, then the IDs of parameters */
+  the derivation IDs of endogenous (in non-decreasing order), then the IDs of parameters */
   map<pair<int,int>, map<vector<int>, expr_t>> params_derivatives;
 
   //! Storage for temporary terms in block/bytecode mode
@@ -144,15 +145,10 @@ protected:
   //! Vector indicating if the equation is linear in endogenous variable (true) or not (false)
   vector<bool> is_equation_linear;
 
-  //! Computes 1st derivatives
-  /*! \param vars the derivation IDs w.r. to which compute the derivatives */
-  void computeJacobian(const set<int> &vars);
-  //! Computes 2nd derivatives
-  /*! \param vars the derivation IDs w.r. to which derive the 1st derivatives */
-  void computeHessian(const set<int> &vars);
-  //! Computes 3rd derivatives
-  /*! \param vars the derivation IDs w.r. to which derive the 2nd derivatives */
-  void computeThirdDerivatives(const set<int> &vars);
+  //! Computes derivatives
+  /*! \param order the derivation order
+      \param vars the derivation IDs w.r.t. which compute the derivatives */
+  void computeDerivatives(int order, const set<int> &vars);
   //! Computes derivatives of the Jacobian and Hessian w.r. to parameters
   void computeParamsDerivatives(int paramsDerivsOrder);
   //! Write derivative of an equation w.r. to a variable

src/StaticModel.cc

@@ -1214,7 +1214,7 @@ StaticModel::collect_first_order_derivatives_endogenous()
 }
 
 void
-StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, bool block, bool bytecode, const bool nopreprocessoroutput)
+StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool bytecode, bool nopreprocessoroutput)
 {
   initializeVariablesAndEquations();
 
@@ -1233,9 +1233,9 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
 
   equations.clear();
   copy(neweqs.begin(), neweqs.end(), back_inserter(equations));
-  // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
-  set<int> vars;
 
+  // Compute derivatives w.r. to all endogenous
+  set<int> vars;
   for (int i = 0; i < symbol_table.endo_nbr(); i++)
     {
       int id = symbol_table.getID(SymbolType::endogenous, i);
@@ -1245,29 +1245,14 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
 
   // Launch computations
   if (!nopreprocessoroutput)
-    cout << "Computing static model derivatives:" << endl
-         << " - order 1" << endl;
-
-  computeJacobian(vars);
-
-  if (hessian)
-    {
-      if (!nopreprocessoroutput)
-        cout << " - order 2" << endl;
-      computeHessian(vars);
-    }
+    cout << "Computing static model derivatives (order " << derivsOrder << ")." << endl;
 
-  if (thirdDerivatives)
-    {
-      if (!nopreprocessoroutput)
-        cout << " - order 3" << endl;
-      computeThirdDerivatives(vars);
-    }
+  computeDerivatives(derivsOrder, vars);
 
   if (paramsDerivsOrder > 0)
     {
       if (!nopreprocessoroutput)
-        cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
+        cout << "Computing static model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl;
       computeParamsDerivatives(paramsDerivsOrder);
     }
 

src/StaticModel.hh

@@ -193,10 +193,10 @@ public:
   /*!
     \param eval_context evaluation context for normalization
     \param no_tmp_terms if true, no temporary terms will be computed in the static files
-    \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed
-    \param paramsDerivsOrder order of derivatives w.r. to a pair (endo/exo/exo_det, parameter) to be computed
+    \param derivsOrder order of derivation with respect to endogenous
+    \param paramsDerivsOrder order of derivatives w.r. to a pair (endogenous, parameter) to be computed
   */
-  void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatices, int paramsDerivsOrder, bool block, bool bytecode, const bool nopreprocessoroutput);
+  void computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool bytecode, bool nopreprocessoroutput);
 
   //! Adds informations for simulation in a binary file for a block decomposed model
   void Write_Inf_To_Bin_File_Block(const string &basename, const int &num,