DynamicModel.cc 305 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2019 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
#include <iostream>
sebastien's avatar
sebastien committed
21
#include <cmath>
22
#include <cstdlib>
23
#include <cassert>
24
25
#include <cstdio>
#include <cerrno>
26
#include <algorithm>
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
27
#include <iterator>
28
#include <numeric>
sebastien's avatar
sebastien committed
29

30
#include "DynamicModel.hh"
sebastien's avatar
sebastien committed
31

32
33
34
void
DynamicModel::copyHelper(const DynamicModel &m)
{
35
  auto f = [this](const ExprNode *e) { return e->clone(*this); };
36
37
38
39

  for (const auto &it : m.static_only_equations)
    static_only_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));

40
  auto convert_vector_tt = [f](vector<temporary_terms_t> vtt)
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
    {
      vector<temporary_terms_t> vtt2;
      for (const auto &tt : vtt)
        {
          temporary_terms_t tt2;
          for (const auto &it : tt)
            tt2.insert(f(it));
          vtt2.push_back(tt2);
        }
      return vtt2;
    };

  for (const auto &it : m.v_temporary_terms)
    v_temporary_terms.push_back(convert_vector_tt(it));

  for (const auto &it : m.first_chain_rule_derivatives)
    first_chain_rule_derivatives[it.first] = f(it.second);

  for (const auto &it : m.equation_type_and_normalized_equation)
60
    equation_type_and_normalized_equation.push_back({it.first, f(it.second)});
61
62
63
64
65

  for (const auto &it : m.blocks_derivatives)
    {
      block_derivatives_equation_variable_laglead_nodeid_t v;
      for (const auto &it2 : it)
66
        v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2)));
67
68
69
70
71
72
      blocks_derivatives.push_back(v);
    }

  for (const auto &it : m.dynamic_jacobian)
    dynamic_jacobian[it.first] = f(it.second);

73
  auto convert_derivative_t = [f](derivative_t dt)
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
    {
      derivative_t dt2;
      for (const auto &it : dt)
        dt2[it.first] = f(it.second);
      return dt2;
    };
  for (const auto &it : m.derivative_endo)
    derivative_endo.push_back(convert_derivative_t(it));
  for (const auto &it : m.derivative_other_endo)
    derivative_other_endo.push_back(convert_derivative_t(it));
  for (const auto &it : m.derivative_exo)
    derivative_exo.push_back(convert_derivative_t(it));
  for (const auto &it : m.derivative_exo_det)
    derivative_exo_det.push_back(convert_derivative_t(it));
}


sebastien's avatar
sebastien committed
91
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
92
                           NumericalConstants &num_constants_arg,
93
                           ExternalFunctionsTable &external_functions_table_arg,
Houtan Bastani's avatar
Houtan Bastani committed
94
95
                           TrendComponentModelTable &trend_component_model_table_arg,
                           VarModelTable &var_model_table_arg) :
96
  ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
97
98
  trend_component_model_table{trend_component_model_table_arg},
  var_model_table{var_model_table_arg}
sebastien's avatar
sebastien committed
99
100
101
{
}

102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
DynamicModel::DynamicModel(const DynamicModel &m) :
  ModelTree {m},
  trend_component_model_table {m.trend_component_model_table},
  var_model_table {m.var_model_table},
  static_only_equations_lineno {m.static_only_equations_lineno},
  static_only_equations_equation_tags {m.static_only_equations_equation_tags},
  deriv_id_table {m.deriv_id_table},
  inv_deriv_id_table {m.inv_deriv_id_table},
  dyn_jacobian_cols_table {m.dyn_jacobian_cols_table},
  max_lag {m.max_lag},
  max_lead {m.max_lead},
  max_endo_lag {m.max_endo_lag},
  max_endo_lead {m.max_endo_lead},
  max_exo_lag {m.max_exo_lag},
  max_exo_lead {m.max_exo_lead},
  max_exo_det_lag {m.max_exo_det_lag},
  max_exo_det_lead {m.max_exo_det_lead},
  max_lag_orig {m.max_lag_orig},
  max_lead_orig {m.max_lead_orig},
121
  max_lag_with_diffs_expanded_orig {m.max_lag_with_diffs_expanded_orig},
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
  max_endo_lag_orig {m.max_endo_lag_orig},
  max_endo_lead_orig {m.max_endo_lead_orig},
  max_exo_lag_orig {m.max_exo_lag_orig},
  max_exo_lead_orig {m.max_exo_lead_orig},
  max_exo_det_lag_orig {m.max_exo_det_lag_orig},
  max_exo_det_lead_orig {m.max_exo_det_lead_orig},
  xrefs {m.xrefs},
  xref_param  {m.xref_param},
  xref_endo {m.xref_endo},
  xref_exo {m.xref_exo},
  xref_exo_det {m.xref_exo_det},
  nonzero_hessian_eqs {m.nonzero_hessian_eqs},
  v_temporary_terms_inuse {m.v_temporary_terms_inuse},
  map_idx {m.map_idx},
  global_temporary_terms {m.global_temporary_terms},
  block_type_firstequation_size_mfs {m.block_type_firstequation_size_mfs},
  blocks_linear {m.blocks_linear},
  other_endo_block {m.other_endo_block},
  exo_block {m.exo_block},
  exo_det_block {m.exo_det_block},
  block_var_exo {m.block_var_exo},
  block_exo_index {m.block_exo_index},
  block_det_exo_index {m.block_det_exo_index},
  block_other_endo_index {m.block_other_endo_index},
  block_col_type {m.block_col_type},
  variable_block_lead_lag {m.variable_block_lead_lag},
  equation_block {m.equation_block},
  var_expectation_functions_to_write {m.var_expectation_functions_to_write},
  endo_max_leadlag_block {m.endo_max_leadlag_block},
  other_endo_max_leadlag_block {m.other_endo_max_leadlag_block},
  exo_max_leadlag_block {m.exo_max_leadlag_block},
  exo_det_max_leadlag_block {m.exo_det_max_leadlag_block},
  max_leadlag_block {m.max_leadlag_block}
sebastien's avatar
sebastien committed
155
{
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
  copyHelper(m);
}

DynamicModel &
DynamicModel::operator=(const DynamicModel &m)
{
  ModelTree::operator=(m);

  assert(&trend_component_model_table == &m.trend_component_model_table);
  assert(&var_model_table == &m.var_model_table);

  static_only_equations_lineno = m.static_only_equations_lineno;
  static_only_equations_equation_tags = m.static_only_equations_equation_tags;
  deriv_id_table = m.deriv_id_table;
  inv_deriv_id_table = m.inv_deriv_id_table;
  dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
  max_lag = m.max_lag;
  max_lead = m.max_lead;
  max_endo_lag = m.max_endo_lag;
  max_endo_lead = m.max_endo_lead;
  max_exo_lag = m.max_exo_lag;
  max_exo_lead = m.max_exo_lead;
  max_exo_det_lag = m.max_exo_det_lag;
  max_exo_det_lead = m.max_exo_det_lead;
  max_lag_orig = m.max_lag_orig;
  max_lead_orig = m.max_lead_orig;
182
  max_lag_with_diffs_expanded_orig = m.max_lag_with_diffs_expanded_orig;
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
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
231
232
233
234
235
236
237
238
239
  max_endo_lag_orig = m.max_endo_lag_orig;
  max_endo_lead_orig = m.max_endo_lead_orig;
  max_exo_lag_orig = m.max_exo_lag_orig;
  max_exo_lead_orig = m.max_exo_lead_orig;
  max_exo_det_lag_orig = m.max_exo_det_lag_orig;
  max_exo_det_lead_orig = m.max_exo_det_lead_orig;
  xrefs = m.xrefs;
  xref_param  = m.xref_param;
  xref_endo = m.xref_endo;
  xref_exo = m.xref_exo;
  xref_exo_det = m.xref_exo_det;
  nonzero_hessian_eqs = m.nonzero_hessian_eqs;

  v_temporary_terms.clear();

  v_temporary_terms_inuse = m.v_temporary_terms_inuse;

  first_chain_rule_derivatives.clear();

  map_idx = m.map_idx;
  global_temporary_terms = m.global_temporary_terms;

  equation_type_and_normalized_equation.clear();

  block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;

  blocks_derivatives.clear();
  dynamic_jacobian.clear();

  blocks_linear = m.blocks_linear;

  derivative_endo.clear();
  derivative_other_endo.clear();
  derivative_exo.clear();
  derivative_exo_det.clear();

  other_endo_block = m.other_endo_block;
  exo_block = m.exo_block;
  exo_det_block = m.exo_det_block;
  block_var_exo = m.block_var_exo;
  block_exo_index = m.block_exo_index;
  block_det_exo_index = m.block_det_exo_index;
  block_other_endo_index = m.block_other_endo_index;
  block_col_type = m.block_col_type;
  variable_block_lead_lag = m.variable_block_lead_lag;
  equation_block = m.equation_block;
  var_expectation_functions_to_write = m.var_expectation_functions_to_write;

  endo_max_leadlag_block = m.endo_max_leadlag_block;
  other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
  exo_max_leadlag_block = m.exo_max_leadlag_block;
  exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
  max_leadlag_block = m.max_leadlag_block;

  copyHelper(m);

  return *this;
sebastien's avatar
sebastien committed
240
241
}

sebastien's avatar
sebastien committed
242
void
243
DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
244
{
245
246
  auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
  if (it != derivatives[1].end())
247
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
248
249
250
  else
    {
      FLDZ_ fldz;
251
      fldz.write(code_file, instruction_number);
252
253
    }
}
254
255

void
256
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
257
{
258
  auto it = first_chain_rule_derivatives.find({ eqr, varr, lag });
259
  if (it != first_chain_rule_derivatives.end())
260
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
261
  else
262
263
    {
      FLDZ_ fldz;
264
      fldz.write(code_file, instruction_number);
265
    }
266
267
}

sebastien's avatar
sebastien committed
268
void
269
DynamicModel::computeTemporaryTermsOrdered()
sebastien's avatar
sebastien committed
270
{
271
  map<expr_t, pair<int, int>> first_occurence;
272
  map<expr_t, int> reference_count;
sebastien's avatar
sebastien committed
273
  BinaryOpNode *eq_node;
274
  first_chain_rule_derivatives_t::const_iterator it_chr;
sebastien's avatar
sebastien committed
275
  ostringstream tmp_s;
276
277
  v_temporary_terms.clear();
  map_idx.clear();
sebastien's avatar
sebastien committed
278

279
  unsigned int nb_blocks = getNbBlocks();
280
  v_temporary_terms = vector<vector<temporary_terms_t>>(nb_blocks);
281
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
sebastien's avatar
sebastien committed
282
  temporary_terms.clear();
283

284
  if (!global_temporary_terms)
285
286
    {
      for (unsigned int block = 0; block < nb_blocks; block++)
sebastien's avatar
sebastien committed
287
        {
288
289
290
291
292
          reference_count.clear();
          temporary_terms.clear();
          unsigned int block_size = getBlockSize(block);
          unsigned int block_nb_mfs = getBlockMfs(block);
          unsigned int block_nb_recursives = block_size - block_nb_mfs;
293
          v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
294
          for (unsigned int i = 0; i < block_size; i++)
sebastien's avatar
sebastien committed
295
            {
296
              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
297
                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
298
              else
sebastien's avatar
sebastien committed
299
                {
300
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
301
                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
sebastien's avatar
sebastien committed
302
303
                }
            }
304
          for (const auto &it : blocks_derivatives[block])
305
            {
306
              expr_t id = get<3>(it);
307
308
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
            }
309
310
311
312
313
          for (const auto &it : derivative_endo[block])
            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
          for (const auto &it : derivative_other_endo[block])
            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
          v_temporary_terms_inuse[block] = {};
sebastien's avatar
sebastien committed
314
315
        }
    }
316
  else
sebastien's avatar
sebastien committed
317
    {
318
      for (unsigned int block = 0; block < nb_blocks; block++)
sebastien's avatar
sebastien committed
319
        {
320
321
322
323
          // Compute the temporary terms reordered
          unsigned int block_size = getBlockSize(block);
          unsigned int block_nb_mfs = getBlockMfs(block);
          unsigned int block_nb_recursives = block_size - block_nb_mfs;
324
          v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
325
          for (unsigned int i = 0; i < block_size; i++)
326
            {
327
              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
328
                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
329
330
              else
                {
331
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
332
333
                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
                }
334
            }
335
          for (const auto &it : blocks_derivatives[block])
sebastien's avatar
sebastien committed
336
            {
337
              expr_t id = get<3>(it);
338
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
sebastien's avatar
sebastien committed
339
            }
340
341
342
343
          for (const auto &it : derivative_endo[block])
            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
          for (const auto &it : derivative_other_endo[block])
            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
344
        }
345
      for (unsigned int block = 0; block < nb_blocks; block++)
346
        {
347
348
349
350
351
352
          // Collect the temporary terms reordered
          unsigned int block_size = getBlockSize(block);
          unsigned int block_nb_mfs = getBlockMfs(block);
          unsigned int block_nb_recursives = block_size - block_nb_mfs;
          set<int> temporary_terms_in_use;
          for (unsigned int i = 0; i < block_size; i++)
sebastien's avatar
sebastien committed
353
            {
354
              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
355
                getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
356
              else
sebastien's avatar
sebastien committed
357
                {
358
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
359
                  eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
sebastien's avatar
sebastien committed
360
361
                }
            }
362
          for (const auto &it : blocks_derivatives[block])
363
            {
364
              expr_t id = get<3>(it);
365
366
              id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
            }
367
368
369
370
371
372
373
374
          for (const auto &it : derivative_endo[block])
            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          for (const auto &it : derivative_other_endo[block])
            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          for (const auto &it : derivative_exo[block])
            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          for (const auto &it : derivative_exo_det[block])
            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
375
          v_temporary_terms_inuse[block] = temporary_terms_in_use;
sebastien's avatar
sebastien committed
376
        }
377
      computeTemporaryTermsMapping();
sebastien's avatar
sebastien committed
378
379
380
    }
}

381
382
383
384
385
void
DynamicModel::computeTemporaryTermsMapping()
{
  // Add a mapping form node ID to temporary terms order
  int j = 0;
386
387
  for (auto temporary_term : temporary_terms)
    map_idx[temporary_term->idx] = j++;
388
389
}

sebastien's avatar
sebastien committed
390
void
391
DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
392
393
394
{
  string tmp_s, sps;
  ostringstream tmp_output, tmp1_output, global_output;
395
  expr_t lhs = nullptr, rhs = nullptr;
396
  BinaryOpNode *eq_node;
397
398
  ostringstream Ufoss;
  vector<string> Uf(symbol_table.endo_nbr(), "");
399
  map<expr_t, int> reference_count;
400
  temporary_terms_t local_temporary_terms;
401
  ofstream  output;
402
  int nze, nze_exo, nze_exo_det, nze_other_endo;
403
404
  vector<int> feedback_variables;
  ExprNodeOutputType local_output_type;
405
  Ufoss.str("");
sebastien's avatar
sebastien committed
406

407
  local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
408
  if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
409
    local_temporary_terms = temporary_terms;
410
411
412
413
414
415
416
417
418
419
420
421

  //----------------------------------------------------------------------
  //For each block
  for (unsigned int block = 0; block < getNbBlocks(); block++)
    {

      //recursive_variables.clear();
      feedback_variables.clear();
      //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
      nze = derivative_endo[block].size();
      nze_other_endo = derivative_other_endo[block].size();
      nze_exo = derivative_exo[block].size();
422
      nze_exo_det = derivative_exo_det[block].size();
423
424
425
426
      BlockSimulationType simulation_type = getBlockSimulationType(block);
      unsigned int block_size = getBlockSize(block);
      unsigned int block_mfs = getBlockMfs(block);
      unsigned int block_recursive = block_size - block_mfs;
427
      deriv_node_temp_terms_t tef_terms;
428
      local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
429
      if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
430
        local_temporary_terms = temporary_terms;
431

432
433
      int prev_lag;
      unsigned int prev_var, count_col, count_col_endo, count_col_exo, count_col_exo_det, count_col_other_endo;
434
435
436
      map<tuple<int, int, int>, expr_t> tmp_block_endo_derivative;
      for (const auto &it : blocks_derivatives[block])
        tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
437
438
439
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_endo = 0;
440
      for (const auto &it : tmp_block_endo_derivative)
441
        {
442
443
          int lag = get<0>(it.first);
          unsigned int var = get<1>(it.first);
444
445
446
447
448
449
450
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_endo++;
            }
        }
451
452
453
      map<tuple<int, int, int>, expr_t> tmp_block_exo_derivative;
      for (const auto &it : derivative_exo[block])
        tmp_block_exo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
454
455
456
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_exo = 0;
457
      for (const auto &it : tmp_block_exo_derivative)
458
        {
459
460
          int lag = get<0>(it.first);
          unsigned int var = get<1>(it.first);
461
462
463
464
465
466
467
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_exo++;
            }
        }
468
469
470
      map<tuple<int, int, int>, expr_t> tmp_block_exo_det_derivative;
      for (const auto &it : derivative_exo_det[block])
        tmp_block_exo_det_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
471
472
473
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_exo_det = 0;
474
      for (const auto &it : tmp_block_exo_det_derivative)
475
        {
476
477
          int lag = get<0>(it.first);
          unsigned int var = get<1>(it.first);
478
479
480
481
482
483
484
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_exo_det++;
            }
        }
485
486
487
      map<tuple<int, int, int>, expr_t> tmp_block_other_endo_derivative;
      for (const auto &it : derivative_other_endo[block])
        tmp_block_other_endo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
488
489
490
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_other_endo = 0;
491
      for (const auto &it : tmp_block_other_endo_derivative)
492
        {
493
494
          int lag = get<0>(it.first);
          unsigned int var = get<1>(it.first);
495
496
497
498
499
500
501
502
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_other_endo++;
            }
        }

503
      tmp1_output.str("");
504
      tmp1_output << packageDir(basename + ".block") << "/dynamic_" << block+1 << ".m";
505
      output.open(tmp1_output.str(), ios::out | ios::binary);
506
507
508
509
510
511
      output << "%" << endl
             << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare" << endl
             << "%" << endl
             << "% Warning : this file is generated automatically by Dynare" << endl
             << "%           from model file (.mod)" << endl << endl
             << "%/" << endl;
512
513
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        {
514
          output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl;
515
516
        }
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
517
        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
518
      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
519
        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
520
      else
521
        output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl;
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
      BlockType block_type;
      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        block_type = SIMULTAN;
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
        block_type = SIMULTANS;
      else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
                || simulation_type == EVALUATE_BACKWARD    || simulation_type == EVALUATE_FORWARD)
               && getBlockFirstEquation(block) < prologue)
        block_type = PROLOGUE;
      else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
                || simulation_type == EVALUATE_BACKWARD    || simulation_type == EVALUATE_FORWARD)
               && getBlockFirstEquation(block) >= equations.size() - epilogue)
        block_type = EPILOGUE;
      else
        block_type = SIMULTANS;
      output << "  % ////////////////////////////////////////////////////////////////////////" << endl
             << "  % //" << string("                     Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
             << "          //" << endl
             << "  % //                     Simulation type "
             << BlockSim(simulation_type) << "  //" << endl
             << "  % ////////////////////////////////////////////////////////////////////////" << endl;
      //The Temporary terms
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        {
546
547
548
549
550
551
          output << "  if(jacobian_eval)" << endl
                 << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");" << endl
                 << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");" << endl
                 << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");" << endl
                 << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");" << endl
                 << "  end;" << endl;
552
553
554
        }
      else
        {
555
556
557
558
559
560
          output << "  if(jacobian_eval)" << endl
                 << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");" << endl
                 << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");" << endl
                 << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");" << endl
                 << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");" << endl
                 << "  else" << endl;
561
          if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
562
563
564
            output << "    g1 = spalloc(" << block_mfs << "*Periods, "
                   << block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
                   << ", " << nze << "*Periods);" << endl;
ferhat's avatar
ferhat committed
565
          else
566
567
568
            output << "    g1 = spalloc(" << block_mfs
                   << ", " << block_mfs << ", " << nze << ");" << endl;
          output << "  end;" << endl;
569
        }
570

571
      output << "  g2=0;g3=0;" << endl;
572
573
574
      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
575
576
          for (int it : v_temporary_terms_inuse[block])
            tmp_output << " T" << it;
577
          output << "  global" << tmp_output.str() << ";" << endl;
578
579
580
        }
      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        {
581
          temporary_terms_t tt2;
582
583
584
585
586
587
          for (int i = 0; i < (int) block_size; i++)
            {
              if (v_temporary_terms[block][i].size() && global_temporary_terms)
                {
                  output << "  " << "% //Temporary variables initialization" << endl
                         << "  " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
588
                  for (auto it : v_temporary_terms[block][i])
589
590
                    {
                      output << "  ";
591
                      // In the following, "Static" is used to avoid getting the "(it_)" subscripting
592
                      it->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, local_temporary_terms, {});
593
594
595
596
597
598
                      output << " = T_zeros;" << endl;
                    }
                }
            }
        }
      if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
599
        output << "  residual=zeros(" << block_mfs << ",1);" << endl;
600
      else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
601
        output << "  residual=zeros(" << block_mfs << ",y_kmin+periods);" << endl;
602
      if (simulation_type == EVALUATE_BACKWARD)
603
        output << "  for it_ = (y_kmin+periods):y_kmin+1" << endl;
604
      if (simulation_type == EVALUATE_FORWARD)
605
        output << "  for it_ = y_kmin+1:(y_kmin+periods)" << endl;
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621

      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        {
          output << "  b = zeros(periods*y_size,1);" << endl
                 << "  for it_ = y_kmin+1:(periods+y_kmin)" << endl
                 << "    Per_y_=it_*y_size;" << endl
                 << "    Per_J_=(it_-y_kmin-1)*y_size;" << endl
                 << "    Per_K_=(it_-1)*y_size;" << endl;
          sps = "  ";
        }
      else
        if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
          sps = "  ";
        else
          sps = "";
      // The equations
Houtan Bastani's avatar
Houtan Bastani committed
622
      temporary_terms_idxs_t temporary_terms_idxs;
623
624
      for (unsigned int i = 0; i < block_size; i++)
        {
625
          temporary_terms_t tt2;
626
627
628
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
629
              for (auto it : v_temporary_terms[block][i])
630
                {
631
                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
632
                    it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
633

634
                  output << "  " <<  sps;
635
                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
636
                  output << " = ";
637
                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
638
                  // Insert current node into tt2
639
                  tt2.insert(it);
640
641
642
643
644
645
646
                  output << ";" << endl;
                }
            }

          int variable_ID = getBlockVariableID(block, i);
          int equation_ID = getBlockEquationID(block, i);
          EquationType equ_type = getBlockEquationType(block, i);
647
          string sModel = symbol_table.getName(symbol_table.getID(SymbolType::endogenous, variable_ID));
648
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
649
650
          lhs = eq_node->arg1;
          rhs = eq_node->arg2;
651
          tmp_output.str("");
652
          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
653
654
655
656
657
658
659
660
661
662
663
664
          switch (simulation_type)
            {
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
            evaluation:     if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
                output << "    % equation " << getBlockEquationID(block, i)+1 << " variable : " << sModel
                       << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
              output << "    ";
              if (equ_type == E_EVALUATE)
                {
                  output << tmp_output.str();
                  output << " = ";
665
                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
666
667
668
669
670
671
672
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  output << "%" << tmp_output.str();
                  output << " = ";
                  if (isBlockEquationRenormalized(block, i))
                    {
673
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
674
                      output << endl << "    ";
675
                      tmp_output.str("");
676
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
677
678
                      lhs = eq_node->arg1;
                      rhs = eq_node->arg2;
679
                      lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
680
                      output << " = ";
681
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
682
683
684
685
                    }
                }
              else
                {
686
                  cerr << "Type mismatch for equation " << equation_ID+1  << endl;
687
688
                  exit(EXIT_FAILURE);
                }
689
              output << ";" << endl;
690
691
692
693
694
695
696
697
698
              break;
            case SOLVE_BACKWARD_SIMPLE:
            case SOLVE_FORWARD_SIMPLE:
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              if (i < block_recursive)
                goto evaluation;
              feedback_variables.push_back(variable_ID);
              output << "  % equation " << equation_ID+1 << " variable : " << sModel
699
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl;
700
701
702
703
704
705
706
707
              output << "  " << "residual(" << i+1-block_recursive << ") = (";
              goto end;
            case SOLVE_TWO_BOUNDARIES_COMPLETE:
            case SOLVE_TWO_BOUNDARIES_SIMPLE:
              if (i < block_recursive)
                goto evaluation;
              feedback_variables.push_back(variable_ID);
              output << "    % equation " << equation_ID+1 << " variable : " << sModel
708
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl;
709
              Ufoss << "    b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
710
              Uf[equation_ID] += Ufoss.str();
711
              Ufoss.str("");
712
713
714
715
716
717
              output << "    residual(" << i+1-block_recursive << ", it_) = (";
              goto end;
            default:
            end:
              output << tmp_output.str();
              output << ") - (";
718
              rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
719
              output << ");" << endl;
sebastien's avatar
sebastien committed
720
#ifdef CONDITION
721
              if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
722
                output << "  condition(" << i+1 << ")=0;" << endl;
sebastien's avatar
sebastien committed
723
#endif
724
725
726
727
            }
        }
      // The Jacobian if we have to solve the block
      if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
728
        output << "  " << sps << "% Jacobian  " << endl << "    if jacobian_eval" << endl;
729
730
731
732
      else
        if (simulation_type == SOLVE_BACKWARD_SIMPLE   || simulation_type == SOLVE_FORWARD_SIMPLE
            || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
          output << "  % Jacobian  " << endl << "  if jacobian_eval" << endl;
sebastien's avatar
sebastien committed
733
        else
734
          output << "    % Jacobian  " << endl << "    if jacobian_eval" << endl;
735
736
737
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
738
      for (const auto &it : tmp_block_endo_derivative)
739
        {
740
741
742
          int lag;
          unsigned int var, eq;
          tie(lag, var, eq) = it.first;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
743
744
          int eqr = getBlockEquationID(block, eq);
          int varr = getBlockVariableID(block, var);
745
          if (var != prev_var || lag != prev_lag)
746
            {
747
748
749
750
              prev_var = var;
              prev_lag = lag;
              count_col++;
            }
751

752
          expr_t id = it.second;
753

754
          output << "      g1(" << eq+1 << ", " << count_col << ") = ";
755
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
756
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
757
                 << "(" << lag
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
758
759
                 << ") " << varr+1 << ", " << var+1
                 << ", equation=" << eqr+1 << ", " << eq+1 << endl;
760
761
762
763
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
764
      for (const auto &it : tmp_block_exo_derivative)
765
        {
766
767
768
          int lag;
          unsigned int var, eq;
          tie(lag, var, eq) = it.first;
769
770
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
771
            {
772
773
774
              prev_var = var;
              prev_lag = lag;
              count_col++;
775
            }
776
          expr_t id = it.second;
777
          output << "      g1_x(" << eqr+1 << ", " << count_col << ") = ";
778
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
779
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::exogenous, var))
780
781
782
783
784
785
786
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
787
      for (const auto &it : tmp_block_exo_det_derivative)
788
        {
789
790
791
          int lag;
          unsigned int var, eq;
          tie(lag, var, eq) = it.first;
792
793
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
794
            {
795
796
797
              prev_var = var;
              prev_lag = lag;
              count_col++;
798
            }
799
          expr_t id = it.second;
800
          output << "      g1_xd(" << eqr+1 << ", " << count_col << ") = ";
801
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
802
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::exogenous, var))
803
804
805
806
807
808
809
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
810
      for (const auto &it : tmp_block_other_endo_derivative)
811
        {
812
813
814
          int lag;
          unsigned int var, eq;
          tie(lag, var, eq) = it.first;
815
816
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
817
            {
818
819
820
              prev_var = var;
              prev_lag = lag;
              count_col++;
821
            }
822
          expr_t id = it.second;
823

824
          output << "      g1_o(" << eqr+1 << ", " << /*var+1+(lag+block_max_lag)*block_size*/ count_col << ") = ";
825
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
826
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, var))
827
828
829
830
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
831
832
833
      output << "      varargout{1}=g1_x;" << endl
             << "      varargout{2}=g1_xd;" << endl
             << "      varargout{3}=g1_o;" << endl;
834
835
836
837
838

      switch (simulation_type)
        {
        case EVALUATE_FORWARD:
        case EVALUATE_BACKWARD:
839
840
          output << "    end;" << endl
                 << "  end;" << endl;
841
842
843
844
845
846
          break;
        case SOLVE_BACKWARD_SIMPLE:
        case SOLVE_FORWARD_SIMPLE:
        case SOLVE_BACKWARD_COMPLETE:
        case SOLVE_FORWARD_COMPLETE:
          output << "  else" << endl;
847
          for (const auto &it : blocks_derivatives[block])
848
            {
849
850
851
852
              unsigned int eq, var;
              expr_t id;
              int lag;
              tie(eq, var, lag, id) = it;
853
854
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
855
856
857
              if (lag == 0)
                {
                  output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
858
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
859
                  output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
860
861
862
863
864
                         << "(" << lag
                         << ") " << varr+1
                         << ", equation=" << eqr+1 << endl;
                }

865
            }
866
          output << "  end;" << endl;
867
868
869
          break;
        case SOLVE_TWO_BOUNDARIES_SIMPLE:
        case SOLVE_TWO_BOUNDARIES_COMPLETE:
870
          output << "    else" << endl;
871
          for (const auto &it : blocks_derivatives[block])
872
            {
873
874
875
876
              unsigned int eq, var;
              int lag;
              expr_t id;
              tie(eq, var, lag, id) = it;
877
878
879
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
              ostringstream tmp_output;
880
              if (eq >= block_recursive && var >= block_recursive)
881
882
                {
                  if (lag == 0)
883
884
885
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+Per_K_)*y(it_, " << varr+1 << ")";
886
                  else if (lag == 1)
887
888
889
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+Per_y_)*y(it_+1, " << varr+1 << ")";
890
                  else if (lag > 0)
891
892
893
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")";
894
                  else
895
896
897
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")";
898
                  Uf[eqr] += Ufoss.str();
899
900
                  Ufoss.str("");

901
902
903
904
905
906
907
908
909
910
911
912
913
                  if (lag == 0)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+Per_K_) = ";
                  else if (lag == 1)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+Per_y_) = ";
                  else if (lag > 0)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+y_size*(it_+" << lag-1 << ")) = ";
                  else if (lag < 0)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+y_size*(it_" << lag-1 << ")) = ";
                  output << " " << tmp_output.str();
914
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
915
                  output << ";";
916
                  output << " %2 variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
917
918
919
                         << "(" << lag << ") " << varr+1
                         << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
                }
920

sebastien's avatar
sebastien committed
921
#ifdef CONDITION
922
923
              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))" << endl
                     << "    condition(" << eqr << ")=u(" << u << "+Per_u_);" << endl;
sebastien's avatar
sebastien committed
924
#endif
925
926
927
928
            }
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i >= block_recursive)
929
                output << "  " << Uf[getBlockEquationID(block, i)] << ";" << endl;
sebastien's avatar
sebastien committed
930
#ifdef CONDITION
931
932
              output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))" << endl
                     << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);" << endl;
sebastien's avatar
sebastien committed
933
#endif
934
            }
sebastien's avatar
sebastien committed
935
#ifdef CONDITION
936
937
938
939
940
941
942
943
944
          for (m = 0; m <= ModelBlock->Block_List[block].Max_Lead+ModelBlock->Block_List[block].Max_Lag; m++)
            {
              k = m-ModelBlock->Block_List[block].Max_Lag;
              for (i = 0; i < ModelBlock->Block_List[block].IM_lead_lag[m].size; i++)
                {
                  unsigned int eq = ModelBlock->Block_List[block].IM_lead_lag[m].Equ_Index[i];
                  unsigned int var = ModelBlock->Block_List[block].IM_lead_lag[m].Var_Index[i];
                  unsigned int u = ModelBlock->Block_List[block].IM_lead_lag[m].u[i];
                  unsigned int eqr = ModelBlock->Block_List[block].IM_lead_lag[m].Equ[i];
945
                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");" << endl;
946
947
948
                }
            }
          for (i = 0; i < ModelBlock->Block_List[block].Size; i++)
949
            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");" << endl;
sebastien's avatar
sebastien committed
950
#endif
951
952
          output << "    end;" << endl
                 << "  end;" << endl;
953
954
955
956
          break;
        default:
          break;
        }
957
      output << "end" << endl;
958
959
960
      output.close();
    }
}
sebastien's avatar
sebastien committed
961
962

void
963
DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &map_idx) const
964
{
965

966
967
  ostringstream tmp_output;
  ofstream code_file;
968
  unsigned int instruction_number = 0;
969
970
  bool file_open = false;

971
972
973
  boost::filesystem::create_directories(basename + "/model/bytecode");

  string main_name = basename + "/model/bytecode/dynamic.cod";
974
  code_file.open(main_name, ios::out | ios::binary | ios::ate);
975
976
  if (!code_file.is_open())
    {
977
      cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
978
979
980
981
982
983
984
985
986
987
988
989
990
      exit(EXIT_FAILURE);
    }

  int count_u;
  int u_count_int = 0;
  BlockSimulationType simulation_type;
  if ((max_endo_lag > 0) && (max_endo_lead > 0))
    simulation_type = SOLVE_TWO_BOUNDARIES_COMPLETE;
  else if ((max_endo_lag >= 0) && (max_endo_lead == 0))
    simulation_type = SOLVE_FORWARD_COMPLETE;
  else
    simulation_type = SOLVE_BACKWARD_COMPLETE;

991
  Write_Inf_To_Bin_File(basename + "/model/bytecode/dynamic.bin", u_count_int, file_open, simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE, symbol_table.endo_nbr());
992
993
994
995
  file_open = true;

  //Temporary variables declaration
  FDIMT_ fdimt(temporary_terms.size());
996
997
998
  fdimt.write(code_file, instruction_number);

  vector<unsigned int> exo, exo_det, other_endo;
999

1000
  for (int i = 0; i < symbol_table.exo_det_nbr(); i++)
1001
    exo_det.push_back(i);
1002
  for (int i = 0; i < symbol_table.exo_nbr(); i++)
1003
    exo.push_back(i);
1004

1005
1006
  map<tuple<int, int, int>, expr_t> first_derivatives_reordered_endo;
  map<tuple<int, SymbolType, int, int>, expr_t> first_derivatives_reordered_exo;
1007
  for (const auto & first_derivative : derivatives[1])
1008
    {
1009
1010
      int deriv_id = first_derivative.first[1];
      unsigned int eq = first_derivative.first[0];
1011
1012
1013
      int symb = getSymbIDByDerivID(deriv_id);
      unsigned int var = symbol_table.getTypeSpecificID(symb);
      int lag = getLagByDerivID(deriv_id);
1014
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
1015
        first_derivatives_reordered_endo[{ lag, var, eq }] = first_derivative.second;
1016
      else if (getTypeByDerivID(deriv_id) == SymbolType::exogenous || getTypeByDerivID(deriv_id) == SymbolType::exogenousDet)
1017
        first_derivatives_reordered_exo[{ lag, getTypeByDerivID(deriv_id), var, eq }] = first_derivative.second;
1018
1019
1020
1021
    }
  int prev_var = -1;
  int prev_lag = -999999999;
  int count_col_endo = 0;
1022
  for (const auto &it : first_derivatives_reordered_endo)
1023
    {
1024
1025
      int var, lag;
      tie(lag, var, ignore) = it.first;
1026
      if (prev_var != var || prev_lag != lag)
1027
1028
1029
1030
1031
1032
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
    }
1033
1034
  prev_var = -1;
  prev_lag = -999999999;
1035
  SymbolType prev_type{SymbolType::unusedEndogenous}; // Any non-exogenous type would do here
1036
  int count_col_exo = 0;
1037
  int count_col_det_exo = 0;
1038

1039
  for (const auto & it : first_derivatives_reordered_exo)
1040
    {
1041
1042
1043
      int var, lag;
      SymbolType type;
      tie(lag, type, var, ignore) = it.first;
1044
      if (prev_var != var || prev_lag != lag || prev_type != type)
1045
1046
1047
        {
          prev_var = var;
          prev_lag = lag;
1048
          prev_type = type;
1049
          if (type == SymbolType::exogenous)
1050
            count_col_exo++;
1051
          else if (type == SymbolType::exogenousDet)
1052
            count_col_det_exo++;
1053
1054
        }
    }
1055

1056
1057
1058
1059
1060
1061
1062
1063
  FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
                           simulation_type,
                           0,
                           symbol_table.endo_nbr(),
                           variable_reordered,
                           equation_reordered,
                           false,
                           symbol_table.endo_nbr(),
1064
1065
                           max_endo_lag,
                           max_endo_lead,
1066
                           u_count_int,
1067
                           count_col_endo,
1068
                           symbol_table.exo_det_nbr(),
1069
1070
                           count_col_det_exo,
                           symbol_table.exo_nbr(),
1071
                           count_col_exo,
1072
                           0,
1073
1074
1075
1076
                           0,
                           exo_det,
                           exo,
                           other_endo
1077
                           );
1078
  fbeginblock.write(code_file, instruction_number);
1079

1080
  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
1081

1082
  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, true, false);
1083
1084

  FENDEQU_ fendequ;
1085
  fendequ.write(code_file, instruction_number);
1086
1087
1088
1089
1090
1091
1092

  // Get the current code_file position and jump if eval = true
  streampos pos1 = code_file.tellp();
  FJMPIFEVAL_ fjmp_if_eval(0);
  fjmp_if_eval.write(code_file, instruction_number);
  int prev_instruction_number = instruction_number;

1093
  vector<vector<tuple<int, int, int>>> my_derivatives(symbol_table.endo_nbr());;
1094
  count_u = symbol_table.endo_nbr();
1095
  for (const auto & first_derivative : derivatives[1])
1096
    {