trunk preprocessor:

* created a distinct expression tree for the static model (thus giving better sharing of sub-expressions and better computation of temporary terms for the static model)
* for that purpose, created StaticModel and DynamicModel classes (ModelTree still persists, but only contains code shared between StaticModel and DynamicModel)
* removed sparse static file (to be later replaced by new algorithm for steady state computation on large models)


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2592 ac1d8469-bf42-47a9-8791-bf33cf982152

CodeInterpreter.hh

@@ -49,21 +49,6 @@ enum BlockType
     SIMULTAN = 3   //<! Simultaneous time unseparable block
   };
 
-/*enum BlockSimulationType
-  {
-  UNKNOWN = -1,                      //!< Unknown simulation type
-  EVALUATE_FORWARD = 0,             //!< Simple evaluation, normalized variable on left-hand side, forward
-  EVALUATE_BACKWARD = 1,             //!< Simple evaluation, normalized variable on left-hand side, backward
-  SOLVE_FORWARD_SIMPLE = 2,         //!< Block of one equation, newton solver needed, forward
-  SOLVE_BACKWARD_SIMPLE = 3,         //!< Block of one equation, newton solver needed, backward
-  SOLVE_TWO_BOUNDARIES_SIMPLE = 4,   //!< Block of one equation, newton solver needed, forward & ackward
-  SOLVE_FORWARD_COMPLETE = 5,       //!< Block of several equations, newton solver needed, forward
-  SOLVE_BACKWARD_COMPLETE = 6,       //!< Block of several equations, newton solver needed, backward
-  SOLVE_TWO_BOUNDARIES_COMPLETE = 7, //!< Block of several equations, newton solver needed, forward and backwar
-  EVALUATE_FORWARD_R = 8,           //!< Simple evaluation, normalized variable on right-hand side, forward
-  EVALUATE_BACKWARD_R = 9            //!< Simple evaluation, normalized variable on right-hand side, backward
-  };
-*/
 enum BlockSimulationType
   {
     UNKNOWN,                      //!< Unknown simulation type
@@ -131,4 +116,9 @@ enum BinaryOpcode
     oDifferent
   };
 
+enum TrinaryOpcode
+  {
+    oNormcdf
+  };
+
 #endif

ComputingTasks.cc

@@ -38,19 +38,6 @@ SteadyStatement::writeOutput(ostream &output, const string &basename) const
   output << "steady;\n";
 }
 
-SteadySparseStatement::SteadySparseStatement(const OptionsList &options_list_arg) :
-  options_list(options_list_arg)
-{
-}
-
-void
-SteadySparseStatement::writeOutput(ostream &output, const string &basename) const
-{
-  options_list.writeOutput(output);
-  //output << basename << "_static;\n";
-  output << "steady;\n";
-}
-
 CheckStatement::CheckStatement(const OptionsList &options_list_arg) :
   options_list(options_list_arg)
 {
@@ -914,7 +901,7 @@ ModelComparisonStatement::writeOutput(ostream &output, const string &basename) c
   output << "model_comparison(ModelNames_,ModelPriors_,oo_,options_,M_.fname);" << endl;
 }
 
-PlannerObjectiveStatement::PlannerObjectiveStatement(ModelTree *model_tree_arg) :
+PlannerObjectiveStatement::PlannerObjectiveStatement(StaticModel *model_tree_arg) :
   model_tree(model_tree_arg)
 {
 }
@@ -938,7 +925,7 @@ void
 PlannerObjectiveStatement::computingPass()
 {
   model_tree->computeStaticHessian = true;
-  model_tree->computingPass(eval_context_type(), false);
+  model_tree->computingPass(false);
 }
 
 void

ComputingTasks.hh

@@ -25,7 +25,7 @@
 #include "SymbolList.hh"
 #include "SymbolTable.hh"
 #include "Statement.hh"
-#include "ModelTree.hh"
+#include "StaticModel.hh"
 
 class SteadyStatement : public Statement
 {
@@ -36,15 +36,6 @@ public:
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
-class SteadySparseStatement : public Statement
-{
-private:
-  const OptionsList options_list;
-public:
-  SteadySparseStatement(const OptionsList &options_list_arg);
-  virtual void writeOutput(ostream &output, const string &basename) const;
-};
-
 class CheckStatement : public Statement
 {
 private:
@@ -409,12 +400,12 @@ public:
 class PlannerObjectiveStatement : public Statement
 {
 private:
-  ModelTree *model_tree;
+  StaticModel *model_tree;
 public:
   //! Constructor
   /*! \param model_tree_arg the model tree used to store the objective function.
     It is owned by the PlannerObjectiveStatement, and will be deleted by its destructor */
-  PlannerObjectiveStatement(ModelTree *model_tree_arg);
+  PlannerObjectiveStatement(StaticModel *model_tree_arg);
   virtual ~PlannerObjectiveStatement();
   /*! \todo check there are only endogenous variables at the current period in the objective
     (no exogenous, no lead/lag) */

DataTree.cc

@@ -75,8 +75,8 @@ DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
     {
       // Simplify x+(-y) in x-y
       UnaryOpNode *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
-      if (uarg2 != NULL && uarg2->op_code == oUminus)
+      if (uarg2 != NULL && uarg2->get_op_code() == oUminus)
-        return AddMinus(iArg1, uarg2->arg);
+        return AddMinus(iArg1, uarg2->get_arg());
 
       // To treat commutativity of "+"
       // Nodes iArg1 and iArg2 are sorted by index
@@ -116,8 +116,8 @@ DataTree::AddUMinus(NodeID iArg1)
     {
       // Simplify -(-x) in x
       UnaryOpNode *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
-      if (uarg != NULL && uarg->op_code == oUminus)
+      if (uarg != NULL && uarg->get_op_code() == oUminus)
-        return uarg->arg;
+        return uarg->get_arg();
 
       return AddUnaryOp(oUminus, iArg1);
     }

DynamicModel.hh

@@ -0,0 +1,96 @@
+/*
+ * Copyright (C) 2003-2009 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef _DYNAMICMODEL_HH
+#define _DYNAMICMODEL_HH
+
+using namespace std;
+
+#include "StaticModel.hh"
+#include "BlockTriangular.hh"
+
+//! Stores a dynamic model
+class DynamicModel : public ModelTree
+{
+private:
+  //! Writes dynamic model file (Matlab version)
+  void writeDynamicMFile(const string &dynamic_basename) const;
+  //! Writes dynamic model file (C version)
+  /*! \todo add third derivatives handling */
+  void writeDynamicCFile(const string &dynamic_basename) const;
+  //! Writes dynamic model file when SparseDLL option is on
+  void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const;
+  //! Writes the dynamic model equations and its derivatives
+  /*! \todo add third derivatives handling in C output */
+  void writeDynamicModel(ostream &DynamicOutput) const;
+  //! Writes the Block reordred structure of the model in M output
+  void writeModelEquationsOrdered_M(Model_Block *ModelBlock, const string &dynamic_basename) const;
+  //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
+  void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const;
+  //! Computes jacobian and prepares for equation normalization
+  /*! Using values from initval/endval blocks and parameter initializations:
+    - computes the jacobian for the model w.r. to contemporaneous variables
+    - removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
+  */
+  void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m);
+  void BlockLinear(Model_Block *ModelBlock);
+  string reform(string name) const;
+  map_idx_type map_idx;
+  //! Build The incidence matrix form the modeltree
+  void BuildIncidenceMatrix();
+
+  void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock);
+  //! Write derivative code of an equation w.r. to a variable
+  void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const;
+
+public:
+  DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
+  //! Absolute value under which a number is considered to be zero
+  double cutoff;
+  //! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc)
+  double markowitz;
+  //! the file containing the model and the derivatives code
+  ofstream code_file;
+  //! Whether dynamic Jacobian (w.r. to endogenous) should be written
+  bool computeJacobian;
+  //! Whether dynamic Jacobian (w.r. to endogenous and exogenous) should be written
+  bool computeJacobianExo;
+  //! Whether dynamic Hessian (w.r. to endogenous and exogenous) should be written
+  bool computeHessian;
+  //! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written
+  bool computeThirdDerivatives;
+  //! Execute computations (variable sorting + derivation)
+  /*! You must set computeJacobian, computeJacobianExo, computeHessian and computeThirdDerivatives to correct values before calling this function
+    \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */
+  void computingPass(const eval_context_type &eval_context, bool no_tmp_terms);
+  //! Writes model initialization and lead/lag incidence matrix to output
+  void writeOutput(ostream &output) const;
+  //! Complete set to block decompose the model
+  BlockTriangular block_triangular;
+  //! Adds informations for simulation in a binary file
+  void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename,
+                             const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
+  //! Writes dynamic model file
+  void writeDynamicFile(const string &basename) const;
+  //! Converts to static model (only the equations)
+  /*! It assumes that the static model given in argument has just been allocated */
+  void toStatic(StaticModel &static_model) const;
+};
+
+#endif

ExprNode.cc

@@ -99,6 +99,7 @@ ExprNode::writeOutput(ostream &output)
   writeOutput(output, oMatlabOutsideModel, temporary_terms_type());
 }
 
+
 NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
   ExprNode(datatree_arg),
   id(id_arg)
@@ -115,7 +116,6 @@ NumConstNode::computeDerivative(int varID)
   return datatree.Zero;
 }
 
-
 void
 NumConstNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
 {
@@ -158,7 +158,6 @@ NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou
   CompileCode.write(reinterpret_cast<char *>(&vard),sizeof(vard));
 }
 
-
 void
 NumConstNode::collectEndogenous(set<pair<int, int> > &result) const
 {
@@ -169,6 +168,12 @@ NumConstNode::collectExogenous(set<pair<int, int> > &result) const
 {
 }
 
+NodeID
+NumConstNode::toStatic(DataTree &static_datatree) const
+{
+  return static_datatree.AddNumConstant(datatree.num_constants.get(id));
+}
+
 
 VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg) :
   ExprNode(datatree_arg),
@@ -479,7 +484,6 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou
     }
 }
 
-
 void
 VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
                                     temporary_terms_type &temporary_terms,
@@ -493,9 +497,6 @@ VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
     datatree.local_variables_table[symb_id]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
 }
 
-
-
-
 void
 VariableNode::collectEndogenous(set<pair<int, int> > &result) const
 {
@@ -514,6 +515,12 @@ VariableNode::collectExogenous(set<pair<int, int> > &result) const
     datatree.local_variables_table[symb_id]->collectExogenous(result);
 }
 
+NodeID
+VariableNode::toStatic(DataTree &static_datatree) const
+{
+  return static_datatree.AddVariable(datatree.symbol_table.getName(symb_id));
+}
+
 
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) :
   ExprNode(datatree_arg),
@@ -923,6 +930,49 @@ UnaryOpNode::collectExogenous(set<pair<int, int> > &result) const
   arg->collectExogenous(result);
 }
 
+NodeID
+UnaryOpNode::toStatic(DataTree &static_datatree) const
+{
+  NodeID sarg = arg->toStatic(static_datatree);
+  switch(op_code)
+    {
+    case oUminus:
+      return static_datatree.AddUMinus(sarg);
+    case oExp:
+      return static_datatree.AddExp(sarg);
+    case oLog:
+      return static_datatree.AddLog(sarg);
+    case oLog10:
+      return static_datatree.AddLog10(sarg);
+    case oCos:
+      return static_datatree.AddCos(sarg);
+    case oSin:
+      return static_datatree.AddSin(sarg);
+    case oTan:
+      return static_datatree.AddTan(sarg);
+    case oAcos:
+      return static_datatree.AddACos(sarg);
+    case oAsin:
+      return static_datatree.AddASin(sarg);
+    case oAtan:
+      return static_datatree.AddATan(sarg);
+    case oCosh:
+      return static_datatree.AddCosH(sarg);
+    case oSinh:
+      return static_datatree.AddSinH(sarg);
+    case oTanh:
+      return static_datatree.AddTanH(sarg);
+    case oAcosh:
+      return static_datatree.AddACosH(sarg);
+    case oAsinh:
+      return static_datatree.AddASinH(sarg);
+    case oAtanh:
+      return static_datatree.AddATanH(sarg);
+    case oSqrt:
+      return static_datatree.AddSqRt(sarg);
+    }
+}
+
 
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                            BinaryOpcode op_code_arg, const NodeID arg2_arg) :
@@ -1405,6 +1455,45 @@ BinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
   arg2->collectExogenous(result);
 }
 
+NodeID
+BinaryOpNode::toStatic(DataTree &static_datatree) const
+{
+  NodeID sarg1 = arg1->toStatic(static_datatree);
+  NodeID sarg2 = arg2->toStatic(static_datatree);
+  switch(op_code)
+    {
+    case oPlus:
+      return static_datatree.AddPlus(sarg1, sarg2);
+    case oMinus:
+      return static_datatree.AddMinus(sarg1, sarg2);
+    case oTimes:
+      return static_datatree.AddTimes(sarg1, sarg2);
+    case oDivide:
+      return static_datatree.AddDivide(sarg1, sarg2);
+    case oPower:
+      return static_datatree.AddPower(sarg1, sarg2);
+    case oEqual:
+      return static_datatree.AddEqual(sarg1, sarg2);
+    case oMax:
+      return static_datatree.AddMaX(sarg1, sarg2);
+    case oMin:
+      return static_datatree.AddMin(sarg1, sarg2);
+    case oLess:
+      return static_datatree.AddLess(sarg1, sarg2);
+    case oGreater:
+      return static_datatree.AddGreater(sarg1, sarg2);
+    case oLessEqual:
+      return static_datatree.AddLessEqual(sarg1, sarg2);
+    case oGreaterEqual:
+      return static_datatree.AddGreaterEqual(sarg1, sarg2);
+    case oEqualEqual:
+      return static_datatree.AddEqualEqual(sarg1, sarg2);
+    case oDifferent:
+      return static_datatree.AddDifferent(sarg1, sarg2);
+    }
+}
+
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                              TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) :
   ExprNode(datatree_arg),
@@ -1697,6 +1786,19 @@ TrinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
   arg3->collectExogenous(result);
 }
 
+NodeID
+TrinaryOpNode::toStatic(DataTree &static_datatree) const
+{
+  NodeID sarg1 = arg1->toStatic(static_datatree);
+  NodeID sarg2 = arg2->toStatic(static_datatree);
+  NodeID sarg3 = arg2->toStatic(static_datatree);
+  switch(op_code)
+    {
+    case oNormcdf:
+      return static_datatree.AddNormcdf(sarg1, sarg2, sarg3);
+    }
+}
+
 
 UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
                                          int symb_id_arg,
@@ -1792,3 +1894,13 @@ UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutput
   cerr << "UnknownFunctionNode::compile: operation impossible!" << endl;
   exit(EXIT_FAILURE);
 }
+
+NodeID
+UnknownFunctionNode::toStatic(DataTree &static_datatree) const
+{
+  vector<NodeID> static_arguments;
+  for(vector<NodeID>::const_iterator it = arguments.begin();
+      it != arguments.end(); it++)
+    static_arguments.push_back((*it)->toStatic(static_datatree));
+  return static_datatree.AddUnknownFunction(datatree.symbol_table.getName(symb_id), static_arguments);
+}

ExprNode.hh

@@ -86,7 +86,7 @@ typedef map<int, double> eval_context_type;
 class ExprNode
 {
   friend class DataTree;
-  friend class ModelTree;
+  friend class DynamicModel;
   friend class ExprNodeLess;
   friend class NumConstNode;
   friend class VariableNode;
@@ -159,6 +159,12 @@ public:
 
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException) = 0;
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const = 0;
+  //! Creates a static version of this node
+  /*!
+    This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
+    adds the result in the static_datatree argument (and not in the original datatree), and returns it.
+  */
+  virtual NodeID toStatic(DataTree &static_datatree) const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -186,6 +192,7 @@ public:
   virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+  virtual NodeID toStatic(DataTree &static_datatree) const;
 };
 
 //! Symbol or variable node
@@ -214,12 +221,12 @@ public:
   virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+  virtual NodeID toStatic(DataTree &static_datatree) const;
 };
 
 //! Unary operator node
 class UnaryOpNode : public ExprNode
 {
-  friend class DataTree;
 private:
   const NodeID arg;
   const UnaryOpcode op_code;
@@ -243,12 +250,16 @@ public:
   static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+  //! Returns operand
+  NodeID get_arg() const { return(arg); };
+  //! Returns op code
+  UnaryOpcode get_op_code() const { return(op_code); };
+  virtual NodeID toStatic(DataTree &static_datatree) const;
 };
 
 //! Binary operator node
 class BinaryOpNode : public ExprNode
 {
-  friend class ModelTree;
 private:
   const NodeID arg1, arg2;
   const BinaryOpcode op_code;
@@ -273,15 +284,15 @@ public:
   static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
-  virtual NodeID get_arg1() { return(arg1);};
+  //! Returns first operand
-  virtual NodeID get_arg2() { return(arg2);};
+  NodeID get_arg1() const { return(arg1); };
+  //! Returns second operand
+  NodeID get_arg2() const { return(arg2); };
+  //! Returns op code
+  BinaryOpcode get_op_code() const { return(op_code); };
+  virtual NodeID toStatic(DataTree &static_datatree) const;
 };
 
-enum TrinaryOpcode
-  {
-    oNormcdf
-  };
-
 //! Trinary operator node
 class TrinaryOpNode : public ExprNode
 {
@@ -310,6 +321,7 @@ public:
   static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+  virtual NodeID toStatic(DataTree &static_datatree) const;
 };
 
 //! Unknown function node
@@ -337,6 +349,7 @@ public:
   virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+  virtual NodeID toStatic(DataTree &static_datatree) const;
 };
 
 //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model

Makefile

@@ -38,6 +38,8 @@ MAIN_OBJS = \
 	DynareBison.o \
 	ComputingTasks.o \
 	ModelTree.o \
+	StaticModel.o \
+	DynamicModel.o \
 	NumericalConstants.o \
 	NumericalInitialization.o \
 	Shocks.o \

ModFile.cc

@@ -24,7 +24,8 @@
 #include "ModFile.hh"
 
 ModFile::ModFile() : expressions_tree(symbol_table, num_constants),
-                     model_tree(symbol_table, num_constants),
+                     static_model(symbol_table, num_constants),
+                     dynamic_model(symbol_table, num_constants),
                      linear(false)
 {
 }
@@ -58,7 +59,7 @@ ModFile::evalAllExpressions()
     }
 
   // Evaluate model local variables
-  model_tree.fillEvalContext(global_eval_context);
+  dynamic_model.fillEvalContext(global_eval_context);
 
   cout << "done" << endl;
 
@@ -100,7 +101,7 @@ ModFile::checkPass()
     || mod_file_struct.ramsey_policy_present;
 
   // Allow empty model only when doing a standalone BVAR estimation
-  if (model_tree.equation_number() == 0
+  if (dynamic_model.equation_number() == 0
       && (mod_file_struct.check_present
           || mod_file_struct.simul_present
           || stochastic_statement_present))
@@ -109,7 +110,7 @@ ModFile::checkPass()
       exit(EXIT_FAILURE);
     }
 
-  if (mod_file_struct.simul_present && stochastic_statement_present && model_tree.mode==0)
+  if (mod_file_struct.simul_present && stochastic_statement_present && dynamic_model.mode == 0)
     {
       cerr << "ERROR: A .mod file cannot contain both a simul command and one of {stoch_simul, estimation, forecast, osr, ramsey_policy}" << endl;
       exit(EXIT_FAILURE);
@@ -125,23 +126,29 @@ ModFile::checkPass()
   */
   if (!mod_file_struct.ramsey_policy_present
       && !((mod_file_struct.bvar_density_present || mod_file_struct.bvar_forecast_present)
-           && model_tree.equation_number() == 0)
+           && dynamic_model.equation_number() == 0)
-      && (model_tree.equation_number() != symbol_table.endo_nbr()))
+      && (dynamic_model.equation_number() != symbol_table.endo_nbr()))
     {
-      cerr << "ERROR: There are " << model_tree.equation_number() << " equations but " << symbol_table.endo_nbr() << " endogenous variables!" << endl;
+      cerr << "ERROR: There are " << dynamic_model.equation_number() << " equations but " << symbol_table.endo_nbr() << " endogenous variables!" << endl;
       exit(EXIT_FAILURE);
     }
+
+  cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
 }
 
 void
 ModFile::computingPass(bool no_tmp_terms)
 {
   // Mod file may have no equation (for example in a standalone BVAR estimation)
-  if (model_tree.equation_number() > 0)
+  if (dynamic_model.equation_number() > 0)
     {
-      // Set things to compute
+      // Compute static model and its derivatives
+      dynamic_model.toStatic(static_model);
+      static_model.computingPass(no_tmp_terms);
+
+      // Set things to compute for dynamic model
       if (mod_file_struct.simul_present)
-        model_tree.computeJacobian = true;
+        dynamic_model.computeJacobian = true;
       else
         {
           if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3)
@@ -149,13 +156,13 @@ ModFile::computingPass(bool no_tmp_terms)
               cerr << "Incorrect order option..." << endl;
               exit(EXIT_FAILURE);
             }
-          model_tree.computeJacobianExo = true;
+          dynamic_model.computeJacobianExo = true;
           if (mod_file_struct.order_option >= 2)
-            model_tree.computeHessian = true;
+            dynamic_model.computeHessian = true;
           if (mod_file_struct.order_option == 3)
-            model_tree.computeThirdDerivatives = true;
+            dynamic_model.computeThirdDerivatives = true;
         }
-      model_tree.computingPass(global_eval_context, no_tmp_terms);
+      dynamic_model.computingPass(global_eval_context, no_tmp_terms);
     }
   for(vector<Statement *>::iterator it = statements.begin();
       it != statements.end(); it++)
@@ -209,18 +216,17 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
               << "warning warning_old_state" << endl;
   mOutputFile << "logname_ = '" << basename << ".log';" << endl;
   mOutputFile << "diary " << basename << ".log" << endl;
-  mOutputFile << "options_.model_mode = " << model_tree.mode << ";\n";
+  mOutputFile << "options_.model_mode = " << dynamic_model.mode << ";\n";
-  if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode)
+  if (dynamic_model.mode == eSparseMode || dynamic_model.mode == eSparseDLLMode)
     {
       mOutputFile << "addpath " << basename << ";\n";
-      mOutputFile << "delete('" << basename << "_static.m');\n";
       mOutputFile << "delete('" << basename << "_dynamic.m');\n";
     }
 
 
-  if (model_tree.equation_number() > 0)
+  if (dynamic_model.equation_number() > 0)
     {
-      if (model_tree.mode == eDLLMode || model_tree.mode == eSparseDLLMode)
+      if (dynamic_model.mode == eDLLMode || dynamic_model.mode == eSparseDLLMode)
         {
           mOutputFile << "if exist('" << basename << "_static.c')" << endl;
           mOutputFile << "   clear " << basename << "_static" << endl;
@@ -244,11 +250,11 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
   if (linear == 1)
     mOutputFile << "options_.linear = 1;" << endl;
 
-  if (model_tree.equation_number() > 0)
+  if (dynamic_model.equation_number() > 0)
     {
-      model_tree.writeOutput(mOutputFile);
+      dynamic_model.writeOutput(mOutputFile);
-      model_tree.writeStaticFile(basename);
+      static_model.writeStaticFile(basename);
-      model_tree.writeDynamicFile(basename);
+      dynamic_model.writeDynamicFile(basename);
     }
 
   // Print statements
@@ -259,7 +265,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
   mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl;
   mOutputFile << "diary off" << endl;
 
-  if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode)
+  if (dynamic_model.mode == eSparseMode || dynamic_model.mode == eSparseDLLMode)
     mOutputFile << "rmpath " << basename << ";\n";
 
   mOutputFile << endl << "disp(['Total computing time : ' dynsec2hms(toc) ]);" << endl;

ModFile.hh

@@ -28,7 +28,8 @@ using namespace std;
 #include "SymbolTable.hh"
 #include "NumericalConstants.hh"
 #include "NumericalInitialization.hh"
-#include "ModelTree.hh"
+#include "StaticModel.hh"
+#include "DynamicModel.hh"
 #include "VariableTable.hh"
 #include "Statement.hh"
 
@@ -44,8 +45,10 @@ public:
   NumericalConstants num_constants;
   //! Expressions outside model block
   DataTree expressions_tree;
-  //! Model equations and their derivatives
+  //! Static model
-  ModelTree model_tree;
+  StaticModel static_model;
+  //! Dynamic model
+  DynamicModel dynamic_model;
   //! Option linear
   bool linear;
   //! Global evaluation context
@@ -69,7 +72,7 @@ public:
   //! Execute computations
   /*! \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
   void computingPass(bool no_tmp_terms);
-  //! Writes Matlab/Scilab output files
+  //! Writes Matlab/Octave output files
   /*!
     \param basename The base name used for writing output files. Should be the name of the mod file without its extension
     \param clear_all Should a "clear all" instruction be written to output ?

ModelTree.hh

@@ -26,12 +26,8 @@ using namespace std;
 #include <vector>
 #include <map>
 #include <ostream>
-#include <algorithm>
 
-#include "SymbolTable.hh"
-#include "NumericalConstants.hh"
 #include "DataTree.hh"
-#include "BlockTriangular.hh"
 
 //! The three in which ModelTree can work
 enum ModelTreeMode
@@ -42,10 +38,10 @@ enum ModelTreeMode
     eSparseDLLMode //!< Sparse DLL mode (static file in Matlab, dynamic file in C with block decomposition plus a binary file)
   };
 
-//! Stores a model's equations and derivatives
+//! Shared code for static and dynamic models
 class ModelTree : public DataTree
 {
-private:
+protected:
   //! Stores declared equations
   vector<BinaryOpNode *> equations;
 
@@ -77,19 +73,13 @@ private:
 
   //! Temporary terms (those which will be noted Txxxx)
   temporary_terms_type temporary_terms;
-  map_idx_type map_idx;
 
   //! Computes derivatives of ModelTree
   void derive(int order);
   //! Write derivative of an equation w.r. to a variable
   void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
-  //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const;
   //! Computes temporary terms
   void computeTemporaryTerms(int order);
-  void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock);
-  //! Build The incidence matrix form the modeltree
-  void BuildIncidenceMatrix();
   //! Writes temporary terms
   void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const;
   //! Writes model local variables
@@ -97,39 +87,6 @@ private:
   void writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const;
   //! Writes model equations
   void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const;
-  //! Writes the static model equations and its derivatives
-  /*! \todo handle hessian in C output */
-  void writeStaticModel(ostream &StaticOutput) const;
-  //! Writes the dynamic model equations and its derivatives
-  /*! \todo add third derivatives handling in C output */
-  void writeDynamicModel(ostream &DynamicOutput) const;
-  //! Writes the Block reordred structure of the model in M output
-  void writeModelEquationsOrdered_M(Model_Block *ModelBlock, const string &dynamic_basename) const;
-  //! Writes the Block reordred structure of the static model in M output
-  void writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const string &static_basename) const;
-  //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
-  void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const;
-  //! Writes static model file (Matlab version)
-  void writeStaticMFile(const string &static_basename) const;
-  //! Writes static model file (C version)
-  void writeStaticCFile(const string &static_basename) const;
-  //! Writes static model file when Sparse option is on (Matlab version)
-  void writeSparseStaticMFile(const string &static_basename, const string &basename, const int mode) const;
-  //! Writes dynamic model file (Matlab version)
-  void writeDynamicMFile(const string &dynamic_basename) const;
-  //! Writes dynamic model file (C version)
-  /*! \todo add third derivatives handling */
-  void writeDynamicCFile(const string &dynamic_basename) const;
-  //! Writes dynamic model file when SparseDLL option is on
-  void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const;
-  //! Computes jacobian and prepares for equation normalization
-  /*! Using values from initval/endval blocks and parameter initializations:
-    - computes the jacobian for the model w.r. to contemporaneous variables
-    - removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
-  */
-  void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m);
-  void BlockLinear(Model_Block *ModelBlock);
-  string reform(string name) const;
 
   //! Writes either (i+1,j+1) or [i+j*n_i] whether we are in Matlab or C mode
   void matrixHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
@@ -138,41 +95,8 @@ public:
   ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
   //! Mode in which the ModelTree is supposed to work (Matlab, DLL or SparseDLL)
   ModelTreeMode mode;
-  //! Absolute value under which a number is considered to be zero
-  double cutoff;
-  //! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc)
-  double markowitz;
-  //! Use a graphical and symbolic version of the symbolic gaussian elimination (new_SGE = false) or use direct gaussian elimination (new_SGE = true)
-  bool new_SGE;
-  //! the file containing the model and the derivatives code
-  ofstream code_file;
   //! Declare a node as an equation of the model
   void addEquation(NodeID eq);
-  //! Whether dynamic Jacobian (w.r. to endogenous) should be written
-  bool computeJacobian;
-  //! Whether dynamic Jacobian (w.r. to endogenous and exogenous) should be written
-  bool computeJacobianExo;
-  //! Whether dynamic Hessian (w.r. to endogenous and exogenous) should be written
-  bool computeHessian;
-  //! Whether static Hessian (w.r. to endogenous only) should be written
-  bool computeStaticHessian;
-  //! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written
-  bool computeThirdDerivatives;
-  //! Execute computations (variable sorting + derivation)
-  /*! You must set computeJacobian, computeJacobianExo, computeHessian, computeStaticHessian and computeThirdDerivatives to correct values before calling this function
-    \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
-  void computingPass(const eval_context_type &eval_context, bool no_tmp_terms);
-  //! Writes model initialization and lead/lag incidence matrix to output
-  void writeOutput(ostream &output) const;
-  //! Writes static model file
-  void writeStaticFile(const string &basename) const;
-  //! Writes dynamic model file
-  void writeDynamicFile(const string &basename) const;
-  //! Complete set to block decompose the model
-  BlockTriangular block_triangular;
-  //! Adds informations for simulation in a binary file
-  void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename,
-                             const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
   //! Returns the number of equations in the model
   int equation_number() const;
 };

ParsingDriver.cc

@@ -392,7 +392,7 @@ ParsingDriver::end_homotopy()
 void
 ParsingDriver::begin_model()
 {
-  set_current_data_tree(&mod_file->model_tree);
+  set_current_data_tree(&mod_file->dynamic_model);
 }
 
 void
@@ -585,9 +585,6 @@ ParsingDriver::add_to_row(NodeID v)
 void
 ParsingDriver::steady()
 {
-  if (mod_file->model_tree.mode == eSparseMode || mod_file->model_tree.mode == eSparseDLLMode)
-    mod_file->addStatement(new SteadySparseStatement(options_list));
-  else
   mod_file->addStatement(new SteadyStatement(options_list));
   options_list.clear();
 }
@@ -619,10 +616,10 @@ ParsingDriver::option_num(const string &name_option, const string &opt)
       && (options_list.num_options.find(name_option) != options_list.num_options.end()))
     error("option " + name_option + " declared twice");
 
-  if ((name_option == "periods") && (mod_file->model_tree.mode == eSparseDLLMode || mod_file->model_tree.mode == eSparseMode))
+  if ((name_option == "periods") && (mod_file->dynamic_model.mode == eSparseDLLMode || mod_file->dynamic_model.mode == eSparseMode))
-    mod_file->model_tree.block_triangular.periods = atoi(opt.c_str());
+    mod_file->dynamic_model.block_triangular.periods = atoi(opt.c_str());
   else if (name_option == "cutoff")
-    mod_file->model_tree.cutoff = atof(opt.c_str());
+    mod_file->dynamic_model.cutoff = atof(opt.c_str());
 
   options_list.num_options[name_option] = opt;
 }
@@ -686,7 +683,7 @@ void ParsingDriver::stoch_simul()
 
 void ParsingDriver::simulate()
 {
-  if (mod_file->model_tree.mode == eSparseDLLMode || mod_file->model_tree.mode == eSparseMode)
+  if (mod_file->dynamic_model.mode == eSparseDLLMode || mod_file->dynamic_model.mode == eSparseMode)
     simul_sparse();
   else
     simul();
@@ -695,7 +692,7 @@ void ParsingDriver::simulate()
 void
 ParsingDriver::simul_sparse()
 {
-  mod_file->addStatement(new SimulSparseStatement(options_list, mod_file->model_tree.mode));
+  mod_file->addStatement(new SimulSparseStatement(options_list, mod_file->dynamic_model.mode));
   options_list.clear();
 }
 
@@ -1033,7 +1030,7 @@ ParsingDriver::run_model_comparison()
 void
 ParsingDriver::begin_planner_objective()
 {
-  set_current_data_tree(new ModelTree(mod_file->symbol_table, mod_file->num_constants));
+  set_current_data_tree(new StaticModel(mod_file->symbol_table, mod_file->num_constants));
 }
 
 void
@@ -1043,7 +1040,7 @@ ParsingDriver::end_planner_objective(NodeID expr)
   NodeID eq = model_tree->AddEqual(expr, model_tree->Zero);
   model_tree->addEquation(eq);
 
-  mod_file->addStatement(new PlannerObjectiveStatement(model_tree));
+  mod_file->addStatement(new PlannerObjectiveStatement(dynamic_cast<StaticModel *>(model_tree)));
 
   reset_data_tree();
 }
@@ -1120,7 +1117,7 @@ ParsingDriver::change_type(SymbolType new_type, vector<string *> *var_list)
 
       // Check if symbol already used in a VariableNode
       if (mod_file->expressions_tree.isSymbolUsed(id)
-          || mod_file->model_tree.isSymbolUsed(id))
+          || mod_file->dynamic_model.isSymbolUsed(id))
         error("You cannot modify the type of symbol " + **it + " after having used it in an expression");
 
       mod_file->symbol_table.changeType(id, new_type);

StaticModel.cc

@@ -0,0 +1,289 @@
+/*
+ * Copyright (C) 2003-2009 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "StaticModel.hh"
+
+StaticModel::StaticModel(SymbolTable &symbol_table_arg,
+                         NumericalConstants &num_constants_arg) :
+  ModelTree(symbol_table_arg, num_constants_arg),
+  computeStaticHessian(false)
+{
+}
+
+void
+StaticModel::writeStaticMFile(const string &static_basename) const
+{
+  string filename = static_basename + ".m";
+
+  ofstream mStaticModelFile;
+  mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
+  if (!mStaticModelFile.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  // Writing comments and function definition command
+  mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "(y, x, params)" << endl
+                   << "%" << endl
+                   << "% Status : Computes static model for Dynare" << endl
+                   << "%" << endl
+                   << "% Warning : this file is generated automatically by Dynare" << endl
+                   << "%           from model file (.mod)" << endl << endl;
+
+  writeStaticModel(mStaticModelFile);
+
+  mStaticModelFile.close();
+}
+
+void
+StaticModel::writeStaticCFile(const string &static_basename) const
+{
+  string filename = static_basename + ".c";
+
+  ofstream mStaticModelFile;
+  mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
+  if (!mStaticModelFile.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  mStaticModelFile << "/*" << endl
+                   << " * " << filename << " : Computes static model for Dynare" << 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;
+
+  // Writing the function Static
+  writeStaticModel(mStaticModelFile);
+
+  // Writing the gateway routine
+  mStaticModelFile << "/* 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;" << 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
+                   << "  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
+                   << "  /* Call the C Static. */" << endl
+                   << "  Static(y, x, params, residual, g1);" << endl
+                   << "}" << endl;
+
+  mStaticModelFile.close();
+}
+
+void
+StaticModel::writeStaticModel(ostream &StaticOutput) const
+{
+  ostringstream model_output;    // Used for storing model equations
+  ostringstream jacobian_output; // Used for storing jacobian equations
+  ostringstream hessian_output;
+  ostringstream lsymetric;       // For symmetric elements in hessian
+
+  ExprNodeOutputType output_type = (mode == eDLLMode ? oCStaticModel : oMatlabStaticModel);
+
+  writeModelLocalVariables(model_output, output_type);
+
+  writeTemporaryTerms(model_output, output_type);
+
+  writeModelEquations(model_output, output_type);
+
+  // Write Jacobian w.r. to endogenous only
+  for (first_derivatives_type::const_iterator it = first_derivatives.begin();
+       it != first_derivatives.end(); it++)
+    {
+      int eq = it->first.first;
+      int var = it->first.second;
+      NodeID d1 = it->second;
+
+      if (variable_table.getType(var) == eEndogenous)
+        {
+          ostringstream g1;
+          g1 << "  g1";
+          matrixHelper(g1, eq, symbol_table.getTypeSpecificID(variable_table.getSymbolID(var)), output_type);
+
+          jacobian_output << g1.str() << "=" << g1.str() << "+";
+          d1->writeOutput(jacobian_output, output_type, temporary_terms);
+          jacobian_output << ";" << endl;
+        }
+    }
+
+  // Write Hessian w.r. to endogenous only
+  if (computeStaticHessian)
+    for (second_derivatives_type::const_iterator it = second_derivatives.begin();
+         it != second_derivatives.end(); it++)
+      {
+        int eq = it->first.first;
+        int var1 = it->first.second.first;
+        int var2 = it->first.second.second;
+        NodeID d2 = it->second;
+
+        // Keep only derivatives w.r. to endogenous variables
+        if (variable_table.getType(var1) == eEndogenous
+            && variable_table.getType(var2) == eEndogenous)
+          {
+            int id1 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var1));
+            int id2 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var2));
+
+            int col_nb = id1*symbol_table.endo_nbr()+id2;
+            int col_nb_sym = id2*symbol_table.endo_nbr()+id1;
+
+            hessian_output << "  g2";
+            matrixHelper(hessian_output, eq, col_nb, output_type);
+            hessian_output << " = ";
+            d2->writeOutput(hessian_output, output_type, temporary_terms);
+            hessian_output << ";" << endl;
+
+            // Treating symetric elements
+            if (var1 != var2)
+              {
+                lsymetric <<  "  g2";
+                matrixHelper(lsymetric, eq, col_nb_sym, output_type);
+                lsymetric << " = " <<  "g2";
+                matrixHelper(lsymetric, eq, col_nb, output_type);
+                lsymetric << ";" << endl;
+              }
+          }
+      }
+
+  // Writing ouputs
+  if (mode != eDLLMode)
+    {
+      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)+2*imag(g1);" << endl
+                   << "  end" << endl
+                   << "end" << endl;
+      if (computeStaticHessian)
+        {
+          StaticOutput << "if nargout >= 3,\n";
+          // Writing initialization instruction for matrix g2
+          int ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
+          StaticOutput << "  g2 = sparse([],[],[], " << equations.size() << ", " << ncols << ", " << 5*ncols << ");" << endl
+                       << endl
+                       << "%" << endl
+                       << "% Hessian matrix" << endl
+                       << "%" << endl
+                       << endl
+                       << hessian_output.str()
+                       << lsymetric.str()
+                       << "end;" << endl;
+        }
+    }
+  else
+    {
+      StaticOutput << "void Static(double *y, double *x, double *params, double *residual, double *g1)" << endl
+                   << "{" << endl
+                   << "  double lhs, rhs;" << endl
+        // Writing residual equations
+                   << "  /* Residual equations */" << endl
+                   << "  if (residual == NULL)" << endl
+                   << "    return;" << endl
+                   << "  else" << endl
+                   << "    {" << endl
+                   << model_output.str()
+        // Writing Jacobian
+                   << "     /* Jacobian for endogenous variables without lag */" << endl
+                   << "     if (g1 == NULL)" << endl
+                   << "       return;" << endl
+                   << "     else" << endl
+                   << "       {" << endl
+                   << jacobian_output.str()
+                   << "       }" << endl
+                   << "    }" << endl
+                   << "}" << endl << endl;
+    }
+}
+
+void
+StaticModel::writeStaticFile(const string &basename) const
+{
+  switch (mode)
+    {
+    case eStandardMode:
+    case eSparseDLLMode:
+    case eSparseMode:
+      writeStaticMFile(basename + "_static");
+      break;
+    case eDLLMode:
+      writeStaticCFile(basename + "_static");
+      break;
+    }
+}
+
+void
+StaticModel::computingPass(bool no_tmp_terms)
+{
+  // Determine derivation order
+  int order = 1;
+  if (computeStaticHessian)
+    order = 2;
+
+  // Launch computations
+  cout << "Computing static model derivatives at order " << order << "...";
+  derive(order);
+  cout << " done" << endl;
+
+  if (!no_tmp_terms)
+    computeTemporaryTerms(order);
+}

StaticModel.hh

@@ -0,0 +1,51 @@
+/*
+ * Copyright (C) 2003-2009 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef _STATICMODEL_HH
+#define _STATICMODEL_HH
+
+using namespace std;
+
+#include "ModelTree.hh"
+
+//! Stores a static model
+class StaticModel : public ModelTree
+{
+private:
+  //! Writes the static model equations and its derivatives
+  /*! \todo handle hessian in C output */
+  void writeStaticModel(ostream &StaticOutput) const;
+  //! Writes static model file (Matlab version)
+  void writeStaticMFile(const string &static_basename) const;
+  //! Writes static model file (C version)
+  void writeStaticCFile(const string &static_basename) const;
+
+public:
+  StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
+  //! Whether static Hessian (w.r. to endogenous only) should be written
+  bool computeStaticHessian;
+  //! Execute computations (derivation)
+  /*! You must set computeStaticHessian before calling this function
+    \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
+  void computingPass(bool no_tmp_terms);
+  //! Writes static model file
+  void writeStaticFile(const string &basename) const;
+};
+
+#endif