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
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},
124
  max_lag_with_diffs_expanded_orig {m.max_lag_with_diffs_expanded_orig},
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
  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
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
  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;
185
  max_lag_with_diffs_expanded_orig = m.max_lag_with_diffs_expanded_orig;
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
  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
251
  auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
  if (it != derivatives[1].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
  first_chain_rule_derivatives_t::const_iterator it_chr;
sebastien's avatar
sebastien committed
280
  ostringstream tmp_s;
281
282
  v_temporary_terms.clear();
  map_idx.clear();
sebastien's avatar
sebastien committed
283

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

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

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

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

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

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

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

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

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

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

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

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

757
          expr_t id = it.second;
758

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

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

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

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

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

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