ExprNode.cc 173 KB
Newer Older
1
/*
2
 * Copyright (C) 2007-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/>.
 */

sebastien's avatar
sebastien committed
20
21
22
23
#include <iostream>
#include <iterator>
#include <algorithm>

24
25
26
27
28
// For select1st()
#ifdef __GNUC__
# include <ext/functional>
using namespace __gnu_cxx;
#endif
sebastien's avatar
sebastien committed
29

30
#include <cassert>
31
#include <cmath>
sebastien's avatar
sebastien committed
32

sebastien's avatar
sebastien committed
33
34
#include "ExprNode.hh"
#include "DataTree.hh"
35
#include "ModFile.hh"
sebastien's avatar
sebastien committed
36

sebastien's avatar
sebastien committed
37
ExprNode::ExprNode(DataTree &datatree_arg) : datatree(datatree_arg), preparedForDerivation(false)
sebastien's avatar
sebastien committed
38
39
40
41
42
43
44
45
46
47
48
49
{
  // Add myself to datatree
  datatree.node_list.push_back(this);

  // Set my index and increment counter
  idx = datatree.node_counter++;
}

ExprNode::~ExprNode()
{
}

50
expr_t
51
ExprNode::getDerivative(int deriv_id)
sebastien's avatar
sebastien committed
52
{
sebastien's avatar
sebastien committed
53
54
  if (!preparedForDerivation)
    prepareForDerivation();
55

sebastien's avatar
sebastien committed
56
  // Return zero if derivative is necessarily null (using symbolic a priori)
57
  set<int>::const_iterator it = non_null_derivatives.find(deriv_id);
sebastien's avatar
sebastien committed
58
59
  if (it == non_null_derivatives.end())
    return datatree.Zero;
60

sebastien's avatar
sebastien committed
61
  // If derivative is stored in cache, use the cached value, otherwise compute it (and cache it)
62
  map<int, expr_t>::const_iterator it2 = derivatives.find(deriv_id);
sebastien's avatar
sebastien committed
63
64
65
66
  if (it2 != derivatives.end())
    return it2->second;
  else
    {
67
      expr_t d = computeDerivative(deriv_id);
68
      derivatives[deriv_id] = d;
sebastien's avatar
sebastien committed
69
70
71
72
      return d;
    }
}

ferhat's avatar
ferhat committed
73
int
74
ExprNode::precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const
75
76
77
78
{
  // For a constant, a variable, or a unary op, the precedence is maximal
  return 100;
}
ferhat's avatar
ferhat committed
79

sebastien's avatar
sebastien committed
80
int
81
ExprNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
82
83
84
85
{
  // For a terminal node, the cost is null
  return 0;
}
sebastien's avatar
sebastien committed
86

sebastien's avatar
sebastien committed
87
88
89
90
91
void
ExprNode::collectEndogenous(set<pair<int, int> > &result) const
{
  set<pair<int, int> > symb_ids;
  collectVariables(eEndogenous, symb_ids);
92
93
  for (set<pair<int, int> >::const_iterator it = symb_ids.begin();
       it != symb_ids.end(); it++)
sebastien's avatar
sebastien committed
94
95
96
97
98
99
100
101
    result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(it->first), it->second));
}

void
ExprNode::collectExogenous(set<pair<int, int> > &result) const
{
  set<pair<int, int> > symb_ids;
  collectVariables(eExogenous, symb_ids);
102
103
  for (set<pair<int, int> >::const_iterator it = symb_ids.begin();
       it != symb_ids.end(); it++)
sebastien's avatar
sebastien committed
104
105
106
107
108
109
110
111
112
    result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(it->first), it->second));
}

void
ExprNode::collectModelLocalVariables(set<int> &result) const
{
  set<pair<int, int> > symb_ids;
  collectVariables(eModelLocalVariable, symb_ids);
  transform(symb_ids.begin(), symb_ids.end(), inserter(result, result.begin()),
113
            select1st<pair<int, int> >());
sebastien's avatar
sebastien committed
114
115
}

sebastien's avatar
sebastien committed
116
void
117
ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
118
                                temporary_terms_t &temporary_terms,
sebastien's avatar
sebastien committed
119
                                bool is_matlab) const
120
121
122
{
  // Nothing to do for a terminal node
}
sebastien's avatar
sebastien committed
123

sebastien's avatar
sebastien committed
124
void
125
ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
126
                                temporary_terms_t &temporary_terms,
127
                                map<expr_t, pair<int, int> > &first_occurence,
sebastien's avatar
sebastien committed
128
                                int Curr_block,
129
                                vector<vector<temporary_terms_t> > &v_temporary_terms,
130
                                int equation) const
131
132
133
{
  // Nothing to do for a terminal node
}
sebastien's avatar
sebastien committed
134

135
136
pair<int, expr_t >
ExprNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
137
{
138
  /* nothing to do */
139
  return (make_pair(0, (expr_t) NULL));
140
}
141

sebastien's avatar
sebastien committed
142
void
143
ExprNode::writeOutput(ostream &output) const
sebastien's avatar
sebastien committed
144
{
145
  writeOutput(output, oMatlabOutsideModel, temporary_terms_t());
sebastien's avatar
sebastien committed
146
147
}

148
149
150
void
ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type) const
{
151
  writeOutput(output, output_type, temporary_terms_t());
152
153
}

154
void
155
ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const
156
{
157
  deriv_node_temp_terms_t tef_terms;
158
159
160
  writeOutput(output, output_type, temporary_terms, tef_terms);
}

161
162
void
ExprNode::compile(ostream &CompileCode, unsigned int &instruction_number,
163
164
                  bool lhs_rhs, const temporary_terms_t &temporary_terms,
                  const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
165
166
167
168
169
{
  deriv_node_temp_terms_t tef_terms;
  compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
}

170
171
void
ExprNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
172
173
                                      const temporary_terms_t &temporary_terms,
                                      deriv_node_temp_terms_t &tef_terms) const
174
175
176
177
{
  // Nothing to do
}

178
179
void
ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
180
181
182
                                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
                                        const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
                                        deriv_node_temp_terms_t &tef_terms) const
183
184
185
186
{
  // Nothing to do
}

sebastien's avatar
sebastien committed
187
VariableNode *
188
ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
sebastien's avatar
sebastien committed
189
190
191
192
193
194
195
196
{
  int n = maxEndoLead();
  assert(n >= 2);

  subst_table_t::const_iterator it = subst_table.find(this);
  if (it != subst_table.end())
    return const_cast<VariableNode *>(it->second);

197
  expr_t substexpr = decreaseLeadsLags(n-1);
sebastien's avatar
sebastien committed
198
199
200
201
  int lag = n-2;

  // Each iteration tries to create an auxvar such that auxvar(+1)=expr(-lag)
  // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to expr(-lag-1) (resp. expr(-lag))
202
  while (lag >= 0)
sebastien's avatar
sebastien committed
203
    {
204
      expr_t orig_expr = decreaseLeadsLags(lag);
sebastien's avatar
sebastien committed
205
206
207
      it = subst_table.find(orig_expr);
      if (it == subst_table.end())
        {
208
          int symb_id = datatree.symbol_table.addEndoLeadAuxiliaryVar(orig_expr->idx);
sebastien's avatar
sebastien committed
209
210
211
212
213
214
215
216
217
218
219
220
221
222
          neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
          substexpr = datatree.AddVariable(symb_id, +1);
          assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
          subst_table[orig_expr] = dynamic_cast<VariableNode *>(substexpr);
        }
      else
        substexpr = const_cast<VariableNode *>(it->second);

      lag--;
    }

  return dynamic_cast<VariableNode *>(substexpr);
}

sebastien's avatar
sebastien committed
223
224
225
226
227
228
229
230
231
232
VariableNode *
ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
{
  int n = maxExoLead();
  assert(n >= 1);

  subst_table_t::const_iterator it = subst_table.find(this);
  if (it != subst_table.end())
    return const_cast<VariableNode *>(it->second);

233
  expr_t substexpr = decreaseLeadsLags(n);
sebastien's avatar
sebastien committed
234
235
236
237
  int lag = n-1;

  // Each iteration tries to create an auxvar such that auxvar(+1)=expr(-lag)
  // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to expr(-lag-1) (resp. expr(-lag))
238
  while (lag >= 0)
sebastien's avatar
sebastien committed
239
    {
240
      expr_t orig_expr = decreaseLeadsLags(lag);
sebastien's avatar
sebastien committed
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
      it = subst_table.find(orig_expr);
      if (it == subst_table.end())
        {
          int symb_id = datatree.symbol_table.addExoLeadAuxiliaryVar(orig_expr->idx);
          neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
          substexpr = datatree.AddVariable(symb_id, +1);
          assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
          subst_table[orig_expr] = dynamic_cast<VariableNode *>(substexpr);
        }
      else
        substexpr = const_cast<VariableNode *>(it->second);

      lag--;
    }

  return dynamic_cast<VariableNode *>(substexpr);
}

259
260
261
262
263
264
265
266
267
268
269
270
bool
ExprNode::isNumConstNodeEqualTo(double value) const
{
  return false;
}

bool
ExprNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
  return false;
}

sebastien's avatar
sebastien committed
271
NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
272
273
  ExprNode(datatree_arg),
  id(id_arg)
sebastien's avatar
sebastien committed
274
275
276
{
  // Add myself to the num const map
  datatree.num_const_node_map[id] = this;
sebastien's avatar
sebastien committed
277
}
sebastien's avatar
sebastien committed
278

sebastien's avatar
sebastien committed
279
280
281
282
void
NumConstNode::prepareForDerivation()
{
  preparedForDerivation = true;
sebastien's avatar
sebastien committed
283
284
285
  // All derivatives are null, so non_null_derivatives is left empty
}

286
expr_t
287
NumConstNode::computeDerivative(int deriv_id)
sebastien's avatar
sebastien committed
288
289
290
291
{
  return datatree.Zero;
}

292
void
293
NumConstNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
294
{
295
  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
296
297
298
  if (it != temporary_terms.end())
    temporary_terms_inuse.insert(idx);
}
299

sebastien's avatar
sebastien committed
300
void
sebastien's avatar
sebastien committed
301
NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
302
303
                          const temporary_terms_t &temporary_terms,
                          deriv_node_temp_terms_t &tef_terms) const
304
{
305
  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
306
307
308
  if (it != temporary_terms.end())
    if (output_type == oMatlabDynamicModelSparse)
      output << "T" << idx << "(it_)";
ferhat's avatar
ferhat committed
309
    else
310
311
312
313
      output << "T" << idx;
  else
    output << datatree.num_constants.get(id);
}
sebastien's avatar
sebastien committed
314

sebastien's avatar
sebastien committed
315
double
316
NumConstNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
sebastien's avatar
sebastien committed
317
{
318
  return (datatree.num_constants.getDouble(id));
sebastien's avatar
sebastien committed
319
320
}

321
void
322
323
324
325
NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number,
                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
                      deriv_node_temp_terms_t &tef_terms) const
326
327
{
  FLDC_ fldc(datatree.num_constants.getDouble(id));
328
  fldc.write(CompileCode, instruction_number);
329
}
330

sebastien's avatar
sebastien committed
331
void
sebastien's avatar
sebastien committed
332
333
334
NumConstNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const
{
}
ferhat's avatar
ferhat committed
335

336
337
pair<int, expr_t >
NumConstNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
338
{
339
  /* return the numercial constant */
340
  return (make_pair(0, datatree.AddNonNegativeConstant(datatree.num_constants.get(id))));
341
}
ferhat's avatar
ferhat committed
342

343
344
expr_t
NumConstNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
ferhat's avatar
ferhat committed
345
{
ferhat's avatar
ferhat committed
346
  return datatree.Zero;
ferhat's avatar
ferhat committed
347
348
}

349
expr_t
sebastien's avatar
sebastien committed
350
NumConstNode::toStatic(DataTree &static_datatree) const
351
{
352
  return static_datatree.AddNonNegativeConstant(datatree.num_constants.get(id));
353
}
sebastien's avatar
sebastien committed
354

355
356
357
expr_t
NumConstNode::cloneDynamic(DataTree &dynamic_datatree) const
{
358
  return dynamic_datatree.AddNonNegativeConstant(datatree.num_constants.get(id));
359
360
}

sebastien's avatar
sebastien committed
361
362
363
364
365
int
NumConstNode::maxEndoLead() const
{
  return 0;
}
ferhat's avatar
ferhat committed
366

sebastien's avatar
sebastien committed
367
368
369
370
371
372
int
NumConstNode::maxExoLead() const
{
  return 0;
}

373
374
375
376
377
378
379
380
381
382
383
384
int
NumConstNode::maxEndoLag() const
{
  return 0;
}

int
NumConstNode::maxExoLag() const
{
  return 0;
}

385
expr_t
sebastien's avatar
sebastien committed
386
387
388
389
390
NumConstNode::decreaseLeadsLags(int n) const
{
  return const_cast<NumConstNode *>(this);
}

391
expr_t
sebastien's avatar
sebastien committed
392
NumConstNode::decreaseLeadsLagsPredeterminedVariables() const
393
{
sebastien's avatar
sebastien committed
394
  return const_cast<NumConstNode *>(this);
395
396
}

397
expr_t
398
NumConstNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
sebastien's avatar
sebastien committed
399
400
401
402
{
  return const_cast<NumConstNode *>(this);
}

403
expr_t
404
405
406
407
408
NumConstNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
{
  return const_cast<NumConstNode *>(this);
}

409
expr_t
410
NumConstNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
sebastien's avatar
sebastien committed
411
412
413
414
{
  return const_cast<NumConstNode *>(this);
}

415
expr_t
sebastien's avatar
sebastien committed
416
NumConstNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
sebastien's avatar
sebastien committed
417
418
419
420
{
  return const_cast<NumConstNode *>(this);
}

421
expr_t
422
423
424
425
426
NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
{
  return const_cast<NumConstNode *>(this);
}

427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
bool
NumConstNode::isNumConstNodeEqualTo(double value) const
{
  if (datatree.num_constants.getDouble(id) == value)
    return true;
  else
    return false;
}

bool
NumConstNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
  return false;
}

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
expr_t
NumConstNode::replaceTrendVar() const
{
  return const_cast<NumConstNode *>(this);
}

expr_t
NumConstNode::detrend(int symb_id, expr_t trend) const
{
  return const_cast<NumConstNode *>(this);
}

expr_t
NumConstNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
{
  return const_cast<NumConstNode *>(this);
}

Houtan Bastani's avatar
Houtan Bastani committed
460
461
462
463
464
465
466
467
468
469
470
471
472
473
VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg) :
  ExprNode(datatree_arg),
  symb_id(symb_id_arg),
  type(datatree.symbol_table.getType(symb_id_arg)),
  lag(lag_arg)
{
  // Add myself to the variable map
  datatree.variable_node_map[make_pair(symb_id, lag)] = this;

  // It makes sense to allow a lead/lag on parameters: during steady state calibration, endogenous and parameters can be swapped
  assert(type != eExternalFunction
         && (lag == 0 || (type != eModelLocalVariable && type != eModFileLocalVariable)));
}

sebastien's avatar
sebastien committed
474
475
476
477
478
void
VariableNode::prepareForDerivation()
{
  if (preparedForDerivation)
    return;
479

sebastien's avatar
sebastien committed
480
  preparedForDerivation = true;
sebastien's avatar
sebastien committed
481

sebastien's avatar
sebastien committed
482
  // Fill in non_null_derivatives
ferhat's avatar
ferhat committed
483
  switch (type)
sebastien's avatar
sebastien committed
484
485
486
487
488
    {
    case eEndogenous:
    case eExogenous:
    case eExogenousDet:
    case eParameter:
489
    case eTrend:
sebastien's avatar
sebastien committed
490
      // For a variable or a parameter, the only non-null derivative is with respect to itself
491
      non_null_derivatives.insert(datatree.getDerivID(symb_id, lag));
sebastien's avatar
sebastien committed
492
      break;
sebastien's avatar
sebastien committed
493
    case eModelLocalVariable:
sebastien's avatar
sebastien committed
494
      datatree.local_variables_table[symb_id]->prepareForDerivation();
sebastien's avatar
sebastien committed
495
      // Non null derivatives are those of the value of the local parameter
sebastien's avatar
sebastien committed
496
497
498
499
      non_null_derivatives = datatree.local_variables_table[symb_id]->non_null_derivatives;
      break;
    case eModFileLocalVariable:
      // Such a variable is never derived
sebastien's avatar
sebastien committed
500
      break;
501
    case eExternalFunction:
sebastien's avatar
sebastien committed
502
      cerr << "VariableNode::prepareForDerivation: impossible case" << endl;
503
      exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
504
505
506
    }
}

507
expr_t
sebastien's avatar
sebastien committed
508
VariableNode::computeDerivative(int deriv_id)
sebastien's avatar
sebastien committed
509
{
ferhat's avatar
ferhat committed
510
  switch (type)
sebastien's avatar
sebastien committed
511
512
513
514
    {
    case eEndogenous:
    case eExogenous:
    case eExogenousDet:
sebastien's avatar
sebastien committed
515
    case eParameter:
516
    case eTrend:
sebastien's avatar
sebastien committed
517
      if (deriv_id == datatree.getDerivID(symb_id, lag))
sebastien's avatar
sebastien committed
518
519
520
        return datatree.One;
      else
        return datatree.Zero;
sebastien's avatar
sebastien committed
521
    case eModelLocalVariable:
sebastien's avatar
sebastien committed
522
      return datatree.local_variables_table[symb_id]->getDerivative(deriv_id);
sebastien's avatar
sebastien committed
523
524
    case eModFileLocalVariable:
      cerr << "ModFileLocalVariable is not derivable" << endl;
525
      exit(EXIT_FAILURE);
526
    case eExternalFunction:
527
      cerr << "Impossible case!" << endl;
528
      exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
529
    }
sebastien's avatar
sebastien committed
530
  // Suppress GCC warning
531
  exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
532
533
}

534
void
535
VariableNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
536
{
537
  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this));
538
539
540
541
542
  if (it != temporary_terms.end())
    temporary_terms_inuse.insert(idx);
  if (type == eModelLocalVariable)
    datatree.local_variables_table[symb_id]->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
}
543

sebastien's avatar
sebastien committed
544
void
sebastien's avatar
sebastien committed
545
VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
546
547
                          const temporary_terms_t &temporary_terms,
                          deriv_node_temp_terms_t &tef_terms) const
548
549
{
  // If node is a temporary term
550
  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this));
551
552
553
554
555
556
557
558
  if (it != temporary_terms.end())
    {
      if (output_type == oMatlabDynamicModelSparse)
        output << "T" << idx << "(it_)";
      else
        output << "T" << idx;
      return;
    }
sebastien's avatar
sebastien committed
559

560
561
562
563
564
565
  if (IS_LATEX(output_type))
    {
      if (output_type == oLatexDynamicSteadyStateOperator)
        output << "\\bar{";
      output << datatree.symbol_table.getTeXName(symb_id);
      if (output_type == oLatexDynamicModel
566
          && (type == eEndogenous || type == eExogenous || type == eExogenousDet || type == eModelLocalVariable || type == eTrend))
567
568
569
570
571
572
573
574
        {
          output << "_{t";
          if (lag != 0)
            {
              if (lag > 0)
                output << "+";
              output << lag;
            }
575
          output << "}";
576
577
578
579
580
        }
      else if (output_type == oLatexDynamicSteadyStateOperator)
        output << "}";
      return;
    }
sebastien's avatar
sebastien committed
581

582
583
584
585
586
  int i;
  int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
  switch (type)
    {
    case eParameter:
587
      if (output_type == oMatlabOutsideModel || output_type == oSteadyStateFile)
588
589
590
591
        output << "M_.params" << "(" << tsid + 1 << ")";
      else
        output << "params" << LEFT_ARRAY_SUBSCRIPT(output_type) << tsid + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type);
      break;
ferhat's avatar
ferhat committed
592

593
    case eModelLocalVariable:
Sébastien Villemot's avatar
Sébastien Villemot committed
594
      if (output_type == oMatlabDynamicModelSparse || output_type == oMatlabStaticModelSparse)
595
596
        {
          output << "(";
597
          datatree.local_variables_table[symb_id]->writeOutput(output, output_type, temporary_terms, tef_terms);
598
599
600
          output << ")";
        }
      else
601
602
603
        /* We append underscores to avoid name clashes with "g1" or "oo_" (see
           also ModelTree::writeModelLocalVariables) */
        output << datatree.symbol_table.getName(symb_id) << "__";
604
      break;
ferhat's avatar
ferhat committed
605

606
607
608
609
    case eModFileLocalVariable:
      output << datatree.symbol_table.getName(symb_id);
      break;

610
611
612
613
614
615
616
617
    case eEndogenous:
      switch (output_type)
        {
        case oMatlabDynamicModel:
        case oCDynamicModel:
          i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
          output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          break;
618
        case oCStaticModel:
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
        case oMatlabStaticModel:
        case oMatlabStaticModelSparse:
          i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
          output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          break;
        case oMatlabDynamicModelSparse:
          i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
          if (lag > 0)
            output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          else if (lag < 0)
            output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_" << lag << ", " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          else
            output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_, " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          break;
        case oMatlabOutsideModel:
          output << "oo_.steady_state(" << tsid + 1 << ")";
          break;
        case oMatlabDynamicSteadyStateOperator:
637
        case oMatlabDynamicSparseSteadyStateOperator:
638
          output << "steady_state(" << tsid + 1 << ")";
639
          break;
640
641
642
        case oCDynamicSteadyStateOperator:
          output << "steady_state[" << tsid << "]";
          break;
643
644
645
        case oSteadyStateFile:
          output << "ys_(" << tsid + 1 << ")";
          break;
646
        default:
647
648
          cerr << "VariableNode::writeOutput: should not reach this point" << endl;
          exit(EXIT_FAILURE);
649
650
        }
      break;
ferhat's avatar
ferhat committed
651

652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
    case eExogenous:
      i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
      switch (output_type)
        {
        case oMatlabDynamicModel:
        case oMatlabDynamicModelSparse:
          if (lag > 0)
            output <<  "x(it_+" << lag << ", " << i << ")";
          else if (lag < 0)
            output <<  "x(it_" << lag << ", " << i << ")";
          else
            output <<  "x(it_, " << i << ")";
          break;
        case oCDynamicModel:
          if (lag == 0)
            output <<  "x[it_+" << i << "*nb_row_x]";
          else if (lag > 0)
            output <<  "x[it_+" << lag << "+" << i << "*nb_row_x]";
          else
            output <<  "x[it_" << lag << "+" << i << "*nb_row_x]";
          break;
673
        case oCStaticModel:
674
675
676
677
678
679
680
681
682
683
684
        case oMatlabStaticModel:
        case oMatlabStaticModelSparse:
          output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          break;
        case oMatlabOutsideModel:
          assert(lag == 0);
          output <<  "oo_.exo_steady_state(" << i << ")";
          break;
        case oMatlabDynamicSteadyStateOperator:
          output <<  "oo_.exo_steady_state(" << i << ")";
          break;
685
686
687
        case oSteadyStateFile:
          output << "exo_(" << i << ")";
          break;
688
        default:
689
690
          cerr << "VariableNode::writeOutput: should not reach this point" << endl;
          exit(EXIT_FAILURE);
691
692
        }
      break;
ferhat's avatar
ferhat committed
693

694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
    case eExogenousDet:
      i = tsid + datatree.symbol_table.exo_nbr() + ARRAY_SUBSCRIPT_OFFSET(output_type);
      switch (output_type)
        {
        case oMatlabDynamicModel:
        case oMatlabDynamicModelSparse:
          if (lag > 0)
            output <<  "x(it_+" << lag << ", " << i << ")";
          else if (lag < 0)
            output <<  "x(it_" << lag << ", " << i << ")";
          else
            output <<  "x(it_, " << i << ")";
          break;
        case oCDynamicModel:
          if (lag == 0)
709
            output <<  "x[it_+" << i << "*nb_row_x]";
710
          else if (lag > 0)
711
            output <<  "x[it_+" << lag << "+" << i << "*nb_row_x]";
712
          else
713
            output <<  "x[it_" << lag << "+" << i << "*nb_row_x]";
714
          break;
715
        case oCStaticModel:
716
717
718
719
720
721
722
723
724
725
726
        case oMatlabStaticModel:
        case oMatlabStaticModelSparse:
          output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
          break;
        case oMatlabOutsideModel:
          assert(lag == 0);
          output <<  "oo_.exo_det_steady_state(" << tsid + 1 << ")";
          break;
        case oMatlabDynamicSteadyStateOperator:
          output <<  "oo_.exo_det_steady_state(" << tsid + 1 << ")";
          break;
727
728
729
        case oSteadyStateFile:
          output << "exo_(" << i << ")";
          break;
730
        default:
731
732
          cerr << "VariableNode::writeOutput: should not reach this point" << endl;
          exit(EXIT_FAILURE);
733
734
        }
      break;
ferhat's avatar
ferhat committed
735

736
    case eExternalFunction:
737
    case eTrend:
738
739
740
741
      cerr << "Impossible case" << endl;
      exit(EXIT_FAILURE);
    }
}
sebastien's avatar
sebastien committed
742

sebastien's avatar
sebastien committed
743
double
744
VariableNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
sebastien's avatar
sebastien committed
745
{
746
  eval_context_t::const_iterator it = eval_context.find(symb_id);
sebastien's avatar
sebastien committed
747
  if (it == eval_context.end())
sebastien's avatar
sebastien committed
748
    throw EvalException();
sebastien's avatar
sebastien committed
749

sebastien's avatar
sebastien committed
750
  return it->second;
sebastien's avatar
sebastien committed
751
752
}

753
void
754
755
756
757
VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
                      deriv_node_temp_terms_t &tef_terms) const
758
759
{
  if (type == eModelLocalVariable || type == eModFileLocalVariable)
760
    datatree.local_variables_table[symb_id]->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
761
762
763
764
765
766
767
768
769
770
771
772
  else
    {
      int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
      if (type == eExogenousDet)
        tsid += datatree.symbol_table.exo_nbr();
      if (!lhs_rhs)
        {
          if (dynamic)
            {
              if (steady_dynamic)  // steady state values in a dynamic model
                {
                  FLDVS_ fldvs(type, tsid);
773
                  fldvs.write(CompileCode, instruction_number);
774
775
776
777
778
779
                }
              else
                {
                  if (type == eParameter)
                    {
                      FLDV_ fldv(type, tsid);
780
                      fldv.write(CompileCode, instruction_number);
781
782
783
784
                    }
                  else
                    {
                      FLDV_ fldv(type, tsid, lag);
785
                      fldv.write(CompileCode, instruction_number);
786
787
788
789
790
791
                    }
                }
            }
          else
            {
              FLDSV_ fldsv(type, tsid);
792
              fldsv.write(CompileCode, instruction_number);
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
            }
        }
      else
        {
          if (dynamic)
            {
              if (steady_dynamic)  // steady state values in a dynamic model
                {
                  cerr << "Impossible case: steady_state in rhs of equation" << endl;
                  exit(EXIT_FAILURE);
                }
              else
                {
                  if (type == eParameter)
                    {
                      FSTPV_ fstpv(type, tsid);
809
                      fstpv.write(CompileCode, instruction_number);
810
811
812
813
                    }
                  else
                    {
                      FSTPV_ fstpv(type, tsid, lag);
814
                      fstpv.write(CompileCode, instruction_number);
815
816
817
818
819
820
                    }
                }
            }
          else
            {
              FSTPSV_ fstpsv(type, tsid);
821
              fstpsv.write(CompileCode, instruction_number);
822
823
824
825
            }
        }
    }
}
826

827
void
828
VariableNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
829
                                    temporary_terms_t &temporary_terms,
830
                                    map<expr_t, pair<int, int> > &first_occurence,
831
                                    int Curr_block,
832
                                    vector<vector<temporary_terms_t> > &v_temporary_terms,
833
                                    int equation) const
834
835
836
837
{
  if (type == eModelLocalVariable)
    datatree.local_variables_table[symb_id]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
}
838

sebastien's avatar
sebastien committed
839
void
sebastien's avatar
sebastien committed
840
841
842
843
844
845
846
VariableNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const
{
  if (type == type_arg)
    result.insert(make_pair(symb_id, lag));
  if (type == eModelLocalVariable)
    datatree.local_variables_table[symb_id]->collectVariables(type_arg, result);
}
ferhat's avatar
ferhat committed
847

848
849
pair<int, expr_t>
VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
850
{
851
852
853
854
855
856
857
858
859
860
861
862
  /* The equation has to be normalized with respect to the current endogenous variable ascribed to it. 
     The two input arguments are : 
        - The ID of the endogenous variable associated to the equation.
        - The list of operators and operands needed to normalize the equation*
        
     The pair returned by NormalizeEquation is composed of 
      - a flag indicating if the expression returned contains (flag = 1) or not (flag = 0) 
        the endogenous variable related to the equation.
        If the expression contains more than one occurence of the associated endogenous variable, 
        the flag is equal to 2. 
      - an expression equal to the RHS if flag = 0 and equal to NULL elsewhere
  */
863
864
865
  if (type == eEndogenous)
    {
      if (datatree.symbol_table.getTypeSpecificID(symb_id) == var_endo && lag == 0)
866
        /* the endogenous variable */
867
        return (make_pair(1, (expr_t) NULL));
868
869
870
871
872
873
874
875
876
877
878
      else
        return (make_pair(0, datatree.AddVariableInternal(symb_id, lag)));
    }
  else
    {
      if (type == eParameter)
        return (make_pair(0, datatree.AddVariableInternal(symb_id, 0)));
      else
        return (make_pair(0, datatree.AddVariableInternal(symb_id, lag)));
    }
}
ferhat's avatar
ferhat committed
879

880
881
expr_t
VariableNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
sebastien's avatar
sebastien committed
882
{
ferhat's avatar
ferhat committed
883
884
885
886
887
888
  switch (type)
    {
    case eEndogenous:
    case eExogenous:
    case eExogenousDet:
    case eParameter:
889
    case eTrend:
sebastien's avatar
sebastien committed
890
      if (deriv_id == datatree.getDerivID(symb_id, lag))
ferhat's avatar
ferhat committed
891
892
893
894
        return datatree.One;
      else
        {
          //if there is in the equation a recursive variable we could use a chaine rule derivation
895
          map<int, expr_t>::const_iterator it = recursive_variables.find(datatree.getDerivID(symb_id, lag));
896
          if (it != recursive_variables.end())
ferhat's avatar
ferhat committed
897
            {
898
              map<int, expr_t>::const_iterator it2 = derivatives.find(deriv_id);
899
900
901
902
              if (it2 != derivatives.end())
                return it2->second;
              else
                {
903
                  map<int, expr_t> recursive_vars2(recursive_variables);
904
                  recursive_vars2.erase(it->first);
905
                  //expr_t c = datatree.AddNonNegativeConstant("1");
906
                  expr_t d = datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id, recursive_vars2));
907
                  //d = datatree.AddTimes(c, d);
sebastien's avatar
sebastien committed
908
                  derivatives[deriv_id] = d;
909
910
                  return d;
                }
ferhat's avatar
ferhat committed
911
912
913
914
915
            }
          else
            return datatree.Zero;
        }
    case eModelLocalVariable:
sebastien's avatar
sebastien committed
916
      return datatree.local_variables_table[symb_id]->getChainRuleDerivative(deriv_id, recursive_variables);
ferhat's avatar
ferhat committed
917
918
919
    case eModFileLocalVariable:
      cerr << "ModFileLocalVariable is not derivable" << endl;
      exit(EXIT_FAILURE);
920
    case eExternalFunction:
ferhat's avatar
ferhat committed
921
922
923
924
925
      cerr << "Impossible case!" << endl;
      exit(EXIT_FAILURE);
    }
  // Suppress GCC warning
  exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
926
927
}

928
expr_t
ferhat's avatar
ferhat committed
929
VariableNode::toStatic(DataTree &static_datatree) const
930
{
931
  return static_datatree.AddVariable(symb_id);
932
}
ferhat's avatar
ferhat committed
933

934
935
936
937
938
939
expr_t
VariableNode::cloneDynamic(DataTree &dynamic_datatree) const
{
  return dynamic_datatree.AddVariable(symb_id, lag);
}

sebastien's avatar
sebastien committed
940
941
942
int
VariableNode::maxEndoLead() const
{
943
  switch (type)
sebastien's avatar
sebastien committed
944
945
946
947
948
949
950
951
952
953
    {
    case eEndogenous:
      return max(lag, 0);
    case eModelLocalVariable:
      return datatree.local_variables_table[symb_id]->maxEndoLead();
    default:
      return 0;
    }
}

sebastien's avatar
sebastien committed
954
955
956
int
VariableNode::maxExoLead() const
{
957
  switch (type)
sebastien's avatar
sebastien committed
958
959
960
961
962
963
964
965
966
967
    {
    case eExogenous:
      return max(lag, 0);
    case eModelLocalVariable:
      return datatree.local_variables_table[symb_id]->maxExoLead();
    default:
      return 0;
    }
}

968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
int
VariableNode::maxEndoLag() const
{
  switch (type)
    {
    case eEndogenous:
      return max(-lag, 0);
    case eModelLocalVariable:
      return datatree.local_variables_table[symb_id]->maxEndoLag();
    default:
      return 0;
    }
}

int
VariableNode::maxExoLag() const
{
  switch (type)
    {
    case eExogenous:
      return max(-lag, 0);
    case eModelLocalVariable:
      return datatree.local_variables_table[symb_id]->maxExoLag();
    default:
      return 0;
    }
}

996
expr_t
sebastien's avatar
sebastien committed
997
998
VariableNode::decreaseLeadsLags(int n) const
{
999
  switch (type)
sebastien's avatar
sebastien committed
1000
1001
1002
1003
    {
    case eEndogenous:
    case eExogenous:
    case eExogenousDet:
1004
    case eTrend:
sebastien's avatar
sebastien committed
1005
1006
1007
1008
1009
1010
1011
1012
      return datatree.AddVariable(symb_id, lag-n);
    case eModelLocalVariable:
      return datatree.local_variables_table[symb_id]->decreaseLeadsLags(n);
    default:
      return const_cast<VariableNode *>(this);
    }
}

1013
expr_t
sebastien's avatar
sebastien committed
1014
VariableNode::decreaseLeadsLagsPredeterminedVariables() const
1015
{
sebastien's avatar
sebastien committed
1016
  if (datatree.symbol_table.isPredetermined(symb_id))
1017
1018
1019
1020
1021
    return decreaseLeadsLags(1);
  else
    return const_cast<VariableNode *>(this);
}

1022
expr_t
1023
VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
sebastien's avatar
sebastien committed
1024
{
1025
  expr_t value;
1026
  switch (type)
sebastien's avatar
sebastien committed
1027
1028
1029
1030
1031
    {
    case eEndogenous:
      if (lag <= 1)
        return const_cast<VariableNode *>(this);
      else
1032
        return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
sebastien's avatar
sebastien committed
1033
1034
1035
1036
1037
    case eModelLocalVariable:
      value = datatree.local_variables_table[symb_id];
      if (value->maxEndoLead() <= 1)
        return const_cast<VariableNode *>(this);
      else
1038
        return value->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
sebastien's avatar
sebastien committed
1039
1040
1041
1042
1043
    default:
      return const_cast<VariableNode *>(this);
    }
}

1044
expr_t
sebastien's avatar