From eb4efa59acd82e94385d54f159f882bba2a335f5 Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtanb@gmail.com>
Date: Tue, 16 Mar 2010 10:15:19 +0100
Subject: [PATCH] second derivatives of model equations w.r.t. parameters.

---
 preprocessor/DynamicModel.cc | 36 +++++++++++++++++++++++++++++++++++-
 preprocessor/DynamicModel.hh |  7 +++++++
 2 files changed, 42 insertions(+), 1 deletion(-)

diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 0b4c51380..28812b20e 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -2644,6 +2644,19 @@ DynamicModel::computeParamsDerivatives()
           residuals_params_derivatives[make_pair(eq, param)] = d1;
         }
 
+      for (first_derivatives_type::const_iterator it2 = residuals_params_derivatives.begin();
+           it2 != residuals_params_derivatives.end(); it2++)
+        {
+          int eq = it2->first.first;
+          int param1 = it2->first.second;
+          NodeID d1 = it2->second;
+
+          NodeID d2 = d1->getDerivative(param);
+          if (d2 == Zero)
+            continue;
+          residuals_params_second_derivatives[make_pair(eq, make_pair(param1, param))] = d2;
+        }
+
       for (first_derivatives_type::const_iterator it2 = first_derivatives.begin();
            it2 != first_derivatives.end(); it2++)
         {
@@ -2706,6 +2719,7 @@ void
 DynamicModel::writeParamsDerivativesFile(const string &basename) const
 {
   if (!residuals_params_derivatives.size()
+      && !residuals_params_second_derivatives.size()
       && !jacobian_params_derivatives.size()
       && !jacobian_params_second_derivatives.size()
       && !hessian_params_derivatives.size())
@@ -2720,7 +2734,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  paramsDerivsFile << "function [rp, gp, gpp, hp] = " << basename << "_params_derivs(y, x, params, it_)" << endl
+  paramsDerivsFile << "function [rp, rpp, gp, gpp, hp] = " << basename << "_params_derivs(y, x, params, it_)" << endl
                    << "%" << endl
                    << "% Warning : this file is generated automatically by Dynare" << endl
                    << "%           from model file (.mod)" << endl << endl;
@@ -2745,6 +2759,26 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       paramsDerivsFile << ";" << endl;
     }
 
+  // Write parameter second derivatives
+  paramsDerivsFile << "rpp = zeros(" << equation_number() << ", " << symbol_table.param_nbr() << ", "
+                   << symbol_table.param_nbr() << ");" << endl;
+
+  for (second_derivatives_type::const_iterator it = residuals_params_second_derivatives.begin();
+       it != residuals_params_second_derivatives.end(); it++)
+    {
+      int eq = it->first.first;
+      int param1 = it->first.second.first;
+      int param2 = it->first.second.second;
+      NodeID d2 = it->second;
+
+      int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
+      int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
+
+      paramsDerivsFile << "rpp(" << eq+1 << ", " << param1_col << ", " << param2_col << ") = ";
+      d2->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms);
+      paramsDerivsFile << ";" << endl;
+    }
+
   // Write jacobian derivatives
   paramsDerivsFile << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", "
                    << symbol_table.param_nbr() << ");" << endl;
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index 93fbc2f21..f143355a6 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -64,6 +64,13 @@ private:
   */
   first_derivatives_type residuals_params_derivatives;
 
+  //! Second derivatives of the residuals w.r. to parameters
+  /*! First index is equation number, second and third indeces are parameters.
+    Only non-null derivatives are stored in the map.
+    Parameter indices are those of the getDerivID() method.
+  */
+  second_derivatives_type residuals_params_second_derivatives;
+
   //! Derivatives of the jacobian w.r. to parameters
   /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter.
     Only non-null derivatives are stored in the map.