ModelTree.hh 20.5 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2011 Dynare Team
sebastien's avatar
sebastien committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * 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/>.
 */

20
21
#ifndef _MODELTREE_HH
#define _MODELTREE_HH
sebastien's avatar
sebastien committed
22
23
24

using namespace std;

25
26
#include <string>
#include <vector>
sebastien's avatar
sebastien committed
27
#include <deque>
sebastien's avatar
sebastien committed
28
#include <map>
sebastien's avatar
sebastien committed
29
#include <ostream>
sebastien's avatar
sebastien committed
30

31
#include "DataTree.hh"
sebastien's avatar
sebastien committed
32

33
34
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
typedef vector<pair<EquationType, expr_t > > equation_type_and_normalized_equation_t;
35
36

//! Vector describing variables: max_lag in the block, max_lead in the block
37
typedef vector<pair< int, int> > lag_lead_vector_t;
38
39

//! for each block contains pair< pair<Simulation_Type, first_equation>, pair < Block_Size, Recursive_part_Size > >
40
typedef vector<pair< pair< BlockSimulationType, int>, pair<int, int> > > block_type_firstequation_size_mfs_t;
41

42
43
//! for a block contains derivatives pair< pair<block_equation_number, block_variable_number> , pair<lead_lag, expr_t> >
typedef vector< pair<pair<int, int>, pair< int, expr_t > > > block_derivatives_equation_variable_laglead_nodeid_t;
44
45

//! for all blocks derivatives description
46
typedef vector<block_derivatives_equation_variable_laglead_nodeid_t> blocks_derivatives_t;
47

sebastien's avatar
sebastien committed
48
//! Shared code for static and dynamic models
49
50
class ModelTree : public DataTree
{
51
52
  friend class DynamicModel;
  friend class StaticModel;
sebastien's avatar
sebastien committed
53
protected:
sebastien's avatar
sebastien committed
54
  //! Stores declared and generated auxiliary equations
sebastien's avatar
sebastien committed
55
56
  vector<BinaryOpNode *> equations;

sebastien's avatar
sebastien committed
57
58
59
  //! Only stores generated auxiliary equations, in an order meaningful for evaluation
  deque<BinaryOpNode *> aux_equations;

60
61
62
  //! Stores equation tags
  vector<pair<int, pair<string, string> > > equation_tags;

63
64
65
  //! Number of non-zero derivatives
  int NNZDerivatives[3];

66
  typedef map<pair<int, int>, expr_t> first_derivatives_t;
sebastien's avatar
sebastien committed
67
68
69
  //! First order derivatives
  /*! First index is equation number, second is variable w.r. to which is computed the derivative.
    Only non-null derivatives are stored in the map.
70
    Variable indices are those of the getDerivID() method.
71
  */
72
  first_derivatives_t first_derivatives;
73

74
  typedef map<pair<int, pair<int, int> >, expr_t> second_derivatives_t;
sebastien's avatar
sebastien committed
75
76
77
78
  //! Second order derivatives
  /*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
    Only non-null derivatives are stored in the map.
    Contains only second order derivatives where var1 >= var2 (for obvious symmetry reasons).
79
    Variable indices are those of the getDerivID() method.
sebastien's avatar
sebastien committed
80
  */
81
  second_derivatives_t second_derivatives;
sebastien's avatar
sebastien committed
82

83
  typedef map<pair<int, pair<int, pair<int, int> > >, expr_t> third_derivatives_t;
sebastien's avatar
sebastien committed
84
85
86
87
  //! Third order derivatives
  /*! First index is equation number, second, third and fourth are variables w.r. to which is computed the derivative.
    Only non-null derivatives are stored in the map.
    Contains only third order derivatives where var1 >= var2 >= var3 (for obvious symmetry reasons).
88
    Variable indices are those of the getDerivID() method.
sebastien's avatar
sebastien committed
89
  */
90
  third_derivatives_t third_derivatives;
sebastien's avatar
sebastien committed
91

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
  //! Derivatives of the residuals w.r. to parameters
  /*! First index is equation number, second is parameter.
    Only non-null derivatives are stored in the map.
    Parameter indices are those of the getDerivID() method.
  */
  first_derivatives_t residuals_params_derivatives;

  //! Second derivatives of the residuals w.r. to parameters
  /*! First index is equation number, second and third indeces are parameters.
    Only non-null derivatives are stored in the map.
    Parameter indices are those of the getDerivID() method.
  */
  second_derivatives_t residuals_params_second_derivatives;

  //! Derivatives of the jacobian w.r. to parameters
  /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter.
    Only non-null derivatives are stored in the map.
    Variable and parameter indices are those of the getDerivID() method.
  */
  second_derivatives_t jacobian_params_derivatives;

  //! Second derivatives of the jacobian w.r. to parameters
  /*! First index is equation number, second is endo/exo/exo_det variable, and third and fourth are parameters.
    Only non-null derivatives are stored in the map.
    Variable and parameter indices are those of the getDerivID() method.
  */
  third_derivatives_t jacobian_params_second_derivatives;

  //! Derivatives of the hessian w.r. to parameters
  /*! First index is equation number, first and second are endo/exo/exo_det variable, and third is parameter.
    Only non-null derivatives are stored in the map.
    Variable and parameter indices are those of the getDerivID() method.
  */
  third_derivatives_t hessian_params_derivatives;


  //! Temporary terms for the static/dynamic file (those which will be noted Txxxx)
129
  temporary_terms_t temporary_terms;
130
131
132
133
134

  //! Temporary terms for the file containing parameters derivatives
  temporary_terms_t params_derivs_temporary_terms;


135
  //! Trend variables and their growth factors
136
137
138
139
140
  map<int, expr_t> trend_symbols_map;

  //! for all trends; the boolean is true if this is a log-trend, false otherwise
  typedef map<int, pair<bool, expr_t> > nonstationary_symbols_map_t;

141
  //! Nonstationary variables and their deflators
142
  nonstationary_symbols_map_t nonstationary_symbols_map;
sebastien's avatar
sebastien committed
143

144
145
146
147
148
149
  //! vector of block reordered variables and equations
  vector<int> equation_reordered, variable_reordered, inv_equation_reordered, inv_variable_reordered;

  //! the file containing the model and the derivatives code
  ofstream code_file;

150
151
152
153
154
155
156
157
158
  //! Computes 1st derivatives
  /*! \param vars the derivation IDs w.r. to which compute the derivatives */
  void computeJacobian(const set<int> &vars);
  //! Computes 2nd derivatives
  /*! \param vars the derivation IDs w.r. to which derive the 1st derivatives */
  void computeHessian(const set<int> &vars);
  //! Computes 3rd derivatives
  /*! \param vars the derivation IDs w.r. to which derive the 2nd derivatives */
  void computeThirdDerivatives(const set<int> &vars);
159
160
  //! Computes derivatives of the Jacobian and Hessian w.r. to parameters
  void computeParamsDerivatives();
161

sebastien's avatar
sebastien committed
162
  //! Write derivative of an equation w.r. to a variable
163
  void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
164
  //! Computes temporary terms (for all equations and derivatives)
165
  void computeTemporaryTerms(bool is_matlab);
166
167
168
  //! Computes temporary terms for the file containing parameters derivatives
  void computeParamsDerivativesTemporaryTerms();
//! Writes temporary terms
169
  void writeTemporaryTerms(const temporary_terms_t &tt, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
170
  //! Compiles temporary terms
171
  void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
172
173
174
  //! Adds informations for simulation in a binary file
  void Write_Inf_To_Bin_File(const string &basename, int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const;

sebastien's avatar
sebastien committed
175
  //! Writes model local variables
176
  /*! No temporary term is used in the output, so that local parameters declarations can be safely put before temporary terms declaration in the output files */
177
  void writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
sebastien's avatar
sebastien committed
178
  //! Writes model equations
sebastien's avatar
sebastien committed
179
  void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const;
180
  //! Compiles model equations
181
  void compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
sebastien's avatar
sebastien committed
182

183
184
185
  //! Writes LaTeX model file
  void writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const;

186
187
  //! Sparse matrix of double to store the values of the Jacobian
  /*! First index is equation number, second index is endogenous type specific ID */
188
  typedef map<pair<int, int>, double> jacob_map_t;
189
190
191

  //! Sparse matrix of double to store the values of the Jacobian
  /*! First index is lag, second index is equation number, third index is endogenous type specific ID */
192
  typedef map<pair<int, pair<int, int> >, expr_t> dynamic_jacob_map_t;
193
194
195
196
197
198
199
200
201

  //! Normalization of equations
  /*! Maps endogenous type specific IDs to equation numbers */
  vector<int> endo2eq;

  //! number of equation in the prologue and in the epilogue
  unsigned int epilogue, prologue;

  //! for each block contains pair< max_lag, max_lead>
202
  lag_lead_vector_t block_lag_lead;
203
204

  //! Compute the matching between endogenous and variable using the jacobian contemporaneous_jacobian
sebastien's avatar
sebastien committed
205
206
207
208
  /*!
    \param contemporaneous_jacobian Jacobian used as an incidence matrix: all elements declared in the map (even if they are zero), are used as vertices of the incidence matrix
    \return True if a complete normalization has been achieved
  */
209
  bool computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose);
sebastien's avatar
sebastien committed
210

211
  //! Try to compute the matching between endogenous and variable using a decreasing cutoff
sebastien's avatar
sebastien committed
212
213
214
215
216
  /*!
    Applied to the jacobian contemporaneous_jacobian and stop when a matching is found.
    If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e. use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to reflect all the edges in the symbolic incidence matrix.
    If no matching is found with a zero cutoff close to zero an error message is printout.
  */
217
  void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian);
sebastien's avatar
sebastien committed
218

219
220
221
  //! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
  void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
  //! Evaluate the jacobian and suppress all the elements below the cutoff
222
  void evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose);
223
  //! Search the equations and variables belonging to the prologue and the epilogue of the model
224
  void computePrologueAndEpilogue(const jacob_map_t &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered);
225
  //! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
226
  equation_type_and_normalized_equation_t equationTypeDetermination(const map<pair<int, pair<int, int> >, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
227
  //! Compute the block decomposition and for a non-recusive block find the minimum feedback set
228
  void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const;
229
  //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
230
  block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int> > > &block_col_type);
231
  //! Determine the maximum number of lead and lag for the endogenous variable in a bloc
232
  void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const;
233
  //! Print an abstract of the block structure of the model
234
  void printBlockDecomposition(const vector<pair<int, int> > &blocks) const;
235
  //! Determine for each block if it is linear or not
236
  vector<bool> BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector<int> &variable_reordered) const;
237
238
239
240
241
242
243
244
245

  //! Determine the simulation type of each block
  virtual BlockSimulationType getBlockSimulationType(int block_number) const = 0;
  //! Return the number of blocks
  virtual unsigned int getNbBlocks() const = 0;
  //! Return the first equation number of a block
  virtual unsigned int getBlockFirstEquation(int block_number) const = 0;
  //! Return the size of the block block_number
  virtual unsigned int getBlockSize(int block_number) const = 0;
246
247
248
249
  //! Return the number of exogenous variable in the block block_number
  virtual unsigned int getBlockExoSize(int block_number) const = 0;
  //! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
  virtual unsigned int getBlockExoColSize(int block_number) const = 0;
250
251
252
253
254
255
  //! Return the number of feedback variable of the block block_number
  virtual unsigned int getBlockMfs(int block_number) const = 0;
  //! Return the maximum lag in a block
  virtual unsigned int getBlockMaxLag(int block_number) const = 0;
  //! Return the maximum lead in a block
  virtual unsigned int getBlockMaxLead(int block_number) const = 0;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
256
257
258
259
260
  inline void setBlockLeadLag(int block, int max_lag, int max_lead) 
    {
       block_lag_lead[block] = make_pair(max_lag, max_lead);
    };
  
261
262
263
264
  //! Return the type of equation (equation_number) belonging to the block block_number
  virtual EquationType getBlockEquationType(int block_number, int equation_number) const = 0;
  //! Return true if the equation has been normalized
  virtual bool isBlockEquationRenormalized(int block_number, int equation_number) const = 0;
265
266
267
268
  //! Return the expr_t of the equation equation_number belonging to the block block_number
  virtual expr_t getBlockEquationExpr(int block_number, int equation_number) const = 0;
  //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
  virtual expr_t getBlockEquationRenormalizedExpr(int block_number, int equation_number) const = 0;
269
270
271
272
  //! Return the original number of equation equation_number belonging to the block block_number
  virtual int getBlockEquationID(int block_number, int equation_number) const = 0;
  //! Return the original number of variable variable_number belonging to the block block_number
  virtual int getBlockVariableID(int block_number, int variable_number) const = 0;
273
274
  //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
  virtual int getBlockVariableExoID(int block_number, int variable_number) const = 0;
275
276
277
278
  //! Return the position of equation_number in the block number belonging to the block block_number
  virtual int getBlockInitialEquationID(int block_number, int equation_number) const = 0;
  //! Return the position of variable_number in the block number belonging to the block block_number
  virtual int getBlockInitialVariableID(int block_number, int variable_number) const = 0;
279
280
281
282
283
284
  //! Return the position of variable_number in the block number belonging to the block block_number
  virtual int getBlockInitialExogenousID(int block_number, int variable_number) const = 0;
  //! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
  virtual int getBlockInitialDetExogenousID(int block_number, int variable_number) const = 0;
  //! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
  virtual int getBlockInitialOtherEndogenousID(int block_number, int variable_number) const = 0;
285
286
  //! Initialize equation_reordered & variable_reordered
  void initializeVariablesAndEquations();
sebastien's avatar
sebastien committed
287
public:
288
  ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_arg);
289
290
291
292
293
294
295
296
297
  //! Absolute value under which a number is considered to be zero
  double cutoff;
  //! Compute the minimum feedback set
  /*!   0 : all endogenous variables are considered as feedback variables
    1 : the variables belonging to non normalized equation are considered as feedback variables
    2 : the variables belonging to a non linear equation are considered as feedback variables
    3 : the variables belonging to a non normalizable non linear equation are considered as feedback variables
    default value = 0 */
  int mfs;
sebastien's avatar
sebastien committed
298
  //! Declare a node as an equation of the model
299
  void addEquation(expr_t eq);
300
301
  //! Declare a node as an equation of the model, also giving its tags
  void addEquation(expr_t eq, vector<pair<string, string> > &eq_tags);
sebastien's avatar
sebastien committed
302
  //! Declare a node as an auxiliary equation of the model, adding it at the end of the list of auxiliary equations
303
  void addAuxEquation(expr_t eq);
304
  //! Returns the number of equations in the model
sebastien's avatar
sebastien committed
305
  int equation_number() const;
306
307
  //! Adds a trend variable with its growth factor
  void addTrendVariables(vector<int> trend_vars, expr_t growth_factor) throw (TrendException);
308
309
  //! Adds a nonstationary variables with their (common) deflator
  void addNonstationaryVariables(vector<int> nonstationary_vars, bool log_deflator, expr_t deflator) throw (TrendException);
310
  void set_cutoff_to_zero();
311
312
313
314
315
316
317
318
  //! Helper for writing the Jacobian elements in MATLAB and C
  /*! Writes either (i+1,j+1) or [i+j*no_eq] */
  void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
  //! Helper for writing the sparse Hessian or third derivatives in MATLAB and C
  /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[1]]
    If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[2]] */
  void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;

319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
  inline static std::string
  c_Equation_Type(int type)
  {
    char c_Equation_Type[4][13] =
      {
        "E_UNKNOWN   ",
        "E_EVALUATE  ",
        "E_EVALUATE_S",
        "E_SOLVE     "
      };
    return (c_Equation_Type[type]);
  };

  inline static std::string
  BlockType0(BlockType type)
  {
    switch (type)
      {
      case SIMULTANS:
        return ("SIMULTANEOUS TIME SEPARABLE  ");
        break;
      case PROLOGUE:
        return ("PROLOGUE                     ");
        break;
      case EPILOGUE:
        return ("EPILOGUE                     ");
        break;
      case SIMULTAN:
        return ("SIMULTANEOUS TIME UNSEPARABLE");
        break;
      default:
        return ("UNKNOWN                      ");
        break;
      }
  };

  inline static std::string
  BlockSim(int type)
  {
    switch (type)
      {
      case EVALUATE_FORWARD:
        return ("EVALUATE FORWARD             ");
        break;
      case EVALUATE_BACKWARD:
        return ("EVALUATE BACKWARD            ");
        break;
      case SOLVE_FORWARD_SIMPLE:
        return ("SOLVE FORWARD SIMPLE         ");
        break;
      case SOLVE_BACKWARD_SIMPLE:
        return ("SOLVE BACKWARD SIMPLE        ");
        break;
      case SOLVE_TWO_BOUNDARIES_SIMPLE:
        return ("SOLVE TWO BOUNDARIES SIMPLE  ");
        break;
      case SOLVE_FORWARD_COMPLETE:
        return ("SOLVE FORWARD COMPLETE       ");
        break;
      case SOLVE_BACKWARD_COMPLETE:
        return ("SOLVE BACKWARD COMPLETE      ");
        break;
      case SOLVE_TWO_BOUNDARIES_COMPLETE:
        return ("SOLVE TWO BOUNDARIES COMPLETE");
        break;
      default:
        return ("UNKNOWN                      ");
        break;
      }
  };
389
};
sebastien's avatar
sebastien committed
390

391
#endif