output mex file for static model (closes #183)

ComputingTasks.cc

@@ -869,7 +869,7 @@ PlannerObjectiveStatement::computingPass()
 void
 PlannerObjectiveStatement::writeOutput(ostream &output, const string &basename) const
 {
-  model_tree->writeStaticFile(basename + "_objective", false, false);
+  model_tree->writeStaticFile(basename + "_objective", false, false, false);
 }
 
 BVARDensityStatement::BVARDensityStatement(int maxnlags_arg, const OptionsList &options_list_arg) :

ExprNode.cc

@@ -615,6 +615,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           i = datatree.getDynJacobianCol(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 oCStaticModel:
         case oMatlabStaticModel:
         case oMatlabStaticModelSparse:
           i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
@@ -669,6 +670,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           else
             output <<  "x[it_" << lag << "+" << i << "*nb_row_x]";
           break;
+        case oCStaticModel:
         case oMatlabStaticModel:
         case oMatlabStaticModelSparse:
           output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -710,6 +712,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           else
             output <<  "x[it_" << lag << "+" << i << "*nb_row_x]";
           break;
+        case oCStaticModel:
         case oMatlabStaticModel:
         case oMatlabStaticModelSparse:
           output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);

ExprNode.hh

@@ -65,6 +65,7 @@ enum ExprNodeOutputType
     oMatlabStaticModelSparse,                     //!< Matlab code, static block decomposed model
     oMatlabDynamicModelSparse,                    //!< Matlab code, dynamic block decomposed model
     oCDynamicModel,                               //!< C code, dynamic model
+    oCStaticModel,                                //!< C code, static model
     oMatlabOutsideModel,                          //!< Matlab code, outside model block (for example in initval)
     oLatexStaticModel,                            //!< LaTeX code, static model
     oLatexDynamicModel,                           //!< LaTeX code, dynamic model
@@ -84,7 +85,7 @@ enum ExprNodeOutputType
                                 || (output_type) == oMatlabDynamicSparseSteadyStateOperator \
                                 || (output_type) == oSteadyStateFile)
 
-#define IS_C(output_type) ((output_type) == oCDynamicModel || (output_type) == oCDynamicSteadyStateOperator)
+#define IS_C(output_type) ((output_type) == oCDynamicModel || (output_type) == oCStaticModel || (output_type) == oCDynamicSteadyStateOperator)
 
 #define IS_LATEX(output_type) ((output_type) == oLatexStaticModel       \
                                || (output_type) == oLatexDynamicModel   \

ModFile.cc

@@ -488,7 +488,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
   unlink(statfile.c_str());
 
   if (!use_dll)
-    mOutputFile << "erase_compiled_function('" + basename + "_dynamic');" << endl;
+    {
+      mOutputFile << "erase_compiled_function('" + basename + "_static');" << endl;
+      mOutputFile << "erase_compiled_function('" + basename + "_dynamic');" << endl;
+    }
 
 #if defined(_WIN32) || defined(__CYGWIN32__)
   // If using USE_DLL with MSVC, check that the user didn't use a function not supported by MSVC (because MSVC doesn't comply with C99 standard)
@@ -524,16 +527,24 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
       // Some mex commands are enclosed in an eval(), because otherwise it will make Octave fail
 #if defined(_WIN32) || defined(__CYGWIN32__)
       if (msvc)
-        mOutputFile << "    eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_dynamic.c')" << endl; // MATLAB/Windows + Microsoft Visual C++
+        // MATLAB/Windows + Microsoft Visual C++
+        mOutputFile << "    eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_dynamic.c')" << endl
+                    << "    eval('mex -O LINKFLAGS=\"$LINKFLAGS /export:Dynamic\" " << basename << "_static.c')" << endl;
       else if (cygwin)
-        mOutputFile << "    eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_dynamic.c')" << endl; // MATLAB/Windows + Cygwin g++
+        // MATLAB/Windows + Cygwin g++
+        mOutputFile << "    eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_dynamic.c')" << endl
+                    << "    eval('mex -O PRELINK_CMDS1=\"echo EXPORTS > mex.def & echo mexFunction >> mex.def & echo Dynamic >> mex.def\" " << basename << "_static.c')" << endl;
       else
         mOutputFile << "    error('When using the USE_DLL option, you must give either ''cygwin'' or ''msvc'' option to the ''dynare'' command')" << endl;
 #else
 # ifdef __linux__
-      mOutputFile << "    eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_dynamic.c')" << endl; // MATLAB/Linux
+      // MATLAB/Linux
+      mOutputFile << "    eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_dynamic.c')" << endl
+                  << "    eval('mex -O LDFLAGS=''-pthread -shared -Wl,--no-undefined'' " << basename << "_static.c')" << endl;
 # else // MacOS
-      mOutputFile << "    eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' " << basename << "_dynamic.c')" << endl; // MATLAB/MacOS
+      // MATLAB/MacOS
+      mOutputFile << "    eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' " << basename << "_dynamic.c')" << endl
+                  << "    eval('mex -O LDFLAGS=''-Wl,-twolevel_namespace -undefined error -arch \\$ARCHS -Wl,-syslibroot,\\$SDKROOT -mmacosx-version-min=\\$MACOSX_DEPLOYMENT_TARGET -bundle'' " << basename << "_static.c')" << endl;
 # endif
 #endif
       mOutputFile << "else" << endl // Octave
@@ -541,6 +552,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
                   << "        sleep(2)" << endl
                   << "    end" << endl
                   << "    mex " << basename << "_dynamic.c" << endl
+                  << "    mex " << basename << "_static.c" << endl
                   << "end" << endl;
     }
 
@@ -595,7 +607,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
   if (dynamic_model.equation_number() > 0)
     {
       if (!no_static)
-        static_model.writeStaticFile(basename, block, byte_code);
+        static_model.writeStaticFile(basename, block, byte_code, use_dll);
 
       dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option);
       dynamic_model.writeParamsDerivativesFile(basename);

StaticModel.cc

@@ -1140,32 +1140,27 @@ StaticModel::writeStaticMFile(const string &func_name) const
          << "% Status : Computes static model for Dynare" << endl
          << "%" << endl
          << "% Warning : this file is generated automatically by Dynare" << endl
-         << "%           from model file (.mod)" << endl
-         << endl
-         << "residual = zeros( " << equations.size() << ", 1);" << endl
-         << endl
-         << "%" << endl
-         << "% Model equations" << endl
-         << "%" << endl
-         << endl;
+         << "%           from model file (.mod)" << endl << endl;
 
-  deriv_node_temp_terms_t tef_terms;
-  writeModelLocalVariables(output, oMatlabStaticModel, tef_terms);
+  writeStaticModel(output, false);
+  output << "end" << endl;
+  output.close();
+}
 
-  writeTemporaryTerms(temporary_terms, output, oMatlabStaticModel, tef_terms);
+void
+StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
+{
+  ostringstream model_output;    // Used for storing model equations
+  ostringstream jacobian_output; // Used for storing jacobian equations
+  ostringstream hessian_output;  // Used for storing Hessian equations
+  ExprNodeOutputType output_type = (use_dll ? oCStaticModel : oMatlabStaticModel);
 
-  writeModelEquations(output, oMatlabStaticModel);
+  deriv_node_temp_terms_t tef_terms;
+  writeModelLocalVariables(model_output, output_type, tef_terms);
 
-  output << "if ~isreal(residual)" << endl
-         << "  residual = real(residual)+imag(residual).^2;" << endl
-         << "end" << endl
-         << "if nargout >= 2," << endl
-         << "  g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl
-         << endl
-         << "%" << endl
-         << "% Jacobian matrix" << endl
-         << "%" << endl
-         << endl;
+  writeTemporaryTerms(temporary_terms, model_output, output_type, tef_terms);
+
+  writeModelEquations(model_output, output_type);
 
   // Write Jacobian w.r. to endogenous only
   for (first_derivatives_t::const_iterator it = first_derivatives.begin();
@@ -1175,72 +1170,212 @@ StaticModel::writeStaticMFile(const string &func_name) const
       int symb_id = it->first.second;
       expr_t d1 = it->second;
 
-      output << "  g1(" << eq+1 << "," << symbol_table.getTypeSpecificID(symb_id)+1 << ")=";
-      d1->writeOutput(output, oMatlabStaticModel, temporary_terms, tef_terms);
-      output << ";" << endl;
+      jacobian_output << "  g1";
+      jacobianHelper(jacobian_output, eq, symbol_table.getTypeSpecificID(symb_id), output_type);
+      jacobian_output << "=";
+      d1->writeOutput(jacobian_output, output_type, temporary_terms, tef_terms);
+      jacobian_output << ";" << endl;
     }
 
-  output << "  if ~isreal(g1)" << endl
-         << "    g1 = real(g1)+imag(g1).^2;" << endl
-         << "  end" << endl
-         << "end" << endl
-         << "if nargout >= 3," << endl
-         << "%" << endl
-         << "% Hessian matrix" << endl
-         << "%" << endl
-         << endl;
-
   int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
-  if (second_derivatives.size())
+  // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
+  int k = 0; // Keep the line of a 2nd derivative in v2
+  for (second_derivatives_t::const_iterator it = second_derivatives.begin();
+       it != second_derivatives.end(); it++)
     {
-      output << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl;
+      int eq = it->first.first;
+      int symb_id1 = it->first.second.first;
+      int symb_id2 = it->first.second.second;
+      expr_t d2 = it->second;
 
-      // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
-      int k = 0; // Keep the line of a 2nd derivative in v2
-      for (second_derivatives_t::const_iterator it = second_derivatives.begin();
-           it != second_derivatives.end(); it++)
-        {
-          int eq = it->first.first;
-          int symb_id1 = it->first.second.first;
-          int symb_id2 = it->first.second.second;
-          expr_t d2 = it->second;
+      int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
+      int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
 
-          int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
-          int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
+      int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
+      int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
 
-          int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
-          int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
+      hessianHelper(hessian_output, k, 0, output_type);
+      hessian_output << "=" << eq + 1 << ";" << endl;
 
-          output << "v2(" << k+1 << ",1)=" << eq + 1 << ";" << endl
-                 << "v2(" << k+1 << ",2)=" << col_nb + 1 << ";" << endl
-                 << "v2(" << k+1 << ",3)=";
-          d2->writeOutput(output, oMatlabStaticModel, temporary_terms, tef_terms);
-          output << ";" << endl;
+      hessianHelper(hessian_output, k, 1, output_type);
+      hessian_output << "=" << col_nb + 1 << ";" << endl;
 
-          k++;
+      hessianHelper(hessian_output, k, 2, output_type);
+      hessian_output << "=";
+      d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
+      hessian_output << ";" << endl;
 
-          // Treating symetric elements
-          if (symb_id1 != symb_id2)
-            {
-              output << "v2(" << k+1 << ",1)=" << eq + 1 << ";" << endl
-                     << "v2(" << k+1 << ",2)=" << col_nb_sym + 1 << ";" << endl
-                     << "v2(" << k+1 << ",3)=v2(" << k << ",3);" << endl;
-              k++;
-            }
+      k++;
+
+      // Treating symetric elements
+      if (symb_id1 != symb_id2)
+        {
+          hessianHelper(hessian_output, k, 0, output_type);
+          hessian_output << "=" << eq + 1 << ";" << endl;
+
+          hessianHelper(hessian_output, k, 1, output_type);
+          hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+
+          hessianHelper(hessian_output, k, 2, output_type);
+          hessian_output << "=";
+          hessianHelper(hessian_output, k-1, 2, output_type);
+          hessian_output << ";" << endl;
+
+          k++;
         }
+    }
 
-      output << "  g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");" << endl;
+  if (!use_dll)
+    {
+      StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl
+                   << "%" << endl
+                   << "% Model equations" << endl
+                   << "%" << endl << endl
+                   << model_output.str()
+                   << "if ~isreal(residual)" << endl
+                   << "  residual = real(residual)+imag(residual).^2;" << endl
+                   << "end" << endl
+                   << "if nargout >= 2," << endl
+                   << "  g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl << endl
+                   << "  %" << endl
+                   << "  % Jacobian matrix" << endl
+                   << "  %" << endl << endl
+                   << jacobian_output.str()
+                   << "  if ~isreal(g1)" << endl
+                   << "    g1 = real(g1)+imag(g1).^2;" << endl
+                   << "  end" << endl
+                   << "end" << endl
+                   << "if nargout >= 3," << endl
+                   << "  %" << endl
+                   << "  % Hessian matrix" << endl
+                   << "  %" << endl
+                   << endl;
+
+      if (second_derivatives.size())
+        StaticOutput << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl
+                     << hessian_output.str()
+                     << "  g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");" << endl;
+      else
+        StaticOutput << "  g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl;
+      StaticOutput << "end" << endl;
     }
-  else // Either hessian is all zero, or we didn't compute it
-    output << "  g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl;
+  else
+    {
+      StaticOutput << "void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2)" << endl
+                   << "{" << endl
+                   << "  double lhs, rhs;" << endl
+                   << endl
+                   << "  /* Residual equations */" << endl
+                   << model_output.str()
+                   << "  /* Jacobian  */" << endl
+                   << "  if (g1 == NULL)" << endl
+                   << "    return;" << endl
+                   << "  else" << endl
+                   << "    {" << endl
+                   << jacobian_output.str()
+                   << "    }" << endl;
+
+      if (second_derivatives.size())
+        StaticOutput << "  /* Hessian for endogenous and exogenous variables */" << endl
+                     << "  if (v2 == NULL)" << endl
+                     << "    return;" << endl
+                     << "  else" << endl
+                     << "    {" << endl
+                     << hessian_output.str()
+                     << "    }" << endl;
+    }
+}
 
-  output << "end" << endl; // Close the if nargout >= 3 statement
-  output << "end" << endl; // Close the *_static function
+void
+StaticModel::writeStaticCFile(const string &func_name) const
+{
+  // Writing comments and function definition command
+  string filename = func_name + "_static.c";
+
+  ofstream output;
+  output.open(filename.c_str(), 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
+         << "#include <math.h>" << endl
+         << "#include \"mex.h\"" << endl
+         << endl
+         << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
+         << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
+
+
+  // Write function definition if oPowerDeriv is used
+  writePowerDerivCHeader(output);
+
+  // Writing the function body
+  writeStaticModel(output, true);
+  output << "}" << endl << endl;
+
+  // Writing the gateway routine
+  output << "/* The gateway routine */" << endl
+         << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
+         << "{" << endl
+         << "  double *y, *x, *params;" << endl
+         << "  double *residual, *g1, *v2;" << endl
+         << "  int nb_row_x;" << endl
+         << endl
+         << "  /* Create a pointer to the input matrix y. */" << endl
+         << "  y = mxGetPr(prhs[0]);" << endl
+         << endl
+         << "  /* Create a pointer to the input matrix x. */" << endl
+         << "  x = mxGetPr(prhs[1]);" << endl
+         << endl
+         << "  /* Create a pointer to the input matrix params. */" << endl
+         << "  params = mxGetPr(prhs[2]);" << endl
+         << endl
+         << "  /* Gets number of rows of matrix x. */" << endl
+         << "  nb_row_x = mxGetM(prhs[1]);" << endl
+         << endl
+         << "  residual = NULL;" << endl
+         << "  if (nlhs >= 1)" << endl
+         << "    {" << endl
+         << "      /* Set the output pointer to the output matrix residual. */" << endl
+         << "      plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
+         << "      /* Create a C pointer to a copy of the output matrix residual. */" << endl
+         << "      residual = mxGetPr(plhs[0]);" << endl
+         << "    }" << endl
+         << endl
+         << "  g1 = NULL;" << endl
+         << "  if (nlhs >= 2)" << endl
+         << "    {" << endl
+         << "      /* Set the output pointer to the output matrix g1. */" << endl
+         << "      plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
+         << "      /* Create a C pointer to a copy of the output matrix g1. */" << endl
+         << "      g1 = mxGetPr(plhs[1]);" << endl
+         << "    }" << endl
+         << endl
+         << "  v2 = NULL;" << endl
+         << "  if (nlhs >= 3)" << endl
+         << "    {" << endl
+         << "      /* Set the output pointer to the output matrix v2. */" << endl
+         << "      plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3
+         << ", mxREAL);" << endl
+         << "      v2 = mxGetPr(plhs[2]);" << endl
+         << "    }" << endl
+         << endl
+         << "  /* Call the C subroutines. */" << endl
+         << "  Static(y, x, nb_row_x, params, residual, g1, v2);" << endl
+         << "}" << endl << endl;
+  writePowerDeriv(output, true);
   output.close();
 }
 
 void
-StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode) const
+StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll) const
 {
   int r;
 
@@ -1267,6 +1402,8 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode)
       chdir("..");
       writeStaticBlockMFSFile(basename);
     }
+  else if(use_dll)
+    writeStaticCFile(basename);
   else
     writeStaticMFile(basename);
   writeAuxVarRecursiveDefinitions(basename);
@@ -1570,7 +1707,7 @@ StaticModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutp
 {
   output << LEFT_ARRAY_SUBSCRIPT(output_type);
   if (IS_MATLAB(output_type))
-    output << eq_nb + 1 << ", " << col_nb + 1;
+    output << eq_nb + 1 << "," << col_nb + 1;
   else
     output << eq_nb + col_nb *equations.size();
   output << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -1579,7 +1716,7 @@ StaticModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutp
 void
 StaticModel::hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
 {
-  output << LEFT_ARRAY_SUBSCRIPT(output_type);
+  output << "v2" << LEFT_ARRAY_SUBSCRIPT(output_type);
   if (IS_MATLAB(output_type))
     output << row_nb + 1 << ", " << col_nb + 1;
   else

StaticModel.hh

@@ -53,13 +53,15 @@ private:
   //! Writes static model file (standard Matlab version)
   void writeStaticMFile(const string &static_basename) const;
 
+  //! Writes static model file (C version)
+  void writeStaticCFile(const string &func_name) const;
+
+  //! Writes the static model equations and its derivatives
+  void writeStaticModel(ostream &StaticOutput, bool use_dll) const;
+
   //! Writes the static function calling the block to solve (Matlab version)
   void writeStaticBlockMFSFile(const string &basename) const;
 
-  //! Writes static model file (C version)
-  /*! \todo add third derivatives handling */
-  void writeStaticCFile(const string &static_basename) const;
-
   //! Writes the Block reordred structure of the model in M output
   void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
 
@@ -182,7 +184,7 @@ public:
                                    int &u_count_int, bool &file_open) const;
 
   //! Writes static model file
-  void writeStaticFile(const string &basename, bool block, bool bytecode) const;
+  void writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll) const;
 
   //! Writes LaTeX file with the equations of the static model
   void writeLatexFile(const string &basename) const;