DynamicModel.cc 293 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2018 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
31
32
#include <boost/filesystem.hpp>

#include "DynamicModel.hh"
sebastien's avatar
sebastien committed
33

34
35
36
void
DynamicModel::copyHelper(const DynamicModel &m)
{
37
  auto f = [this](const ExprNode *e) { return e->clone(*this); };
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95

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

  auto convert_vector_tt = [this,f](vector<temporary_terms_t> vtt)
    {
      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)
    equation_type_and_normalized_equation.push_back(make_pair(it.first, f(it.second)));

  for (const auto &it : m.blocks_derivatives)
    {
      block_derivatives_equation_variable_laglead_nodeid_t v;
      for (const auto &it2 : it)
        v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
      blocks_derivatives.push_back(v);
    }

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

  auto convert_derivative_t = [this,f](derivative_t dt)
    {
      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));

  for (const auto &it : m.pac_expectation_info)
    pac_expectation_info.insert(dynamic_cast<const PacExpectationNode *>(f(it)));
}


sebastien's avatar
sebastien committed
96
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
97
                           NumericalConstants &num_constants_arg,
98
                           ExternalFunctionsTable &external_functions_table_arg,
Houtan Bastani's avatar
Houtan Bastani committed
99
100
                           TrendComponentModelTable &trend_component_model_table_arg,
                           VarModelTable &var_model_table_arg) :
101
  ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
102
103
  trend_component_model_table{trend_component_model_table_arg},
  var_model_table{var_model_table_arg}
sebastien's avatar
sebastien committed
104
105
106
{
}

107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
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},
  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
159
{
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
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
240
241
242
243
244
  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;
  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;

  pac_expectation_info.clear();

  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
245
246
}

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

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

sebastien's avatar
sebastien committed
273
void
274
DynamicModel::computeTemporaryTermsOrdered()
sebastien's avatar
sebastien committed
275
{
276
  map<expr_t, pair<int, int>> first_occurence;
277
  map<expr_t, int> reference_count;
sebastien's avatar
sebastien committed
278
  BinaryOpNode *eq_node;
279
280
  first_derivatives_t::const_iterator it;
  first_chain_rule_derivatives_t::const_iterator it_chr;
sebastien's avatar
sebastien committed
281
  ostringstream tmp_s;
282
283
  v_temporary_terms.clear();
  map_idx.clear();
sebastien's avatar
sebastien committed
284

285
  unsigned int nb_blocks = getNbBlocks();
286
  v_temporary_terms = vector<vector<temporary_terms_t>>(nb_blocks);
287
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
sebastien's avatar
sebastien committed
288
  temporary_terms.clear();
289

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

389
390
391
392
393
void
DynamicModel::computeTemporaryTermsMapping()
{
  // Add a mapping form node ID to temporary terms order
  int j = 0;
394
395
  for (auto temporary_term : temporary_terms)
    map_idx[temporary_term->idx] = j++;
396
397
}

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

415
  local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
416
  if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
417
    local_temporary_terms = temporary_terms;
418
419
420
421
422
423
424
425
426
427
428
429

  //----------------------------------------------------------------------
  //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();
430
      nze_exo_det = derivative_exo_det[block].size();
431
432
433
434
      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;
435
      deriv_node_temp_terms_t tef_terms;
436
      local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
437
      if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
438
        local_temporary_terms = temporary_terms;
439

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

511
      tmp1_output.str("");
512
      tmp1_output << packageDir(basename + ".block") << "/dynamic_" << block+1 << ".m";
513
      output.open(tmp1_output.str(), ios::out | ios::binary);
514
515
516
517
518
519
520
521
      output << "%\n";
      output << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare\n";
      output << "%\n";
      output << "% Warning : this file is generated automatically by Dynare\n";
      output << "%           from model file (.mod)\n\n";
      output << "%/\n";
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        {
522
          output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)\n";
523
524
        }
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
525
        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
526
      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
527
        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
528
      else
529
        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)\n";
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
      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)
        {
          output << "  if(jacobian_eval)\n";
555
556
          output << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
557
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
558
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
559
560
561
562
563
          output << "  end;\n";
        }
      else
        {
          output << "  if(jacobian_eval)\n";
564
565
          output << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
566
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
567
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
568
569
570
          output << "  else\n";
          if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
            {
571
572
573
              output << "    g1 = spalloc(" << block_mfs << "*Periods, "
                     << block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
                     << ", " << nze << "*Periods);\n";
574
            }
ferhat's avatar
ferhat committed
575
          else
576
577
578
579
580
581
            {
              output << "    g1 = spalloc(" << block_mfs
                     << ", " << block_mfs << ", " << nze << ");\n";
            }
          output << "  end;\n";
        }
582

583
584
585
586
      output << "  g2=0;g3=0;\n";
      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
587
588
          for (int it : v_temporary_terms_inuse[block])
            tmp_output << " T" << it;
589
590
591
592
          output << "  global" << tmp_output.str() << ";\n";
        }
      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        {
593
          temporary_terms_t tt2;
594
595
596
597
598
599
600
          tt2.clear();
          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;
601
                  for (auto it : v_temporary_terms[block][i])
602
603
                    {
                      output << "  ";
604
                      // In the following, "Static" is used to avoid getting the "(it_)" subscripting
605
                      it->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, local_temporary_terms, {});
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
                      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)
        output << "  residual=zeros(" << block_mfs << ",1);\n";
      else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        output << "  residual=zeros(" << block_mfs << ",y_kmin+periods);\n";
      if (simulation_type == EVALUATE_BACKWARD)
        output << "  for it_ = (y_kmin+periods):y_kmin+1\n";
      if (simulation_type == EVALUATE_FORWARD)
        output << "  for it_ = y_kmin+1:(y_kmin+periods)\n";

      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
635
      temporary_terms_idxs_t temporary_terms_idxs;
636
637
      for (unsigned int i = 0; i < block_size; i++)
        {
638
          temporary_terms_t tt2;
639
640
641
642
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
643
              for (auto it : v_temporary_terms[block][i])
644
                {
645
                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
646
                    it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
647

648
                  output << "  " <<  sps;
649
                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
650
                  output << " = ";
651
                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
652
                  // Insert current node into tt2
653
                  tt2.insert(it);
654
655
656
657
658
659
660
                  output << ";" << endl;
                }
            }

          int variable_ID = getBlockVariableID(block, i);
          int equation_ID = getBlockEquationID(block, i);
          EquationType equ_type = getBlockEquationType(block, i);
661
          string sModel = symbol_table.getName(symbol_table.getID(SymbolType::endogenous, variable_ID));
662
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
663
664
665
          lhs = eq_node->get_arg1();
          rhs = eq_node->get_arg2();
          tmp_output.str("");
666
          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
667
668
669
670
671
672
673
674
675
676
677
678
          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 << " = ";
679
                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
680
681
682
683
684
685
686
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  output << "%" << tmp_output.str();
                  output << " = ";
                  if (isBlockEquationRenormalized(block, i))
                    {
687
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
688
689
                      output << "\n    ";
                      tmp_output.str("");
690
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
691
692
                      lhs = eq_node->get_arg1();
                      rhs = eq_node->get_arg2();
693
                      lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
694
                      output << " = ";
695
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
696
697
698
699
                    }
                }
              else
                {
700
                  cerr << "Type mismatch for equation " << equation_ID+1  << "\n";
701
702
703
704
705
706
707
708
709
710
711
712
                  exit(EXIT_FAILURE);
                }
              output << ";\n";
              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
713
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl;
714
715
716
717
718
719
720
721
              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
722
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl;
723
              Ufoss << "    b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
724
              Uf[equation_ID] += Ufoss.str();
725
              Ufoss.str("");
726
727
728
729
730
731
              output << "    residual(" << i+1-block_recursive << ", it_) = (";
              goto end;
            default:
            end:
              output << tmp_output.str();
              output << ") - (";
732
              rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
733
              output << ");\n";
sebastien's avatar
sebastien committed
734
#ifdef CONDITION
735
736
              if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
                output << "  condition(" << i+1 << ")=0;\n";
sebastien's avatar
sebastien committed
737
#endif
738
739
740
741
            }
        }
      // The Jacobian if we have to solve the block
      if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
742
        output << "  " << sps << "% Jacobian  " << endl << "    if jacobian_eval" << endl;
743
744
745
746
      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
747
        else
748
          output << "    % Jacobian  " << endl << "    if jacobian_eval" << endl;
749
750
751
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
752
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
753
        {
754
755
756
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
757
758
          int eqr = getBlockEquationID(block, eq);
          int varr = getBlockVariableID(block, var);
759
          if (var != prev_var || lag != prev_lag)
760
            {
761
762
763
764
              prev_var = var;
              prev_lag = lag;
              count_col++;
            }
765

766
          expr_t id = it->second;
767

768
          output << "      g1(" << eq+1 << ", " << count_col << ") = ";
769
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
770
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
771
                 << "(" << lag
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
772
773
                 << ") " << varr+1 << ", " << var+1
                 << ", equation=" << eqr+1 << ", " << eq+1 << endl;
774
775
776
777
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
778
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
779
780
781
782
783
784
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
785
            {
786
787
788
              prev_var = var;
              prev_lag = lag;
              count_col++;
789
            }
790
791
          expr_t id = it->second;
          output << "      g1_x(" << eqr+1 << ", " << count_col << ") = ";
792
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
793
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::exogenous, var))
794
795
796
797
798
799
800
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
801
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_det_derivative.begin(); it != tmp_block_exo_det_derivative.end(); it++)
802
803
804
805
806
807
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
808
            {
809
810
811
              prev_var = var;
              prev_lag = lag;
              count_col++;
812
            }
813
814
          expr_t id = it->second;
          output << "      g1_xd(" << eqr+1 << ", " << count_col << ") = ";
815
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
816
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::exogenous, var))
817
818
819
820
821
822
823
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
824
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_other_endo_derivative.begin(); it != tmp_block_other_endo_derivative.end(); it++)
825
826
827
828
829
830
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
831
            {
832
833
834
              prev_var = var;
              prev_lag = lag;
              count_col++;
835
            }
836
837
          expr_t id = it->second;

838
          output << "      g1_o(" << eqr+1 << ", " << /*var+1+(lag+block_max_lag)*block_size*/ count_col << ") = ";
839
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
840
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, var))
841
842
843
844
845
846
847
848
849
850
851
852
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      output << "      varargout{1}=g1_x;\n";
      output << "      varargout{2}=g1_xd;\n";
      output << "      varargout{3}=g1_o;\n";

      switch (simulation_type)
        {
        case EVALUATE_FORWARD:
        case EVALUATE_BACKWARD:
853
854
855
856
857
858
859
860
          output << "    end;" << endl;
          output << "  end;" << endl;
          break;
        case SOLVE_BACKWARD_SIMPLE:
        case SOLVE_FORWARD_SIMPLE:
        case SOLVE_BACKWARD_COMPLETE:
        case SOLVE_FORWARD_COMPLETE:
          output << "  else" << endl;
861
          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
862
863
864
865
866
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
867
              expr_t id = it->second.second;
868
              int lag = it->second.first;
869
870
871
              if (lag == 0)
                {
                  output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
872
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
873
                  output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
874
875
876
877
878
                         << "(" << lag
                         << ") " << varr+1
                         << ", equation=" << eqr+1 << endl;
                }

879
880
881
882
883
            }
          output << "  end;\n";
          break;
        case SOLVE_TWO_BOUNDARIES_SIMPLE:
        case SOLVE_TWO_BOUNDARIES_COMPLETE:
884
          output << "    else" << endl;
885
          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
886
887
888
889
890
891
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
              ostringstream tmp_output;
892
              expr_t id = it->second.second;
893
              int lag = it->second.first;
894
              if (eq >= block_recursive && var >= block_recursive)
895
896
                {
                  if (lag == 0)
897
898
899
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+Per_K_)*y(it_, " << varr+1 << ")";
900
                  else if (lag == 1)
901
902
903
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+Per_y_)*y(it_+1, " << varr+1 << ")";
904
                  else if (lag > 0)
905
906
907
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")";
908
                  else
909
910
911
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")";
912
                  Uf[eqr] += Ufoss.str();
913
914
                  Ufoss.str("");

915
916
917
918
919
920
921
922
923
924
925
926
927
                  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();
928
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
929
                  output << ";";
930
                  output << " %2 variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
931
932
933
                         << "(" << lag << ") " << varr+1
                         << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
                }
934

sebastien's avatar
sebastien committed
935
#ifdef CONDITION
936
937
              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
              output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
sebastien's avatar
sebastien committed
938
#endif
939
940
941
942
            }
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i >= block_recursive)
943
                output << "  " << Uf[getBlockEquationID(block, i)] << ";\n";
sebastien's avatar
sebastien committed
944
#ifdef CONDITION
945
946
              output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
              output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
sebastien's avatar
sebastien committed
947
#endif
948
            }
sebastien's avatar
sebastien committed
949
#ifdef CONDITION
950
951
952
953
954
955
956
957
958
959
960
961
962
963
          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];
                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
                }
            }
          for (i = 0; i < ModelBlock->Block_List[block].Size; i++)
            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
sebastien's avatar
sebastien committed
964
#endif
965
966
          output << "    end;" << endl;
          output << "  end;" << endl;
967
968
969
970
          break;
        default:
          break;
        }
971
      output << "end" << endl;
972
973
974
      output.close();
    }
}
sebastien's avatar
sebastien committed
975
976

void
977
DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &map_idx) const
978
{
979

980
981
  ostringstream tmp_output;
  ofstream code_file;
982
  unsigned int instruction_number = 0;
983
984
  bool file_open = false;

985
986
987
  boost::filesystem::create_directories(basename + "/model/bytecode");

  string main_name = basename + "/model/bytecode/dynamic.cod";
988
  code_file.open(main_name, ios::out | ios::binary | ios::ate);
989
990
  if (!code_file.is_open())
    {
991
      cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
      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;

1005
  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());
1006
1007
1008
1009
  file_open = true;

  //Temporary variables declaration
  FDIMT_ fdimt(temporary_terms.size());
1010
1011
1012
  fdimt.write(code_file, instruction_number);

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

1014
  for (int i = 0; i < symbol_table.exo_det_nbr(); i++)
1015
    exo_det.push_back(i);
1016
  for (int i = 0; i < symbol_table.exo_nbr(); i++)
1017
    exo.push_back(i);
1018

1019
  map<pair< int, pair<int, int>>, expr_t> first_derivatives_reordered_endo;
1020
  map<pair< pair<int, SymbolType>, pair<int, int>>, expr_t>  first_derivatives_reordered_exo;
1021
  for (const auto & first_derivative : first_derivatives)
1022
    {
1023
1024
      int deriv_id = first_derivative.first.second;
      unsigned int eq = first_derivative.first.first;
1025
1026
1027
      int symb = getSymbIDByDerivID(deriv_id);
      unsigned int var = symbol_table.getTypeSpecificID(symb);
      int lag = getLagByDerivID(deriv_id);
1028
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
1029
        first_derivatives_reordered_endo[{ lag, make_pair(var, eq) }] = first_derivative.second;
1030
      else if (getTypeByDerivID(deriv_id) == SymbolType::exogenous || getTypeByDerivID(deriv_id) == SymbolType::exogenousDet)
1031
        first_derivatives_reordered_exo[{ { lag, getTypeByDerivID(deriv_id) }, { var, eq } }] = first_derivative.second;
1032
1033
1034
1035
    }
  int prev_var = -1;
  int prev_lag = -999999999;
  int count_col_endo = 0;
1036
  for (map<pair< int, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
1037
1038
1039
1040
       it != first_derivatives_reordered_endo.end(); it++)
    {
      int var = it->first.second.first;
      int lag = it->first.first;
1041
      if (prev_var != var || prev_lag != lag)
1042
1043
1044
1045
1046
1047
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
    }
1048
1049
  prev_var = -1;
  prev_lag = -999999999;
1050
  SymbolType prev_type{SymbolType::unusedEndogenous}; // Any non-exogenous type would do here
1051
  int count_col_exo = 0;
1052
  int count_col_det_exo = 0;
1053

1054
  for (const auto & it : first_derivatives_reordered_exo)
1055
    {
1056
1057
1058
      int var{it.first.second.first};
      int lag{it.first.first.first};
      SymbolType type{it.first.first.second};
1059
      if (prev_var != var || prev_lag != lag || prev_type != type)
1060
1061
1062
        {
          prev_var = var;
          prev_lag = lag;
1063
          prev_type = type;
1064
          if (type == SymbolType::exogenous)
1065
            count_col_exo++;
1066
          else if (type == SymbolType::exogenousDet)
1067
            count_col_det_exo++;
1068
1069
        }
    }
1070

1071
1072
1073
1074
1075
1076
1077
1078
  FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
                           simulation_type,
                           0,
                           symbol_table.endo_nbr(),
                           variable_reordered,
                           equation_reordered,
                           false,
                           symbol_table.endo_nbr(),
1079
1080
                           max_endo_lag,
                           max_endo_lead,
1081
                           u_count_int,
1082
                           count_col_endo,
1083
                           symbol_table.exo_det_nbr(),
1084
1085
                           count_col_det_exo,
                           symbol_table.exo_nbr(),
1086
                           count_col_exo,
1087
                           0,
1088
1089
1090
1091
                           0,
                           exo_det,
                           exo,
                           other_endo
1092
                           );
1093
  fbeginblock.write(code_file, instruction_number);
1094

1095
  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
1096

1097
  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, true, false);
1098
1099

  FENDEQU_ fendequ;
1100
  fendequ.write(code_file, instruction_number);
1101
1102
1103
1104
1105
1106
1107

  // 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;

1108
  vector<vector<pair<pair<int, int>, int >>> derivatives;
1109
1110
  derivatives.resize(symbol_table.endo_nbr());
  count_u = symbol_table.endo_nbr();
1111
  for (const auto & first_derivative : first_derivatives)
1112
    {
1113
      int deriv_id = first_derivative.first.second;
1114
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
1115
        {
1116
1117
          expr_t d1 = first_derivative.second;
          unsigned int eq = first_derivative.first.first;
1118
1119
1120
          int symb = getSymbIDByDerivID(deriv_id);
          unsigned int var = symbol_table.getTypeSpecificID(symb);
          int lag = getLagByDerivID(deriv_id);
1121
          FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
1122
          fnumexpr.write(code_file, instruction_number);
1123
1124
          if (!derivatives[eq].size())
            derivatives[eq].clear();
1125
          derivatives[eq].emplace_back(make_pair(var, lag), count_u);
1126
          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
1127
1128

          FSTPU_ fstpu(count_u);