From c628f2124572d887183fb8d7aa9c83c8aa6ec1ed Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 18 Apr 2019 17:07:55 +0200
Subject: [PATCH] JSON: output derivatives at an arbitrary order

Backward incompatible change: the temporary terms for 3rd order are now stored
in "temporary_terms_third_derivative" (without the final "s"; same for external
functions), for consistency with the name of the slot for the derivatives
themselves ("third_derivative").

Ref dynare#217
---
 src/DynamicModel.cc | 194 ++++++++++++++----------------------------
 src/StaticModel.cc  | 201 ++++++++++++++------------------------------
 2 files changed, 127 insertions(+), 268 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index e7e62804..9fef7ade 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -6779,163 +6779,93 @@ void
 DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) const
 {
   ostringstream model_local_vars_output;  // Used for storing model local vars
-  ostringstream model_output;             // Used for storing model temp vars and equations
-  ostringstream jacobian_output;          // Used for storing jacobian equations
-  ostringstream hessian_output;           // Used for storing Hessian equations
-  ostringstream third_derivatives_output; // Used for storing third order derivatives equations
+  vector<ostringstream> d_output(derivatives.size());  // Derivatives output (at all orders, including 0=residual)
 
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temp_term_union;
 
-  int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
-
   writeJsonModelLocalVariables(model_local_vars_output, tef_terms);
 
-  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, model_output, tef_terms, "");
-  model_output << ", ";
-  writeJsonModelEquations(model_output, true);
+  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
+  d_output[0] << ", ";
+  writeJsonModelEquations(d_output[0], true);
 
-  // Writing Jacobian
-  writeJsonTemporaryTerms(temporary_terms_derivatives[1], temp_term_union, jacobian_output, tef_terms, "jacobian");
-  jacobian_output << R"(, "jacobian": {)"
+  int ncols = dynJacobianColsNbr;
+  for (size_t i = 1; i < derivatives.size(); i++)
+    {
+      string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative";
+      writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
+      temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
+      d_output[i] << R"(, ")" << matrix_name  << R"(": {)"
                   << R"(  "nrows": )" << equations.size()
-                  << R"(, "ncols": )" << dynJacobianColsNbr
+                  << R"(, "ncols": )" << ncols
                   << R"(, "entries": [)";
-  for (auto it = derivatives[1].begin();
-       it != derivatives[1].end(); it++)
-    {
-      if (it != derivatives[1].begin())
-        jacobian_output << ", ";
-
-      int eq, var;
-      tie(eq, var) = vectorToTuple<2>(it->first);
-      int col =  getDynJacobianCol(var);
-      expr_t d1 = it->second;
-
-      if (writeDetails)
-        jacobian_output << R"({"eq": )" << eq + 1;
-      else
-        jacobian_output << R"({"row": )" << eq + 1;
 
-      jacobian_output << R"(, "col": )" << col + 1;
-
-      if (writeDetails)
-        jacobian_output << R"(, "var": ")" << symbol_table.getName(getSymbIDByDerivID(var)) << R"(")"
-                        << R"(, "shift": )" << getLagByDerivID(var);
-
-      jacobian_output << R"(, "val": ")";
-      d1->writeJsonOutput(jacobian_output, temp_term_union, tef_terms);
-      jacobian_output << R"("})" << endl;
-    }
-  jacobian_output << "]}";
-
-  // Writing Hessian
-  writeJsonTemporaryTerms(temporary_terms_derivatives[2], temp_term_union, hessian_output, tef_terms, "hessian");
-  hessian_output << R"(, "hessian": {)"
-                 << R"(  "nrows": )" << equations.size()
-                 << R"(, "ncols": )" << hessianColsNbr
-                 << R"(, "entries": [)";
-  for (auto it = derivatives[2].begin();
-       it != derivatives[2].end(); it++)
-    {
-      if (it != derivatives[2].begin())
-        hessian_output << ", ";
+      for (auto it = derivatives[i].begin(); it != derivatives[i].end(); ++it)
+        {
+          if (it != derivatives[i].begin())
+            d_output[i] << ", ";
 
-      int eq, var1, var2;
-      tie(eq, var1, var2) = vectorToTuple<3>(it->first);
-      expr_t d2 = it->second;
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
-      int col_nb = id1 * dynJacobianColsNbr + id2;
-      int col_nb_sym = id2 * dynJacobianColsNbr + id1;
+          const vector<int> &vidx = it->first;
+          expr_t d = it->second;
+          int eq = vidx[0];
 
-      if (writeDetails)
-        hessian_output << R"({"eq": )" << eq + 1;
-      else
-        hessian_output << R"({"row": )" << eq + 1;
+          int col_idx = 0;
+          for (size_t j = 1; j < vidx.size(); j++)
+            {
+              col_idx *= dynJacobianColsNbr;
+              col_idx += getDynJacobianCol(vidx[j]);
+            }
 
-      hessian_output << R"(, "col": [)" << col_nb + 1;
-      if (id1 != id2)
-        hessian_output << ", " << col_nb_sym + 1;
-      hessian_output << "]";
+          if (writeDetails)
+            d_output[i] << R"({"eq": )" << eq + 1;
+          else
+            d_output[i] << R"({"row": )" << eq + 1;
 
-      if (writeDetails)
-        hessian_output << R"(, "var1": ")" << symbol_table.getName(getSymbIDByDerivID(var1)) << R"(")"
-                       << R"(, "shift1": )" << getLagByDerivID(var1)
-                       << R"(, "var2": ")" << symbol_table.getName(getSymbIDByDerivID(var2)) << R"(")"
-                       << R"(, "shift2": )" << getLagByDerivID(var2);
+          d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
 
-      hessian_output << R"(, "val": ")";
-      d2->writeJsonOutput(hessian_output, temp_term_union, tef_terms);
-      hessian_output << R"("})" << endl;
-    }
-  hessian_output << "]}";
-
-  // Writing third derivatives
-  writeJsonTemporaryTerms(temporary_terms_derivatives[3], temp_term_union, third_derivatives_output, tef_terms, "third_derivatives");
-  third_derivatives_output << R"(, "third_derivative": {)"
-                           << R"(  "nrows": )" << equations.size()
-                           << R"(, "ncols": )" << hessianColsNbr * dynJacobianColsNbr
-                           << R"(, "entries": [)";
-  for (auto it = derivatives[3].begin();
-       it != derivatives[3].end(); it++)
-    {
-      if (it != derivatives[3].begin())
-        third_derivatives_output << ", ";
+          if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
+            {
+              int col_idx_sym = getDynJacobianCol(vidx[2]) * dynJacobianColsNbr + getDynJacobianCol(vidx[1]);
+              d_output[i] << ", " << col_idx_sym + 1;
+            }
+          if (i == 3) // Symmetric elements in 3rd derivative
+            {
+              vector<int> idx3{getDynJacobianCol(vidx[1]), getDynJacobianCol(vidx[2]), getDynJacobianCol(vidx[3])};
+              sort(idx3.begin(), idx3.end());
+              do
+                {
+                  int col_idx_sym = (idx3[0]*dynJacobianColsNbr + idx3[1])*dynJacobianColsNbr + idx3[2];
+                  if (col_idx_sym != col_idx)
+                    d_output[i] << ", " << col_idx_sym+1;
+                }
+              while (next_permutation(idx3.begin(), idx3.end()));
+            }
+          if (i > 1)
+            d_output[i] << "]";
 
-      int eq, var1, var2, var3;
-      tie(eq, var1, var2, var3) = vectorToTuple<4>(it->first);
-      expr_t d3 = it->second;
+          if (writeDetails)
+            for (size_t j = 1; j < vidx.size(); j++)
+              d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << symbol_table.getName(getSymbIDByDerivID(vidx[j])) << R"(")"
+                                                                                                                                 << R"(, "shift)" << (i > 1 ? to_string(j) : "") << R"(": )" << getLagByDerivID(vidx[j]);
 
-      if (writeDetails)
-        third_derivatives_output << R"({"eq": )" << eq + 1;
-      else
-        third_derivatives_output << R"({"row": )" << eq + 1;
-
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
-      int id3 = getDynJacobianCol(var3);
-      set<int> cols;
-      cols.insert(id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3);
-      cols.insert(id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2);
-      cols.insert(id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3);
-      cols.insert(id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1);
-      cols.insert(id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2);
-      cols.insert(id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1);
-
-      third_derivatives_output << R"(, "col": [)";
-      for (auto it2 = cols.begin(); it2 != cols.end(); it2++)
-        {
-          if (it2 != cols.begin())
-            third_derivatives_output << ", ";
-          third_derivatives_output << *it2 + 1;
+          d_output[i] << R"(, "val": ")";
+          d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
+          d_output[i] << R"("})" << endl;
         }
-      third_derivatives_output << "]";
+     d_output[i] << "]}";
 
-      if (writeDetails)
-        third_derivatives_output << R"(, "var1": ")" << symbol_table.getName(getSymbIDByDerivID(var1)) << R"(")"
-                                 << R"(, "shift1": )" << getLagByDerivID(var1)
-                                 << R"(, "var2": ")" << symbol_table.getName(getSymbIDByDerivID(var2)) << R"(")"
-                                 << R"(, "shift2": )" << getLagByDerivID(var2)
-                                 << R"(, "var3": ")" << symbol_table.getName(getSymbIDByDerivID(var3)) << R"(")"
-                                 << R"(, "shift3": )" << getLagByDerivID(var3);
-
-      third_derivatives_output << R"(, "val": ")";
-      d3->writeJsonOutput(third_derivatives_output, temp_term_union, tef_terms);
-      third_derivatives_output << R"("})" << endl;
+     ncols *= dynJacobianColsNbr;
     }
-  third_derivatives_output << "]}";
 
   if (writeDetails)
     output << R"("dynamic_model": {)";
   else
     output << R"("dynamic_model_simple": {)";
-  output << model_local_vars_output.str()
-         << ", " << model_output.str()
-         << ", " << jacobian_output.str()
-         << ", " << hessian_output.str()
-         << ", " << third_derivatives_output.str()
-         << "}";
+  output << model_local_vars_output.str();
+  for (const auto &it : d_output)
+    output << ", " << it.str();
+  output << "}";
 }
 
 void
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 63a4580d..9d41d45f 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1482,7 +1482,7 @@ StaticModel::writeStaticModel(const string &basename,
   int JacobianColsNbr = symbol_table.endo_nbr();
   int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
 
-  auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); };;
+  auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); };
 
   // Write Jacobian w.r. to endogenous only
   if (!derivatives[1].empty())
@@ -2808,165 +2808,94 @@ StaticModel::writeJsonOutput(ostream &output) const
 void
 StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) const
 {
-  ostringstream model_local_vars_output;   // Used for storing model local vars
-  ostringstream model_output;              // Used for storing model
-  ostringstream jacobian_output;           // Used for storing jacobian equations
-  ostringstream hessian_output;            // Used for storing Hessian equations
-  ostringstream third_derivatives_output;  // Used for storing third order derivatives equations
+  ostringstream model_local_vars_output;  // Used for storing model local vars
+  vector<ostringstream> d_output(derivatives.size());  // Derivatives output (at all orders, including 0=residual)
 
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temp_term_union;
 
   writeJsonModelLocalVariables(model_local_vars_output, tef_terms);
 
-  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, model_output, tef_terms, "");
-  model_output << ", ";
-  writeJsonModelEquations(model_output, true);
-
-  int nrows = equations.size();
-  int JacobianColsNbr = symbol_table.endo_nbr();
-  int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
-
-  // Write Jacobian w.r. to endogenous only
-  writeJsonTemporaryTerms(temporary_terms_derivatives[1], temp_term_union, jacobian_output, tef_terms, "jacobian");
-  jacobian_output << R"(, "jacobian": {)"
-                  << R"(  "nrows": )" << nrows
-                  << R"(, "ncols": )" << JacobianColsNbr
-                  << R"(, "entries": [)";
-  for (auto it = derivatives[1].begin();
-       it != derivatives[1].end(); it++)
-    {
-      if (it != derivatives[1].begin())
-        jacobian_output << ", ";
-
-      int eq, var;
-      tie(eq, var) = vectorToTuple<2>(it->first);
-      int symb_id = getSymbIDByDerivID(var);
-      int col = symbol_table.getTypeSpecificID(symb_id);
-      expr_t d1 = it->second;
-
-      if (writeDetails)
-        jacobian_output << R"({"eq": )" << eq + 1;
-      else
-        jacobian_output << R"({"row": )" << eq + 1;
-
-      jacobian_output << R"(, "col": )" << col + 1;
+  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
+  d_output[0] << ", ";
+  writeJsonModelEquations(d_output[0], true);
 
-      if (writeDetails)
-        jacobian_output << R"(, "var": ")" << symbol_table.getName(symb_id) << R"(")";
-
-      jacobian_output << R"(, "val": ")";
-      d1->writeJsonOutput(jacobian_output, temp_term_union, tef_terms);
-      jacobian_output << R"("})" << endl;
-    }
-  jacobian_output << "]}";
-
-  // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
-  writeJsonTemporaryTerms(temporary_terms_derivatives[2], temp_term_union, hessian_output, tef_terms, "hessian");
-  hessian_output << R"(, "hessian": {)"
-                 << R"(  "nrows": )" << equations.size()
-                 << R"(, "ncols": )" << hessianColsNbr
-                 << R"(, "entries": [)";
-  for (auto it = derivatives[2].begin();
-       it != derivatives[2].end(); it++)
+  auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); };
+  int ncols = symbol_table.endo_nbr();
+  for (size_t i = 1; i < derivatives.size(); i++)
     {
-      if (it != derivatives[2].begin())
-        hessian_output << ", ";
-
-      int eq, var1, var2;
-      tie(eq, var1, var2) = vectorToTuple<3>(it->first);
-      int symb_id1 = getSymbIDByDerivID(var1);
-      int symb_id2 = getSymbIDByDerivID(var2);
-      expr_t d2 = it->second;
-
-      int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
-      int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
+      string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative";
+      writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
+      temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
+      d_output[i] << R"(, ")" << matrix_name  << R"(": {)"
+                  << R"(  "nrows": )" << equations.size()
+                  << R"(, "ncols": )" << ncols
+                  << R"(, "entries": [)";
 
-      int col = tsid1*symbol_table.endo_nbr()+tsid2;
-      int col_sym = tsid2*symbol_table.endo_nbr()+tsid1;
+      for (auto it = derivatives[i].begin(); it != derivatives[i].end(); ++it)
+        {
+          if (it != derivatives[i].begin())
+            d_output[i] << ", ";
 
-      if (writeDetails)
-        hessian_output << R"({"eq": )" << eq + 1;
-      else
-        hessian_output << R"({"row": )" << eq + 1;
+          const vector<int> &vidx = it->first;
+          expr_t d = it->second;
+          int eq = vidx[0];
 
-      hessian_output << R"(, "col": [)" << col + 1;
+          int col_idx = 0;
+          for (size_t j = 1; j < vidx.size(); j++)
+            {
+              col_idx *= symbol_table.endo_nbr();
+              col_idx += getJacobCol(vidx[j]);
+            }
 
-      if (writeDetails)
-        hessian_output << R"(, "var1": ")" << symbol_table.getName(symb_id1) << R"(")"
-                       << R"(, "var2": ")" << symbol_table.getName(symb_id2) << R"(")";
+          if (writeDetails)
+            d_output[i] << R"({"eq": )" << eq + 1;
+          else
+            d_output[i] << R"({"row": )" << eq + 1;
 
-      if (symb_id1 != symb_id2)
-        hessian_output << ", " <<  col_sym + 1;
-      hessian_output << "]"
-                     << R"(, "val": ")";
-      d2->writeJsonOutput(hessian_output, temp_term_union, tef_terms);
-      hessian_output << R"("})" << endl;
-    }
-  hessian_output << "]}";
+          d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
 
-  // Writing third derivatives
-  writeJsonTemporaryTerms(temporary_terms_derivatives[3], temp_term_union, third_derivatives_output, tef_terms, "third_derivatives");
-  third_derivatives_output << R"(, "third_derivative": {)"
-                           << R"(  "nrows": )" << equations.size()
-                           << R"(, "ncols": )" << hessianColsNbr * JacobianColsNbr
-                           << R"(, "entries": [)";
-  for (auto it = derivatives[3].begin();
-       it != derivatives[3].end(); it++)
-    {
-      if (it != derivatives[3].begin())
-        third_derivatives_output << ", ";
+          if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
+            {
+              int col_idx_sym = getJacobCol(vidx[2]) * symbol_table.endo_nbr() + getJacobCol(vidx[1]);
+              d_output[i] << ", " << col_idx_sym + 1;
+            }
+          if (i == 3) // Symmetric elements in 3rd derivative
+            {
+              vector<int> idx3{getJacobCol(vidx[1]), getJacobCol(vidx[2]), getJacobCol(vidx[3])};
+              sort(idx3.begin(), idx3.end());
+              do
+                {
+                  int col_idx_sym = (idx3[0]*symbol_table.endo_nbr() + idx3[1])*symbol_table.endo_nbr() + idx3[2];
+                  if (col_idx_sym != col_idx)
+                    d_output[i] << ", " << col_idx_sym+1;
+                }
+              while (next_permutation(idx3.begin(), idx3.end()));
+            }
+          if (i > 1)
+            d_output[i] << "]";
 
-      int eq, var1, var2, var3;
-      tie(eq, var1, var2, var3) = vectorToTuple<4>(it->first);
-      expr_t d3 = it->second;
+          if (writeDetails)
+            for (size_t j = 1; j < vidx.size(); j++)
+              d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << symbol_table.getName(getSymbIDByDerivID(vidx[j])) << R"(")";
 
-      if (writeDetails)
-        third_derivatives_output << R"({"eq": )" << eq + 1;
-      else
-        third_derivatives_output << R"({"row": )" << eq + 1;
-
-      int id1 = getSymbIDByDerivID(var1);
-      int id2 = getSymbIDByDerivID(var2);
-      int id3 = getSymbIDByDerivID(var3);
-      set<int> cols;
-      cols.insert(id1 * hessianColsNbr + id2 * JacobianColsNbr + id3);
-      cols.insert(id1 * hessianColsNbr + id3 * JacobianColsNbr + id2);
-      cols.insert(id2 * hessianColsNbr + id1 * JacobianColsNbr + id3);
-      cols.insert(id2 * hessianColsNbr + id3 * JacobianColsNbr + id1);
-      cols.insert(id3 * hessianColsNbr + id1 * JacobianColsNbr + id2);
-      cols.insert(id3 * hessianColsNbr + id2 * JacobianColsNbr + id1);
-
-      third_derivatives_output << R"(, "col": [)";
-      for (auto it2 = cols.begin(); it2 != cols.end(); it2++)
-        {
-          if (it2 != cols.begin())
-            third_derivatives_output << ", ";
-          third_derivatives_output << *it2 + 1;
+          d_output[i] << R"(, "val": ")";
+          d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
+          d_output[i] << R"("})" << endl;
         }
-      third_derivatives_output << "]";
-
-      if (writeDetails)
-        third_derivatives_output << R"(, "var1": ")" << symbol_table.getName(getSymbIDByDerivID(var1)) << R"(")"
-                                 << R"(, "var2": ")" << symbol_table.getName(getSymbIDByDerivID(var2)) << R"(")"
-                                 << R"(, "var3": ")" << symbol_table.getName(getSymbIDByDerivID(var3)) << R"(")";
+     d_output[i] << "]}";
 
-      third_derivatives_output << R"(, "val": ")";
-      d3->writeJsonOutput(third_derivatives_output, temp_term_union, tef_terms);
-      third_derivatives_output << R"("})" << endl;
+     ncols *= symbol_table.endo_nbr();
     }
-  third_derivatives_output << "]}";
 
   if (writeDetails)
     output << R"("static_model": {)";
   else
     output << R"("static_model_simple": {)";
-  output << model_local_vars_output.str()
-         << ", " << model_output.str()
-         << ", " << jacobian_output.str()
-         << ", " << hessian_output.str()
-         << ", " << third_derivatives_output.str()
-         << "}";
+  output << model_local_vars_output.str();
+  for (const auto &it : d_output)
+    output << ", " << it.str();
+  output << "}";
 }
 
 void
-- 
GitLab