- new Incidence_Matrix class

- lead and lag on exogenous variables - corrections in dr1_sparse and dr11_sparse - minor corrections in simulate git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2255 ac1d8469-bf42-47a9-8791-bf33cf982152
2008-11-14 16:07:47 +00:00 · 2008-11-14 16:07:47 +00:00 · a67a20bf66
parent 011944221c
commit a67a20bf66
9 changed files with 1177 additions and 671 deletions
--- a/BlockTriangular.cc
+++ b/BlockTriangular.cc
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@ -169,7 +169,11 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename) const
 {
  options_list.writeOutput(output);
  symbol_list.writeOutput("var_list_", output);
-  output << "info = stoch_simul(var_list_);\n";
+  output << "if(options_.model_mode)\n";
+  output << "  info = stoch_simul_sparse(var_list_);\n";
+  output << "else\n";
+  output << "  info = stoch_simul(var_list_);\n";
+  output << "end\n";
 }

 ForecastStatement::ForecastStatement(const SymbolList &symbol_list_arg,
--- a/ExprNode.cc
+++ b/ExprNode.cc
@ -164,6 +164,12 @@ NumConstNode::collectEndogenous(set<pair<int, int> > &result) const
 {
 }

+void
+NumConstNode::collectExogenous(set<pair<int, int> > &result) const
+{
+}
+
+
 VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg) :
  ExprNode(datatree_arg),
  symb_id(symb_id_arg),
@ -473,6 +479,14 @@ VariableNode::collectEndogenous(set<pair<int, int> > &result) const
    result.insert(make_pair(symb_id, lag));
 }

+void
+VariableNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  if (type == eExogenous)
+    result.insert(make_pair(symb_id, lag));
+}
+
+
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) :
  ExprNode(datatree_arg),
  arg(arg_arg),
@ -864,6 +878,13 @@ UnaryOpNode::collectEndogenous(set<pair<int, int> > &result) const
  arg->collectEndogenous(result);
 }

+void
+UnaryOpNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  arg->collectExogenous(result);
+}
+
+
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                           BinaryOpcode op_code_arg, const NodeID arg2_arg) :
  ExprNode(datatree_arg),
@ -1322,6 +1343,14 @@ BinaryOpNode::collectEndogenous(set<pair<int, int> > &result) const
  arg2->collectEndogenous(result);
 }

+
+void
+BinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  arg1->collectExogenous(result);
+  arg2->collectExogenous(result);
+}
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                           TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) :
  ExprNode(datatree_arg),
@ -1334,7 +1363,7 @@ TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,

  // Non-null derivatives are the union of those of the arguments
  // Compute set union of arg{1,2,3}->non_null_derivatives
-  set<int> non_null_derivatives_tmp;  
+  set<int> non_null_derivatives_tmp;
  set_union(arg1->non_null_derivatives.begin(),
            arg1->non_null_derivatives.end(),
            arg2->non_null_derivatives.begin(),
@ -1590,6 +1619,15 @@ TrinaryOpNode::collectEndogenous(set<pair<int, int> > &result) const
  arg3->collectEndogenous(result);
 }

+void
+TrinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  arg1->collectExogenous(result);
+  arg2->collectExogenous(result);
+  arg3->collectExogenous(result);
+}
+
+
 UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
                                         int symb_id_arg,
                                         const vector<NodeID> &arguments_arg) :
@ -1650,6 +1688,15 @@ UnknownFunctionNode::collectEndogenous(set<pair<int, int> > &result) const
    (*it)->collectEndogenous(result);
 }

+void
+UnknownFunctionNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  for(vector<NodeID>::const_iterator it = arguments.begin();
+      it != arguments.end(); it++)
+    (*it)->collectExogenous(result);
+}
+
+
 double
 UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (EvalException)
 {
--- a/ModFile.cc
+++ b/ModFile.cc
@ -69,7 +69,7 @@ ModFile::checkPass()
      exit(EXIT_FAILURE);
    }

-  if (mod_file_struct.simul_present && stochastic_statement_present)
+  if (mod_file_struct.simul_present && stochastic_statement_present && model_tree.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);
--- a/ModelTree.cc
+++ b/ModelTree.cc
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@ -169,7 +169,7 @@ ParsingDriver::add_model_variable(string *name, string *olag)
  if (type == eUnknownFunction)
    error("Symbol " + *name + " is a function name unknown to Dynare. It cannot be used inside model.");

-  if (type == eExogenous && lag != 0)
+  if (type == eExogenous && lag != 0 && (model_tree->mode != eSparseDLLMode && model_tree->mode != eSparseMode))
    warning("Exogenous variable " + *name + " has lead/lag " + *olag);

  if (type == eModelLocalVariable && lag != 0)
@ -179,9 +179,13 @@ ParsingDriver::add_model_variable(string *name, string *olag)

  NodeID id = model_tree->AddVariable(*name, lag);

-  if ((type == eEndogenous) && (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode))
-    model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
-
+  /*if (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode)
+    {
+      if (type == eEndogenous)
+        model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
+      if (type == eExogenous)
+        model_tree->block_triangular.fill_IM_X(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
+    }*/
  delete name;
  delete olag;
  return id;
@ -360,14 +364,16 @@ void
 ParsingDriver::sparse_dll()
 {
  model_tree->mode = eSparseDLLMode;
-  model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/
 }

 void
 ParsingDriver::sparse()
 {
  model_tree->mode = eSparseMode;
-  model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/
 }

 void
--- a/include/BlockTriangular.hh
+++ b/include/BlockTriangular.hh
@ -36,41 +36,61 @@ struct List_IM
  bool* IM;
 };

+
+//! create and manage the incidence matrix
+class IncidenceMatrix //: SymbolTable
+{
+  //friend class BlockTriangular;
+public:
+const SymbolTable &symbol_table;
+  IncidenceMatrix(const SymbolTable &symbol_table_arg);
+  List_IM* Build_IM(int lead_lag, SymbolType type);
+  List_IM* Get_IM(int lead_lag, SymbolType type) const;
+  bool* bGet_IM(int lead_lag, SymbolType type) const;
+  void fill_IM(int equation, int variable_endo, int lead_lag, SymbolType type);
+  void unfill_IM(int equation, int variable_endo, int lead_lag, SymbolType type);
+  void init_incidence_matrix();
+  void Free_IM() const;
+  List_IM* Get_First(SymbolType type) const;
+  void Print_IM(SymbolType type) const;
+  void Print_SIM(bool* IM, SymbolType type) const;
+
+  void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const;
+private:
+  List_IM *First_IM, *Last_IM, *First_IM_X, *Last_IM_X ;
+public:
+  int Model_Max_Lead, Model_Max_Lag;
+  int Model_Max_Lead_Endo, Model_Max_Lag_Endo, Model_Max_Lead_Exo, Model_Max_Lag_Exo;
+};
+
+
 //! Matrix of doubles for representing jacobian
 typedef map<pair<int ,int >,double> jacob_map;

 //! Create the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
 class BlockTriangular
 {
+  //friend class IncidenceMatrix;
 public:
-  BlockTriangular(const SymbolTable &symbol_table_arg);
  const SymbolTable &symbol_table;
+  BlockTriangular(const SymbolTable &symbol_table_arg);
+  //BlockTriangular(const IncidenceMatrix &incidence_matrix_arg);
+  //const SymbolTable &symbol_table;
  Blocks blocks;
  Normalization normalization;
-  List_IM* Build_IM(int lead_lag);
-  List_IM* Get_IM(int lead_lag);
-  bool* bGet_IM(int lead_lag) const;
-  void fill_IM(int equation, int variable_endo, int lead_lag);
-  void unfill_IM(int equation, int variable_endo, int lead_lag);
-  void init_incidence_matrix(int nb_endo);
-  void Print_IM(int n) const;
-  void Free_IM(List_IM* First_IM) const;
-  void Print_SIM(bool* IM, int n) const;
+  IncidenceMatrix incidencematrix;
  void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m);
  bool Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0 , jacob_map j_m);
  void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0);
-  void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n);
  void Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock);
  void Free_Block(Model_Block* ModelBlock) const;
-  List_IM *First_IM ;
-  List_IM *Last_IM ;
  simple *Index_Equ_IM;
  simple *Index_Var_IM;
  int prologue, epilogue;
-  int Model_Max_Lead, Model_Max_Lag, periods;
  bool bt_verbose;
-  int endo_nbr;
+  //int endo_nbr, exo_nbr;
  Model_Block* ModelBlock;
+  int periods;
  inline static std::string BlockType0(int type)
  {
    switch (type)
@ -98,14 +118,14 @@ public:
      {
      case EVALUATE_FORWARD:
      case EVALUATE_FORWARD_R:
-        return ("EVALUATE FORWARD            ");
+        return ("EVALUATE FORWARD             ");
        break;
      case EVALUATE_BACKWARD:
      case EVALUATE_BACKWARD_R:
        return ("EVALUATE BACKWARD            ");
        break;
      case SOLVE_FORWARD_SIMPLE:
-        return ("SOLVE FORWARD SIMPLE        ");
+        return ("SOLVE FORWARD SIMPLE         ");
        break;
      case SOLVE_BACKWARD_SIMPLE:
        return ("SOLVE BACKWARD SIMPLE        ");
@ -114,7 +134,7 @@ public:
        return ("SOLVE TWO BOUNDARIES SIMPLE  ");
        break;
      case SOLVE_FORWARD_COMPLETE:
-        return ("SOLVE FORWARD COMPLETE      ");
+        return ("SOLVE FORWARD COMPLETE       ");
        break;
      case SOLVE_BACKWARD_COMPLETE:
        return ("SOLVE BACKWARD COMPLETE      ");
--- a/include/ExprNode.hh
+++ b/include/ExprNode.hh
@ -143,6 +143,7 @@ public:
  /*! Endogenous are stored as integer pairs of the form (symb_id, lag)
      They are added to the set given in argument */
  virtual void collectEndogenous(set<pair<int, int> > &result) const = 0;
+  virtual void collectExogenous(set<pair<int, int> > &result) const = 0;
  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
                                     temporary_terms_type &temporary_terms,
                                     map<NodeID, int> &first_occurence,
@ -179,6 +180,7 @@ public:
  NumConstNode(DataTree &datatree_arg, int id_arg);
  virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
  virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) 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;
 };
@ -198,6 +200,7 @@ public:
  VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg);
  virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms = temporary_terms_type()) const;
  virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) 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;
 };
@ -223,6 +226,7 @@ public:
                                     Model_Block *ModelBlock,
                                     map_idx_type &map_idx) const;
  virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
  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;
@ -250,6 +254,7 @@ public:
                                     Model_Block *ModelBlock,
                                     map_idx_type &map_idx) const;
  virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
  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;
@ -282,6 +287,7 @@ public:
                                     Model_Block *ModelBlock,
                                     map_idx_type &map_idx) const;
  virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
  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;
@ -307,6 +313,7 @@ public:
                                     Model_Block *ModelBlock,
                                     map_idx_type &map_idx) const;
  virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) 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;
 };
@ -314,21 +321,22 @@ public:
 //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
 struct IM_compact
 {
-  int size, u_init, u_finish, nb_endo;
-  int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index;
+  int size, u_init, u_finish, nb_endo, size_exo;
+  int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index;
 };

 //! One block of the model
 struct Block
 {
-  int Size, Sized;
+  int Size, Sized, nb_exo, nb_exo_det;
  BlockType Type;
  BlockSimulationType Simulation_Type;
  int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo;
+  int Max_Lag_Endo, Max_Lead_Endo;
+  int Max_Lag_Exo, Max_Lead_Exo;
  bool is_linear;
  int *Equation, *Own_Derivative;
-  int *Variable, *Variable_Sorted, *dVariable;
-  int *variable_dyn_index, *variable_dyn_leadlag;
+  int *Variable, *Exogenous;
  temporary_terms_type *Temporary_terms;
  IM_compact *IM_lead_lag;
  int Code_Start, Code_Length;
@ -339,7 +347,7 @@ struct Model_Block
 {
  int Size, Periods;
  Block* Block_List;
-  int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block;
+  //int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block;
 };

 #endif
--- a/include/ModelTree.hh
+++ b/include/ModelTree.hh
@ -26,6 +26,7 @@ using namespace std;
 #include <vector>
 #include <map>
 #include <ostream>
+#include <algorithm>

 #include "SymbolTable.hh"
 #include "NumericalConstants.hh"
@ -82,12 +83,14 @@ private:
  //! 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;
+  void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, SymbolType type) 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