ModelTree.cc 27 KB
Newer Older
1
#include <iostream>
sebastien's avatar
sebastien committed
2
3
4
#include <fstream>
#include <sstream>

5
6
#include "ModelTree.hh"
#include "Interface.hh"
7

sebastien's avatar
sebastien committed
8
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
sebastien's avatar
sebastien committed
9
10
                     NumericalConstants &num_constants_arg) :
  DataTree(symbol_table_arg, num_constants_arg),
sebastien's avatar
sebastien committed
11
  computeJacobian(false),
sebastien's avatar
sebastien committed
12
  computeJacobianExo(false),
sebastien's avatar
sebastien committed
13
14
  computeHessian(false),
  computeStaticHessian(false)
15
16
17
{
}

sebastien's avatar
sebastien committed
18
void
sebastien's avatar
sebastien committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
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
129
130
131
132
133
134
ModelTree::derive(int order)
{
  cout << "Processing derivation ..." << endl;

  cout << "  Processing Order 1... ";
  for(int var = 0; var < variable_table.size(); var++)
    for(int eq = 0; eq < (int) equations.size(); eq++)
      {
        NodeID d1 = equations[eq]->getDerivative(var);
        if (d1 == Zero)
          continue;
        first_derivatives[make_pair(eq, var)] = d1;
      }
  cout << "done" << endl;

  if (order == 2)
    {
      cout << "  Processing Order 2... ";
      for(first_derivatives_type::const_iterator it = first_derivatives.begin();
          it != first_derivatives.end(); it++)
        {
          int eq = it->first.first;
          int var1 = it->first.second;
          NodeID d1 = it->second;
      
          // Store only second derivatives with var2 <= var1
          for(int var2 = 0; var2 <= var1; var2++)
            {
              NodeID d2 = d1->getDerivative(var2);
              if (d2 == Zero)
                continue;
              second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
            }
        }
      cout << "done" << endl;
    }
}

void
ModelTree::computeTemporaryTerms(int order)
{
  map<NodeID, int> reference_count;
  temporary_terms.clear();

  for(vector<BinaryOpNode *>::iterator it = equations.begin();
      it != equations.end(); it++)
    (*it)->computeTemporaryTerms(reference_count, temporary_terms);

  for(first_derivatives_type::iterator it = first_derivatives.begin();
      it != first_derivatives.end(); it++)
    it->second->computeTemporaryTerms(reference_count, temporary_terms);

  if (order == 2)
    for(second_derivatives_type::iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
      it->second->computeTemporaryTerms(reference_count, temporary_terms);
}

void
ModelTree::writeTemporaryTerms(ostream &output, bool is_dynamic) const
{
  // A copy of temporary terms
  temporary_terms_type tt2;

  for(temporary_terms_type::const_iterator it = temporary_terms.begin();
      it != temporary_terms.end(); it++)
    {
      (*it)->writeOutput(output, is_dynamic, temporary_terms);
      output << " = ";

      (*it)->writeOutput(output, is_dynamic, tt2);

      // Insert current node into tt2
      tt2.insert(*it);

      output << ";" << endl;
    }
}

void
ModelTree::writeLocalParameters(ostream &output, bool is_dynamic) const
{
  for(map<int, NodeID>::const_iterator it = local_parameters_table.begin();
      it != local_parameters_table.end(); it++)
    {
      int id = it->first;
      NodeID value = it->second;
      output << symbol_table.getNameByID(eLocalParameter, id) << " = ";
      value->writeOutput(output, is_dynamic, temporary_terms);
      output << ";" << endl;
    }
}

void
ModelTree::writeModelEquations(ostream &output, bool is_dynamic) const
{
  for(int eq = 0; eq < (int) equations.size(); eq++)
    {
      BinaryOpNode *eq_node = equations[eq];

      NodeID lhs = eq_node->arg1;
      output << "lhs =";
      lhs->writeOutput(output, is_dynamic, temporary_terms);
      output << ";" << endl;

      NodeID rhs = eq_node->arg2;
      output << "rhs =";
      rhs->writeOutput(output, is_dynamic, temporary_terms);
      output << ";" << endl;

      output << "residual" << lpar << eq + 1 << rpar << "= lhs-rhs;" << endl;
    }
}

void
ModelTree::writeStaticMFile(const string &static_basename) const
135
{
sebastien's avatar
sebastien committed
136
  string filename = static_basename + interfaces::function_file_extension();
137

sebastien's avatar
sebastien committed
138
139
140
  ofstream mStaticModelFile;
  mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
  if (!mStaticModelFile.is_open())
141
    {
sebastien's avatar
sebastien committed
142
      cerr << "Error: Can't open file " << filename << " for writing" << endl;
143
144
      exit(-1);
    }
sebastien's avatar
sebastien committed
145
  // Writing comments and function definition command
sebastien's avatar
sebastien committed
146
  mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "( y, x )\n";
sebastien's avatar
sebastien committed
147
148
149
150
151
152
153
154
155
156
157
  mStaticModelFile << interfaces::comment()+"\n"+interfaces::comment();
  mStaticModelFile << "Status : Computes static model for Dynare\n" << interfaces::comment() << "\n";
  mStaticModelFile << interfaces::comment();
  mStaticModelFile << "Warning : this file is generated automatically by Dynare\n";
  mStaticModelFile << interfaces::comment();
  mStaticModelFile << "  from model file (.mod)\n\n";

  writeStaticModel(mStaticModelFile);

  interfaces::function_close();
  mStaticModelFile.close();
158
159
}

sebastien's avatar
sebastien committed
160
161

void
sebastien's avatar
sebastien committed
162
ModelTree::writeDynamicMFile(const string &dynamic_basename) const
163
{
sebastien's avatar
sebastien committed
164
165
166
167
168
  string filename = dynamic_basename + interfaces::function_file_extension();

  ofstream mDynamicModelFile;
  mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
  if (!mDynamicModelFile.is_open())
169
    {
sebastien's avatar
sebastien committed
170
      cerr << "Error: Can't open file " << filename << " for writing" << endl;
171
172
      exit(-1);
    }
sebastien's avatar
sebastien committed
173
174
175
176
177
178
179
180
181
182
183
184
  mDynamicModelFile << "function [residual, g1, g2] = " << dynamic_basename << "(y, x)\n";
  mDynamicModelFile << interfaces::comment()+"\n"+interfaces::comment();
  mDynamicModelFile << "Status : Computes dynamic model for Dynare\n" << interfaces::comment() << "\n";
  mDynamicModelFile << interfaces::comment();
  mDynamicModelFile << "Warning : this file is generated automatically by Dynare\n";
  mDynamicModelFile << interfaces::comment();
  mDynamicModelFile << "  from model file (.mod)\n\n";

  writeDynamicModel(mDynamicModelFile);  

  interfaces::function_close();
  mDynamicModelFile.close();
185
186
}

sebastien's avatar
sebastien committed
187
void
sebastien's avatar
sebastien committed
188
ModelTree::writeStaticCFile(const string &static_basename) const
189
{
sebastien's avatar
sebastien committed
190
191
192
193
194
  string filename = static_basename + ".c";

  ofstream mStaticModelFile;
  mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
  if (!mStaticModelFile.is_open())
195
    {
sebastien's avatar
sebastien committed
196
197
      cerr << "Error: Can't open file " << filename << " for writing" << endl;
      exit(-1);
198
    }
sebastien's avatar
sebastien committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
  mStaticModelFile << "/*\n";
  mStaticModelFile << " * " << filename << " : Computes static model for Dynare\n";
  mStaticModelFile << " * Warning : this file is generated automatically by Dynare\n";
  mStaticModelFile << " *           from model file (.mod)\n\n";
  mStaticModelFile << " */\n";
  mStaticModelFile << "#include <math.h>\n";
  mStaticModelFile << "#include \"mex.h\"\n";
  // A flobal variable for model parameters
  mStaticModelFile << "double *params;\n";

  // Writing the function Static
  writeStaticModel(mStaticModelFile);

  // Writing the gateway routine
  mStaticModelFile << "/* The gateway routine */\n";
  mStaticModelFile << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\n";
  mStaticModelFile << "{\n";
  mStaticModelFile << "  double *y, *x;\n";
  mStaticModelFile << "  double *residual, *g1;\n";
  mStaticModelFile << "  mxArray *M_;\n";
  mStaticModelFile << "\n";
  mStaticModelFile << "  /* Create a pointer to the input matrix y. */\n";
  mStaticModelFile << "  y = mxGetPr(prhs[0]);\n";
  mStaticModelFile << "\n";
  mStaticModelFile << "  /* Create a pointer to the input matrix x. */\n";
  mStaticModelFile << "  x = mxGetPr(prhs[1]);\n";
  mStaticModelFile << "\n";

  mStaticModelFile << "  residual = NULL;\n";
  mStaticModelFile << "  if (nlhs >= 1)\n";
  mStaticModelFile << "  {\n";
  mStaticModelFile << "      /* Set the output pointer to the output matrix residual. */\n";
sebastien's avatar
sebastien committed
231
  mStaticModelFile << "      plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);\n";
sebastien's avatar
sebastien committed
232
233
234
235
236
237
238
  mStaticModelFile << "     /* Create a C pointer to a copy of the output matrix residual. */\n";
  mStaticModelFile << "     residual = mxGetPr(plhs[0]);\n";
  mStaticModelFile << "  }\n\n";
  mStaticModelFile << "  g1 = NULL;\n";
  mStaticModelFile << "  if (nlhs >= 2)\n";
  mStaticModelFile << "  {\n";
  mStaticModelFile << "      /* Set the output pointer to the output matrix g1. */\n";
sebastien's avatar
sebastien committed
239
  mStaticModelFile << "      plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr << ", mxREAL);\n";
sebastien's avatar
sebastien committed
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  mStaticModelFile << "      /* Create a C pointer to a copy of the output matrix g1. */\n";
  mStaticModelFile << "      g1 = mxGetPr(plhs[1]);\n";
  mStaticModelFile << "  }\n\n";
  mStaticModelFile << "  /* Gets model parameters from global workspace of Matlab */\n";
  mStaticModelFile << "  M_ = mexGetVariable(\"global\",\"M_\");\n";
  mStaticModelFile << "  if (M_ == NULL ){\n";
  mStaticModelFile << "	    mexPrintf(\"Global variable not found : \");\n";
  mStaticModelFile << "	    mexErrMsgTxt(\"M_ \\n\");\n";
  mStaticModelFile << "  }\n";
  mStaticModelFile << "  params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));\n";
  mStaticModelFile << "  /* Call the C Static. */\n";
  mStaticModelFile << "  Static(y, x, residual, g1);\n";
  mStaticModelFile << "}\n";
  mStaticModelFile.close();
254
255
}

sebastien's avatar
sebastien committed
256
257
 
void
sebastien's avatar
sebastien committed
258
ModelTree::writeDynamicCFile(const string &dynamic_basename) const
259
{
sebastien's avatar
sebastien committed
260
261
262
263
264
  string filename = dynamic_basename + ".c";

  ofstream mDynamicModelFile;
  mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
  if (!mDynamicModelFile.is_open())
265
    {
sebastien's avatar
sebastien committed
266
267
      cerr << "Error: Can't open file " << filename << " for writing" << endl;
      exit(-1);
268
    }
sebastien's avatar
sebastien committed
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
  mDynamicModelFile << "/*\n";
  mDynamicModelFile << " * " << filename << " : Computes dynamic model for Dynare\n";
  mDynamicModelFile  << " *\n";
  mDynamicModelFile << " * Warning : this file is generated automatically by Dynare\n";
  mDynamicModelFile << " *           from model file (.mod)\n\n";
  mDynamicModelFile << " */\n";
  mDynamicModelFile << "#include <math.h>\n";
  mDynamicModelFile << "#include \"mex.h\"\n";
  // A flobal variable for model parameters
  mDynamicModelFile << "double *params;\n";
  // A global variable for it_
  mDynamicModelFile << "int it_;\n";
  mDynamicModelFile << "int nb_row_x;\n";

  // Writing the function body
  writeDynamicModel(mDynamicModelFile);  

  // Writing the gateway routine
  mDynamicModelFile << "/* The gateway routine */\n";
  mDynamicModelFile << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\n";
  mDynamicModelFile << "{\n";
  mDynamicModelFile << "  double *y, *x;\n";
  mDynamicModelFile << "  double *residual, *g1, *g2;\n";
  mDynamicModelFile << "  mxArray *M_;\n";
  mDynamicModelFile << "\n";
  mDynamicModelFile << "  /* Create a pointer to the input matrix y. */\n";
  mDynamicModelFile << "  y = mxGetPr(prhs[0]);\n";
  mDynamicModelFile << "\n";
  mDynamicModelFile << "  /* Create a pointer to the input matrix x. */\n";
  mDynamicModelFile << "  x = mxGetPr(prhs[1]);\n";
  mDynamicModelFile << "  /* Gets number of rows of matrix x. */\n";
  mDynamicModelFile << "  nb_row_x = mxGetM(prhs[1]);\n";
  mDynamicModelFile << "\n";
  mDynamicModelFile << "  residual = NULL;\n";
  mDynamicModelFile << "  if (nlhs >= 1)\n";
  mDynamicModelFile << "  {\n";
  mDynamicModelFile << "     /* Set the output pointer to the output matrix residual. */\n";
sebastien's avatar
sebastien committed
306
  mDynamicModelFile << "     plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);\n";
sebastien's avatar
sebastien committed
307
308
309
310
311
312
313
314
  mDynamicModelFile << "     /* Create a C pointer to a copy of the output matrix residual. */\n";
  mDynamicModelFile << "     residual = mxGetPr(plhs[0]);\n";
  mDynamicModelFile << "  }\n\n";
  mDynamicModelFile << "  g1 = NULL;\n";
  mDynamicModelFile << "  if (nlhs >= 2)\n";
  mDynamicModelFile << "  {\n";
  mDynamicModelFile << "     /* Set the output pointer to the output matrix g1. */\n";
  if (computeJacobianExo)
sebastien's avatar
sebastien committed
315
    mDynamicModelFile << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr() << ", mxREAL);\n";
sebastien's avatar
sebastien committed
316
  else if (computeJacobian)
sebastien's avatar
sebastien committed
317
    mDynamicModelFile << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.var_endo_nbr << ", mxREAL);\n";
sebastien's avatar
sebastien committed
318
319
320
321
322
323
324
  mDynamicModelFile << "     /* Create a C pointer to a copy of the output matrix g1. */\n";
  mDynamicModelFile << "     g1 = mxGetPr(plhs[1]);\n";
  mDynamicModelFile << "  }\n\n";
  mDynamicModelFile << "  g2 = NULL;\n";
  mDynamicModelFile << " if (nlhs >= 3)\n";
  mDynamicModelFile << "  {\n";
  mDynamicModelFile << "     /* Set the output pointer to the output matrix g2. */\n";
sebastien's avatar
sebastien committed
325
  mDynamicModelFile << "     plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr()*variable_table.get_dyn_var_nbr() << ", mxREAL);\n";
sebastien's avatar
sebastien committed
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
  mDynamicModelFile << "     /* Create a C pointer to a copy of the output matrix g1. */\n";
  mDynamicModelFile << "     g2 = mxGetPr(plhs[2]);\n";
  mDynamicModelFile << "  }\n\n";
  mDynamicModelFile << "  /* Gets model parameters from global workspace of Matlab */\n";
  mDynamicModelFile << "  M_ = mexGetVariable(\"global\",\"M_\");\n";
  mDynamicModelFile << "  if (M_ == NULL )\n";
  mDynamicModelFile << "  {\n";
  mDynamicModelFile << "	    mexPrintf(\"Global variable not found : \");\n";
  mDynamicModelFile << "	    mexErrMsgTxt(\"M_ \\n\");\n";
  mDynamicModelFile << "  }\n";
  mDynamicModelFile << "  params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));\n";
  mDynamicModelFile << "  /* Gets it_ from global workspace of Matlab */\n";
  mDynamicModelFile << "  it_ = (int) floor(mxGetScalar(mexGetVariable(\"global\", \"it_\")))-1;\n";
  mDynamicModelFile << "  /* Call the C subroutines. */\n";
  mDynamicModelFile << "  Dynamic(y, x, residual, g1, g2);\n";
  mDynamicModelFile << "}\n";
  mDynamicModelFile.close();
343
344
}

sebastien's avatar
sebastien committed
345
void
sebastien's avatar
sebastien committed
346
ModelTree::writeStaticModel(ostream &StaticOutput) const
347
348
349
{
  ostringstream model_output;    // Used for storing model equations
  ostringstream jacobian_output; // Used for storing jacobian equations
sebastien's avatar
sebastien committed
350
351
  ostringstream hessian_output;
  ostringstream lsymetric;       // For symmetric elements in hessian
352

sebastien's avatar
sebastien committed
353
  writeTemporaryTerms(model_output, false);
354

sebastien's avatar
sebastien committed
355
356
357
  writeLocalParameters(model_output, false);

  writeModelEquations(model_output, false);
358

sebastien's avatar
sebastien committed
359
  // Write Jacobian w.r. to endogenous only
sebastien's avatar
sebastien committed
360
361
  for(first_derivatives_type::const_iterator it = first_derivatives.begin();
      it != first_derivatives.end(); it++)
362
    {
sebastien's avatar
sebastien committed
363
364
365
366
367
      int eq = it->first.first;
      int var = it->first.second;
      NodeID d1 = it->second;

      if (variable_table.getType(var) == eEndogenous)
368
        {
sebastien's avatar
sebastien committed
369
370
371
372
373
374
          ostringstream g1;
          g1 << "  g1" << lpar << eq + 1 << ", " << variable_table.getSymbolID(var) + 1 << rpar;
          
          jacobian_output << g1.str() << "=" << g1.str() << "+";
          d1->writeOutput(jacobian_output, false, temporary_terms);
          jacobian_output << ";" << endl;
375
376
377
        }
    }

sebastien's avatar
sebastien committed
378
379
  // Write Hessian w.r. to endogenous only
  if (computeStaticHessian)
sebastien's avatar
sebastien committed
380
381
382
383
384
385
386
387
388
389
390
391
392
393
    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 = variable_table.getSymbolID(var1);
            int id2 = variable_table.getSymbolID(var2);
sebastien's avatar
sebastien committed
394

sebastien's avatar
sebastien committed
395
396
            int col_nb = id1*symbol_table.endo_nbr+id2+1;
            int col_nb_sym = id2*symbol_table.endo_nbr+id1+1;
sebastien's avatar
sebastien committed
397

sebastien's avatar
sebastien committed
398
399
400
            hessian_output << "  g2" << lpar << eq+1 << ", " << col_nb << rpar << " = ";
            d2->writeOutput(hessian_output, false, temporary_terms);
            hessian_output << ";" << endl;
sebastien's avatar
sebastien committed
401

sebastien's avatar
sebastien committed
402
403
404
405
406
            // Treating symetric elements
            if (var1 != var2)
              lsymetric <<  "  g2" << lpar << eq+1 << ", " << col_nb_sym << rpar << " = "
                        <<  "g2" << lpar << eq+1 << ", " << col_nb << rpar << ";" << endl;
          }
sebastien's avatar
sebastien committed
407

sebastien's avatar
sebastien committed
408
      }
sebastien's avatar
sebastien committed
409

410
411
412
413
414
415
  // Writing ouputs
  if (offset == 1)
    {
      StaticOutput << "global M_ \n";
      StaticOutput << "if M_.param_nbr > 0\n  params = M_.params;\nend\n";

sebastien's avatar
sebastien committed
416
      StaticOutput << "  residual = zeros( " << equations.size() << ", 1);\n";
417
418
419
420
421
422
423
424
425
      StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
      StaticOutput << "Model equations\n\t";
      StaticOutput << interfaces::comment() + "\n\n";
      StaticOutput << model_output.str();
      StaticOutput << "if ~isreal(residual)\n";
      StaticOutput << "  residual = real(residual)+imag(residual).^2;\n";
      StaticOutput << "end\n";
      StaticOutput << "if nargout >= 2,\n";
      StaticOutput << "  g1 = " <<
sebastien's avatar
sebastien committed
426
        "zeros(" << equations.size() << ", " <<
427
        symbol_table.endo_nbr << ");\n" ;
428
429
430
431
432
433
434
435
      StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
      StaticOutput << "Jacobian matrix\n\t";
      StaticOutput << interfaces::comment() + "\n\n";
      StaticOutput << jacobian_output.str();
      StaticOutput << "  if ~isreal(g1)\n";
      StaticOutput << "    g1 = real(g1)+2*imag(g1);\n";
      StaticOutput << "  end\n";
      StaticOutput << "end\n";
sebastien's avatar
sebastien committed
436
437
438
439
440
441
      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 = " <<
sebastien's avatar
sebastien committed
442
            "sparse([],[],[]," << equations.size() << ", " << ncols << ", " <<
sebastien's avatar
sebastien committed
443
444
445
446
447
448
449
            5*ncols << ");\n";
          StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
          StaticOutput << "Hessian matrix\n\t";
          StaticOutput << interfaces::comment() + "\n\n";
          StaticOutput << hessian_output.str() << lsymetric.str();
          StaticOutput << "end;\n";
        }
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
    }
  else
    {
      StaticOutput << "void Static(double *y, double *x, double *residual, double *g1)\n";
      StaticOutput << "{\n";
      StaticOutput << "  double lhs, rhs;\n\n";
      // Writing residual equations
      StaticOutput << "  /* Residual equations */\n";
      StaticOutput << "  if (residual == NULL) return;\n";
      StaticOutput << " {\n";
      StaticOutput << model_output.str();
      // Writing Jacobian
      StaticOutput << "   /* Jacobian for endogenous variables without lag */\n";
      StaticOutput << "   if (g1 == NULL) return;\n";
      StaticOutput << " {\n";
      StaticOutput << jacobian_output.str();
      StaticOutput << "  }\n";
      StaticOutput << " }\n";
      StaticOutput << "}\n\n";
    }
}

sebastien's avatar
sebastien committed
472
void
sebastien's avatar
sebastien committed
473
ModelTree::writeDynamicModel(ostream &DynamicOutput) const
474
{
sebastien's avatar
sebastien committed
475
  ostringstream lsymetric;       // Used when writing symetric elements in Hessian
476
477
478
479
  ostringstream model_output;    // Used for storing model equations
  ostringstream jacobian_output; // Used for storing jacobian equations
  ostringstream hessian_output;  // Used for storing Hessian equations

sebastien's avatar
sebastien committed
480
  writeTemporaryTerms(model_output, true);
481

sebastien's avatar
sebastien committed
482
483
484
485
486
  writeLocalParameters(model_output, true);

  writeModelEquations(model_output, true);

  // Writing Jacobian
487
  if (computeJacobian || computeJacobianExo)
sebastien's avatar
sebastien committed
488
489
490
491
492
493
    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;
494

sebastien's avatar
sebastien committed
495
496
497
498
499
500
501
502
503
504
        if (computeJacobianExo || variable_table.getType(var) == eEndogenous)
          {
            ostringstream g1;
            g1 << "  g1" << lpar << eq + 1 << ", " << variable_table.getSortID(var) + 1 << rpar;

            jacobian_output << g1.str() << "=" << g1.str() << "+";
            d1->writeOutput(jacobian_output, true, temporary_terms);
            jacobian_output << ";" << endl;
          }
      }
505

sebastien's avatar
sebastien committed
506
  // Writing Hessian
507
  if (computeHessian)
sebastien's avatar
sebastien committed
508
509
510
511
512
513
514
    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;
515

sebastien's avatar
sebastien committed
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
        int id1 = variable_table.getSortID(var1);
        int id2 = variable_table.getSortID(var2);

        int col_nb = id1*variable_table.get_dyn_var_nbr()+id2+1;
        int col_nb_sym = id2*variable_table.get_dyn_var_nbr()+id1+1;

        hessian_output << "  g2" << lpar << eq+1 << ", " << col_nb << rpar << " = ";
        d2->writeOutput(hessian_output, true, temporary_terms);
        hessian_output << ";" << endl;

        // Treating symetric elements
        if (id1 != id2)
          lsymetric <<  "  g2" << lpar << eq+1 << ", " << col_nb_sym << rpar << " = "
                    <<  "g2" << lpar << eq+1 << ", " << col_nb << rpar << ";" << endl;
      }

  int nrows = equations.size();
533
  int nvars = variable_table.var_endo_nbr;
sebastien's avatar
sebastien committed
534
  if (computeJacobianExo)
535
    nvars += symbol_table.exo_nbr + symbol_table.exo_det_nbr;
sebastien's avatar
sebastien committed
536

537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
  if (offset == 1)
    {
      DynamicOutput << "global M_ it_\n";
      DynamicOutput << "if M_.param_nbr > 0\n  params =  M_.params;\nend\n";
      DynamicOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
      DynamicOutput << "Model equations\n\t";
      DynamicOutput << interfaces::comment() + "\n\n";
      DynamicOutput << "residual = zeros(" << nrows << ", 1);\n";

      DynamicOutput << model_output.str();

      if (computeJacobian || computeJacobianExo)
        {
          DynamicOutput << "if nargout >= 2,\n";
          // Writing initialization instruction for matrix g1
          DynamicOutput << "  g1 = " <<
553
            "zeros(" << nrows << ", " << nvars << ");\n" ;
554
555
556
557
558
559
560
561
562
563
          DynamicOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
          DynamicOutput << "Jacobian matrix\n\t";
          DynamicOutput << interfaces::comment()+"\n\n";
          DynamicOutput << jacobian_output.str();
          DynamicOutput << "end\n";
        }
      if (computeHessian)
        {
          DynamicOutput << "if nargout >= 3,\n";
          // Writing initialization instruction for matrix g2
564
          int ncols = nvars*nvars;
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
          DynamicOutput << "  g2 = " <<
            "sparse([],[],[]," << nrows << ", " << ncols << ", " <<
            5*ncols << ");\n";
          DynamicOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
          DynamicOutput << "Hessian matrix\n\t";
          DynamicOutput << interfaces::comment() + "\n\n";
          DynamicOutput << hessian_output.str() << lsymetric.str();
          DynamicOutput << "end;\n";
        }
    }
  else
    {
      DynamicOutput << "void Dynamic(double *y, double *x, double *residual, double *g1, double *g2)\n";
      DynamicOutput << "{\n";
      DynamicOutput << "  double lhs, rhs;\n\n";
      DynamicOutput << "  /* Residual equations */\n";
      DynamicOutput << model_output.str();
      if (computeJacobian || computeJacobianExo)
        {
          DynamicOutput << "  /* Jacobian  */\n";
          DynamicOutput << "  if (g1 == NULL) return;\n";
          DynamicOutput << "  {\n";
          DynamicOutput << jacobian_output.str();
          DynamicOutput << "  }\n";
        }
      if (computeHessian)
        {
          DynamicOutput << "  /* Hessian for endogenous and exogenous variables */\n";
          DynamicOutput << "  if (g2 == NULL) return;\n";
          DynamicOutput << "   {\n";
          DynamicOutput << hessian_output.str() << lsymetric.str();
          DynamicOutput << "   }\n";
        }
      DynamicOutput << "}\n\n";
    }
}

sebastien's avatar
sebastien committed
602
void
sebastien's avatar
sebastien committed
603
ModelTree::writeOutput(ostream &output) const
604
605
606
607
608
609
610
611
612
613
614
{
  /* Writing initialisation for M_.lead_lag_incidence matrix
     M_.lead_lag_incidence is a matrix with as many columns as there are
     endogenous variables and as many rows as there are periods in the
     models (nbr of rows = M_.max_lag+M_.max_lead+1)

     The matrix elements are equal to zero if a variable isn't present in the
     model at a given period.
  */
  output << "M_.lead_lag_incidence = [";
  // Loop on endogenous variables
615
  for (int endoID = 0; endoID < symbol_table.endo_nbr; endoID++)
616
617
618
    {
      output << "\n\t";
      // Loop on periods
619
      for (int lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++)
620
621
        {
          // Getting name of symbol
622
          string name = symbol_table.getNameByID(eEndogenous, endoID);
623
          // and its variableID if exists with current period
624
          int varID = variable_table.getID(name, lag);
sebastien's avatar
sebastien committed
625
626
          if (varID >= 0)
            output << " " << variable_table.getPrintIndex(varID) + 1;
627
          else
sebastien's avatar
sebastien committed
628
            output << " 0";
629
630
631
632
633
634
        }
      output << ";";
    }
  output << "]';\n";

  // Writing initialization for some other variables
635
636
637
638
  output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr << "];\n";
  output << "M_.maximum_lag = " << variable_table.max_lag << ";\n";
  output << "M_.maximum_lead = " << variable_table.max_lead << ";\n";
  if (symbol_table.endo_nbr)
639
    {
640
641
642
      output << "M_.maximum_endo_lag = " << variable_table.max_endo_lag << ";\n";
      output << "M_.maximum_endo_lead = " << variable_table.max_endo_lead << ";\n";
      output << "oo_.steady_state = zeros(" << symbol_table.endo_nbr << ", 1);\n";
643
    }
644
  if (symbol_table.exo_nbr)
645
    {
646
647
648
      output << "M_.maximum_exo_lag = " << variable_table.max_exo_lag << ";\n";
      output << "M_.maximum_exo_lead = " << variable_table.max_exo_lead << ";\n";
      output << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr << ", 1);\n";
649
    }
650
  if (symbol_table.exo_det_nbr)
651
    {
652
653
654
      output << "M_.maximum_exo_det_lag = " << variable_table.max_exo_det_lag << ";\n";
      output << "M_.maximum_exo_det_lead = " << variable_table.max_exo_det_lead << ";\n";
      output << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr << ", 1);\n";
655
    }
656
  if (symbol_table.recur_nbr)
657
    {
658
659
660
      output << "M_.maximum_recur_lag = " << variable_table.max_recur_lag << ";\n";
      output << "M_.maximum_recur_lead = " << variable_table.max_recur_lead << ";\n";
      output << "oo_.recur_steady_state = zeros(" << symbol_table.recur_nbr << ", 1);\n";
661
    }
662
663
  if (symbol_table.parameter_nbr)
    output << "M_.params = zeros(" << symbol_table.parameter_nbr << ", 1);\n";
664
665
}

sebastien's avatar
sebastien committed
666
667
void
ModelTree::addEquation(NodeID eq)
668
{
sebastien's avatar
sebastien committed
669
670
671
  BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);

  if (beq == NULL || beq->op_code != oEqual)
672
    {
sebastien's avatar
sebastien committed
673
674
      cerr << "ModelTree::addEquation: you didn't provide an equal node!" << endl;
      exit(-1);
675
    }
sebastien's avatar
sebastien committed
676
677

  equations.push_back(beq);
678
}
sebastien's avatar
sebastien committed
679
680

void
sebastien's avatar
sebastien committed
681
682
683
ModelTree::checkPass() const
{
  // Exit if there is no equation in model file
sebastien's avatar
sebastien committed
684
  if (equations.size() == 0)
sebastien's avatar
sebastien committed
685
686
687
688
689
690
691
692
    {
      cerr << "No equation found in model file" << endl;
      exit(-1);
    }
}

void
ModelTree::computingPass()
sebastien's avatar
sebastien committed
693
{
sebastien's avatar
sebastien committed
694
695
  cout << equations.size() << " equation(s) found" << endl;

sebastien's avatar
sebastien committed
696
697
698
699
700
  // Sorting variable table
  variable_table.Sort();

  if (offset == 1)
    {
sebastien's avatar
sebastien committed
701
      min_cost = 40 * 90;
sebastien's avatar
sebastien committed
702
703
704
705
706
      lpar = '(';
      rpar = ')';
    }
  else
    {
sebastien's avatar
sebastien committed
707
      min_cost = 40 * 4;
sebastien's avatar
sebastien committed
708
709
710
      lpar = '[';
      rpar = ']';
    }
sebastien's avatar
sebastien committed
711

sebastien's avatar
sebastien committed
712
  if (computeHessian || computeStaticHessian)
sebastien's avatar
sebastien committed
713
714
715
716
    {
      derive(2);
      computeTemporaryTerms(2);
    }
sebastien's avatar
sebastien committed
717
  else
sebastien's avatar
sebastien committed
718
719
720
721
    {
      derive(1);
      computeTemporaryTerms(1);
    }
sebastien's avatar
sebastien committed
722
}
sebastien's avatar
sebastien committed
723

sebastien's avatar
sebastien committed
724
void
sebastien's avatar
sebastien committed
725
ModelTree::writeStaticFile(const string &basename) const
sebastien's avatar
sebastien committed
726
727
728
729
730
731
{
  if (offset)
    writeStaticMFile(basename + "_static");
  else
    writeStaticCFile(basename + "_static");
}
sebastien's avatar
sebastien committed
732

sebastien's avatar
sebastien committed
733
void
sebastien's avatar
sebastien committed
734
ModelTree::writeDynamicFile(const string &basename) const
sebastien's avatar
sebastien committed
735
736
737
{
  if (offset)
    writeDynamicMFile(basename + "_dynamic");
sebastien's avatar
sebastien committed
738
  else
sebastien's avatar
sebastien committed
739
    writeDynamicCFile(basename + "_dynamic");
sebastien's avatar
sebastien committed
740
}