DynamicModel.cc 289 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
#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
60
61
62
63
64
65
    {
      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)
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
91
92
93
    {
      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
94
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
95
                           NumericalConstants &num_constants_arg,
96
                           ExternalFunctionsTable &external_functions_table_arg,
Houtan Bastani's avatar
Houtan Bastani committed
97
98
                           TrendComponentModelTable &trend_component_model_table_arg,
                           VarModelTable &var_model_table_arg) :
99
  ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
100
101
  trend_component_model_table{trend_component_model_table_arg},
  var_model_table{var_model_table_arg}
sebastien's avatar
sebastien committed
102
103
104
{
}

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
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
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
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
  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
243
244
}

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

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

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

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

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

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

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

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

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

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

506
      tmp1_output.str("");
507
      tmp1_output << packageDir(basename + ".block") << "/dynamic_" << block+1 << ".m";
508
      output.open(tmp1_output.str(), ios::out | ios::binary);
509
510
511
512
513
514
      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;
515
516
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        {
517
          output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl;
518
519
        }
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
520
        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
521
      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
522
        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
523
      else
524
        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;
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
      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)
        {
549
550
551
552
553
554
          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;
555
556
557
        }
      else
        {
558
559
560
561
562
563
          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;
564
          if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
565
566
567
            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
568
          else
569
570
571
            output << "    g1 = spalloc(" << block_mfs
                   << ", " << block_mfs << ", " << nze << ");" << endl;
          output << "  end;" << endl;
572
        }
573

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

      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
625
      temporary_terms_idxs_t temporary_terms_idxs;
626
627
      for (unsigned int i = 0; i < block_size; i++)
        {
628
          temporary_terms_t tt2;
629
630
631
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
632
              for (auto it : v_temporary_terms[block][i])
633
                {
634
                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
635
                    it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
636

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

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

755
          expr_t id = it.second;
756

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

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

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

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

904
905
906
907
908
909
910
911
912
913
914
915
916
                  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();
917
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
918
                  output << ";";
919
                  output << " %2 variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
920
921
922
                         << "(" << lag << ") " << varr+1
                         << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
                }
923

sebastien's avatar
sebastien committed
924
#ifdef CONDITION
925
926
              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))" << endl
                     << "    condition(" << eqr << ")=u(" << u << "+Per_u_);" << endl;
sebastien's avatar
sebastien committed
927
#endif
928
929
930
931
            }
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i >= block_recursive)
932
                output << "  " << Uf[getBlockEquationID(block, i)] << ";" << endl;
sebastien's avatar
sebastien committed
933
#ifdef CONDITION
934
935
              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
936
#endif
937
            }
sebastien's avatar
sebastien committed
938
#ifdef CONDITION
939
940
941
942
943
944
945
946
947
          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];
948
                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");" << endl;
949
950
951
                }
            }
          for (i = 0; i < ModelBlock->Block_List[block].Size; i++)
952
            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");" << endl;
sebastien's avatar
sebastien committed
953
#endif
954
955
          output << "    end;" << endl
                 << "  end;" << endl;
956
957
958
959
<