ModelTree.cc 89.7 KB
Newer Older
1
/*
Houtan Bastani's avatar
Houtan Bastani committed
2
 * Copyright (C) 2003-2018 Dynare Team
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 <cstdlib>
21
#include <cassert>
sebastien's avatar
sebastien committed
22
#include <cmath>
23
#include <iostream>
24
#include <fstream>
25
26

#include "ModelTree.hh"
27
28
29
30
31
32
33
34
#include "MinimumFeedbackSet.hh"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>

using namespace MFS;

35
36
37
void
ModelTree::copyHelper(const ModelTree &m)
{
38
  auto f = [this](expr_t e) { return e->clone(*this); };
39
40
41
42
43
44
45

  // Equations
  for (const auto & it : m.equations)
    equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
  for (const auto & it : m.aux_equations)
    aux_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));

46
47
48
49
50
51
52
53
  auto convert_deriv_map = [f](map<vector<int>, expr_t> dm)
    {
      map<vector<int>, expr_t> dm2;
      for (const auto &it : dm)
        dm2.emplace(it.first, f(it.second));
      return dm2;
    };

54
  // Derivatives
55
56
57
58
59
60
61
62
63
64
65
66
  for (const auto &it : m.derivatives)
    derivatives.push_back(convert_deriv_map(it));
  for (const auto &it : m.params_derivatives)
    params_derivatives[it.first] = convert_deriv_map(it.second);

  auto convert_temporary_terms_t = [f](temporary_terms_t tt)
    {
      temporary_terms_t tt2;
      for (const auto &it : tt)
        tt2.insert(f(it));
      return tt2;
    };
67
68
69
70
71
72

  // Temporary terms
  for (const auto & it : m.temporary_terms)
    temporary_terms.insert(f(it));
  for (const auto & it : m.temporary_terms_mlv)
    temporary_terms_mlv[f(it.first)] = f(it.second);
73
74
  for (const auto &it : m.temporary_terms_derivatives)
    temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
75
76
77
  for (const auto & it : m.temporary_terms_idxs)
    temporary_terms_idxs[f(it.first)] = it.second;
  for (const auto & it : m.params_derivs_temporary_terms)
78
    params_derivs_temporary_terms[it.first] = convert_temporary_terms_t(it.second);
79
80
81
82
83
84
85
86
87
88
  for (const auto & it : m.params_derivs_temporary_terms_idxs)
    params_derivs_temporary_terms_idxs[f(it.first)] = it.second;

  // Other stuff
  for (const auto & it : m.trend_symbols_map)
    trend_symbols_map[it.first] = f(it.second);
  for (const auto & it : m.nonstationary_symbols_map)
    nonstationary_symbols_map[it.first] = make_pair(it.second.first, f(it.second.second));
}

89
90
91
92
93
94
95
96
97
98
99
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
                     NumericalConstants &num_constants_arg,
                     ExternalFunctionsTable &external_functions_table_arg,
                     bool is_dynamic_arg) :
  DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
  derivatives(4),
  NNZDerivatives(4, 0),
  temporary_terms_derivatives(4)
{
}

100
101
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
ModelTree::ModelTree(const ModelTree &m) :
  DataTree {m},
  equations_lineno {m.equations_lineno},
  equation_tags {m.equation_tags},
  NNZDerivatives {m.NNZDerivatives},
  equation_reordered {m.equation_reordered},
  variable_reordered {m.variable_reordered},
  inv_equation_reordered {m.inv_equation_reordered},
  inv_variable_reordered {m.inv_variable_reordered},
  is_equation_linear {m.is_equation_linear},
  endo2eq {m.endo2eq},
  epilogue {m.epilogue},
  prologue {m.prologue},
  block_lag_lead {m.block_lag_lead},
  cutoff {m.cutoff},
  mfs {m.mfs}
{
  copyHelper(m);
}

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

  equations.clear();
  equations_lineno = m.equations_lineno;
  aux_equations.clear();
  equation_tags = m.equation_tags;
  NNZDerivatives = m.NNZDerivatives;

131
132
  derivatives.clear();
  params_derivatives.clear();
133
134
135

  temporary_terms.clear();
  temporary_terms_mlv.clear();
136
  temporary_terms_derivatives.clear();
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
  params_derivs_temporary_terms.clear();
  params_derivs_temporary_terms_idxs.clear();

  trend_symbols_map.clear();
  nonstationary_symbols_map.clear();

  equation_reordered = m.equation_reordered;
  variable_reordered = m.variable_reordered;
  inv_equation_reordered = m.inv_equation_reordered;
  inv_variable_reordered = m.inv_variable_reordered;
  is_equation_linear = m.is_equation_linear;
  endo2eq = m.endo2eq;
  epilogue = m.epilogue;
  prologue = m.prologue;
  block_lag_lead = m.block_lag_lead;
  cutoff = m.cutoff;
  mfs = m.mfs;

  copyHelper(m);

  return *this;
}


sebastien's avatar
sebastien committed
161
bool
162
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
163
{
164
  const int n = equations.size();
165
166
167

  assert(n == symbol_table.endo_nbr());

168
  using BipartiteGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
169
170
171
172
173
174
175
176

  /*
    Vertices 0 to n-1 are for endogenous (using type specific ID)
    Vertices n to 2*n-1 are for equations (using equation no.)
  */
  BipartiteGraph g(2 * n);

  // Fill in the graph
177
  set<pair<int, int>> endo;
178

179
180
  for (const auto & it : contemporaneous_jacobian)
    add_edge(it.first.first + n, it.first.second, g);
181
182
183
184
185
186
187
188
189
190
191
192

  // Compute maximum cardinality matching
  vector<int> mate_map(2*n);

#if 1
  bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
#else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
  fill(mate_map.begin(), mate_map.end(), graph_traits<BipartiteGraph>::null_vertex());

  multimap<int, int> natural_endo2eqs;
  computeNormalizedEquations(natural_endo2eqs);

193
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
194
195
196
197
198
199
200
201
202
203
204
205
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

      int j = natural_endo2eqs.find(i)->second;

      put(&mate_map[0], i, n+j);
      put(&mate_map[0], n+j, i);
    }

  edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
  bool not_maximum_yet = true;
206
  while (not_maximum_yet)
207
208
209
210
211
212
213
214
215
216
217
    {
      not_maximum_yet = augmentor.augment_matching();
    }
  augmentor.get_current_matching(&mate_map[0]);

  bool check = maximum_cardinality_matching_verifier<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(vertex_index, g));
#endif

  assert(check);

#ifdef DEBUG
218
  for (int i = 0; i < n; i++)
219
220
221
222
223
    cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
         << " matched with equation " << (mate_map[i]-n+1) << endl;
#endif

  // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
224
  endo2eq.resize(equations.size());
225
226
227
228
229
230
231
232
  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));

#ifdef DEBUG
  multimap<int, int> natural_endo2eqs;
  computeNormalizedEquations(natural_endo2eqs);

  int n1 = 0, n2 = 0;

233
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

      n1++;

      pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
      if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
        cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
             << " not used." << endl;
      else
        n2++;
    }

  cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
#endif

  // Check if all variables are normalized
252
  auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
253
  if (it != mate_map.begin() + n)
sebastien's avatar
sebastien committed
254
255
256
    {
      if (verbose)
        cerr << "ERROR: Could not normalize the model. Variable "
257
             << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
sebastien's avatar
sebastien committed
258
259
260
261
             << " is not in the maximum cardinality matching." << endl;
      check = false;
    }
  return check;
262
263
264
}

void
265
ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian)
266
{
sebastien's avatar
sebastien committed
267
268
  bool check = false;

269
270
  cout << "Normalizing the model..." << endl;

271
  int n = equations.size();
272

sebastien's avatar
sebastien committed
273
274
  // compute the maximum value of each row of the contemporaneous Jacobian matrix
  //jacob_map normalized_contemporaneous_jacobian;
275
  jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
sebastien's avatar
sebastien committed
276
  vector<double> max_val(n, 0.0);
277
  for (jacob_map_t::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++)
sebastien's avatar
sebastien committed
278
279
    if (fabs(iter->second) > max_val[iter->first.first])
      max_val[iter->first.first] = fabs(iter->second);
280

281
282
  for (auto & iter : normalized_contemporaneous_jacobian)
    iter.second /= max_val[iter.first.first];
sebastien's avatar
sebastien committed
283
284
285
286
287
288

  //We start with the highest value of the cutoff and try to normalize the model
  double current_cutoff = 0.99999999;

  int suppressed = 0;
  while (!check && current_cutoff > 1e-19)
289
    {
290
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
sebastien's avatar
sebastien committed
291
      int suppress = 0;
292
293
      for (auto & iter : normalized_contemporaneous_jacobian)
        if (fabs(iter.second) > max(current_cutoff, cutoff))
294
          tmp_normalized_contemporaneous_jacobian[{ iter.first.first, iter.first.second }] = iter.second;
sebastien's avatar
sebastien committed
295
296
297
298
299
300
301
        else
          suppress++;

      if (suppress != suppressed)
        check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
      suppressed = suppress;
      if (!check)
302
        {
sebastien's avatar
sebastien committed
303
304
305
306
          current_cutoff /= 2;
          // In this last case try to normalize with the complete jacobian
          if (current_cutoff <= 1e-19)
            check = computeNormalization(normalized_contemporaneous_jacobian, false);
307
308
309
        }
    }

sebastien's avatar
sebastien committed
310
  if (!check)
311
    {
sebastien's avatar
sebastien committed
312
313
      cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
      //if no non-singular normalization can be found, try to find a normalization even with a potential singularity
314
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
315
      set<pair<int, int>> endo;
sebastien's avatar
sebastien committed
316
      for (int i = 0; i < n; i++)
317
318
319
        {
          endo.clear();
          equations[i]->collectEndogenous(endo);
320
          for (const auto & it : endo)
321
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
322
        }
sebastien's avatar
sebastien committed
323
324
      check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
      if (check)
325
        {
sebastien's avatar
sebastien committed
326
          // Update the jacobian matrix
327
          for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
sebastien's avatar
sebastien committed
328
            {
329
330
331
332
333
334
              if (static_jacobian.find({ it->first.first, it->first.second }) == static_jacobian.end())
                static_jacobian[{ it->first.first, it->first.second }] = 0;
              if (dynamic_jacobian.find({ 0, { it->first.first, it->first.second } }) == dynamic_jacobian.end())
                dynamic_jacobian[{ 0, { it->first.first, it->first.second } }] = nullptr;
              if (contemporaneous_jacobian.find({ it->first.first, it->first.second }) == contemporaneous_jacobian.end())
                contemporaneous_jacobian[{ it->first.first, it->first.second }] = 0;
335
336
              try
                {
337
338
                  if (derivatives[1].find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == derivatives[1].end())
                    derivatives[1][{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
339
                }
Stéphane Adjemian's avatar
Stéphane Adjemian committed
340
              catch (DataTree::UnknownDerivIDException &e)
341
                {
342
                  cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it->first.second))
343
344
345
                       << " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
                  exit(EXIT_FAILURE);
                }
sebastien's avatar
sebastien committed
346
            }
347
348
        }
    }
sebastien's avatar
sebastien committed
349
350
351
352
353
354

  if (!check)
    {
      cerr << "No normalization could be computed. Aborting." << endl;
      exit(EXIT_FAILURE);
    }
355
356
357
358
359
}

void
ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
{
360
  for (size_t i = 0; i < equations.size(); i++)
361
    {
362
      auto *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
363
      if (lhs == nullptr)
364
365
366
        continue;

      int symb_id = lhs->get_symb_id();
367
      if (symbol_table.getType(symb_id) != SymbolType::endogenous)
368
369
        continue;

370
      set<pair<int, int>> endo;
371
      equations[i]->get_arg2()->collectEndogenous(endo);
372
      if (endo.find({ symbol_table.getTypeSpecificID(symb_id), 0 }) != endo.end())
373
374
        continue;

375
      endo2eqs.emplace(symbol_table.getTypeSpecificID(symb_id), (int) i);
376
377
378
379
380
      cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
    }
}

void
381
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose)
382
383
{
  int nb_elements_contemparenous_Jacobian = 0;
384
385
386
  set<vector<int>> jacobian_elements_to_delete;
  for (auto it = derivatives[1].begin();
       it != derivatives[1].end(); it++)
387
    {
388
      int deriv_id = it->first[1];
389
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
390
        {
391
          expr_t Id = it->second;
392
          int eq = it->first[0];
393
394
395
396
397
398
399
400
          int symb = getSymbIDByDerivID(deriv_id);
          int var = symbol_table.getTypeSpecificID(symb);
          int lag = getLagByDerivID(deriv_id);
          double val = 0;
          try
            {
              val = Id->eval(eval_context);
            }
401
402
403
404
          catch (ExprNode::EvalExternalFunctionException &e)
            {
              val = 1;
            }
405
406
          catch (ExprNode::EvalException &e)
            {
407
              cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
408
              Id->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, temporary_terms, {});
409
410
411
412
413
414
415
              cerr << endl;
              exit(EXIT_FAILURE);
            }
          if (fabs(val) < cutoff)
            {
              if (verbose)
                cout << "the coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")" << endl;
416
              jacobian_elements_to_delete.insert({ eq, deriv_id });
417
418
419
420
421
422
            }
          else
            {
              if (lag == 0)
                {
                  nb_elements_contemparenous_Jacobian++;
423
                  contemporaneous_jacobian[{ eq, var }] = val;
424
                }
425
426
              if (static_jacobian.find({ eq, var }) != static_jacobian.end())
                static_jacobian[{ eq, var }] += val;
427
              else
428
429
                static_jacobian[{ eq, var }] = val;
              dynamic_jacobian[{ lag, { eq, var } }] = Id;
430
431
432
433
434
            }
        }
    }

  // Get rid of the elements of the Jacobian matrix below the cutoff
435
  for (const auto & it : jacobian_elements_to_delete)
436
    derivatives[1].erase(it);
437

438
  if (jacobian_elements_to_delete.size() > 0)
439
    {
440
      cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
441
442
443
444
           << "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
    }
}

445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
vector<pair<int, int> >
ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered,
                                                     vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered,
                                                     lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead,
                                                     vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed)
{
  vector<int> eq2endo(equations.size(), 0);
  /*equation_reordered.resize(equations.size());
  variable_reordered.resize(equations.size());*/
  unsigned int num = 0;
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++)
    if (!is_equation_linear[*it])
      num++;
  vector<int> endo2block = vector<int>(endo2eq.size(), 1);
  vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
  int i = 0, j = 0;
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, j++)
    {
      if (!is_equation_linear[*it])
        {
          equation_reordered[i] = *it;
          variable_reordered[i] = j;
          endo2block[j] = 0;
          components_set[endo2block[j]].first.insert(i);
          /*cout << " -----------------------------------------" << endl;
          cout.flush();
          cout << "equation_reordered[" << *it << "] = " << i << endl;
          cout.flush();
          cout << "variable_reordered[" << j << "] = " << i << endl;
          cout.flush();
          cout << "components_set[" << endo2block[j] << "].first[" << i << "] = " << i << endl;
          cout << "endo2block[" << j << "] = " << 0 << endl;
          cout.flush();
          */
          i++;
        }
    }
/*  for (unsigned int j = 0; j < is_equation_linear.size() ; j++)
     cout << "endo2block[" << j << "] = " << endo2block[j] << endl;*/
  /*cout << "before getVariableLeadLAgByBlock\n";
  cout.flush();*/
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, endo2block.size(), equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
  n_static = vector<unsigned int>(endo2eq.size(), 0);
  n_forward = vector<unsigned int>(endo2eq.size(), 0);
  n_backward = vector<unsigned int>(endo2eq.size(), 0);
  n_mixed = vector<unsigned int>(endo2eq.size(), 0);
  for (unsigned int i = 0; i < endo2eq.size(); i++)
    {
      if      (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second != 0)
        n_mixed[i]++;
      else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second != 0)
        n_forward[i]++;
      else if (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second == 0)
        n_backward[i]++;
      else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second == 0)
        n_static[i]++;
    }
  cout.flush();
  int nb_endo = is_equation_linear.size();
  vector<pair<int, int> > blocks = vector<pair<int, int> >(1, make_pair(i, i));
  inv_equation_reordered = vector<int>(nb_endo);
  inv_variable_reordered = vector<int>(nb_endo);
  for (int i = 0; i < nb_endo; i++)
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
  return blocks;
}

bool
ModelTree::computeNaturalNormalization()
{
  bool bool_result = true;
  set<pair<int, int> > result;
  endo2eq.resize(equations.size());
  for (int eq = 0; eq < (int) equations.size(); eq++)
    if (!is_equation_linear[eq])
      {
        BinaryOpNode *eq_node = equations[eq];
        expr_t lhs = eq_node->get_arg1();
        result.clear();
        lhs->collectDynamicVariables(SymbolType::endogenous, result);
        if (result.size() == 1 && result.begin()->second == 0)
          {
            //check if the endogenous variable has not been already used in an other match !
             vector<int>::iterator it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
             if (it == endo2eq.end())
               endo2eq[result.begin()->first] = eq;
             else
              {
                bool_result = false;
                break;
              }
          }
      }
  return bool_result;
}

544
void
545
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
546
{
547
548
549
  vector<int> eq2endo(equations.size(), 0);
  equation_reordered.resize(equations.size());
  variable_reordered.resize(equations.size());
550
  bool *IM;
551
  int n = equations.size();
552
  IM = (bool *) calloc(n*n, sizeof(bool));
553
  int i = 0;
554
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
555
556
557
558
    {
      eq2endo[*it] = i;
      equation_reordered[i] = i;
      variable_reordered[*it] = i;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
559
560
561
    }
  if (cutoff == 0)
    {
562
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
563
      for (int i = 0; i < n; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
564
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
565
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
566
          equations[i]->collectEndogenous(endo);
567
568
          for (const auto & it : endo)
            IM[i * n + endo2eq[it.first]] = true;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
569
570
571
        }
    }
  else
572
573
    for (const auto & it : static_jacobian_arg)
      IM[it.first.first * n + endo2eq[it.first.second]] = true;
574
575
576
577
578
579
580
581
  bool something_has_been_done = true;
  prologue = 0;
  int k = 0;
  // Find the prologue equations and place first the AR(1) shock equations first
  while (something_has_been_done)
    {
      int tmp_prologue = prologue;
      something_has_been_done = false;
582
      for (int i = prologue; i < n; i++)
583
584
        {
          int nze = 0;
585
586
          for (int j = tmp_prologue; j < n; j++)
            if (IM[i * n + j])
587
              {
588
                nze++;
589
590
                k = j;
              }
591
          if (nze == 1)
592
            {
593
              for (int j = 0; j < n; j++)
594
595
596
597
598
599
600
601
                {
                  bool tmp_bool = IM[tmp_prologue * n + j];
                  IM[tmp_prologue * n + j] = IM[i * n + j];
                  IM[i * n + j] = tmp_bool;
                }
              int tmp = equation_reordered[tmp_prologue];
              equation_reordered[tmp_prologue] = equation_reordered[i];
              equation_reordered[i] = tmp;
602
              for (int j = 0; j < n; j++)
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
                {
                  bool tmp_bool = IM[j * n + tmp_prologue];
                  IM[j * n + tmp_prologue] = IM[j * n + k];
                  IM[j * n + k] = tmp_bool;
                }
              tmp = variable_reordered[tmp_prologue];
              variable_reordered[tmp_prologue] = variable_reordered[k];
              variable_reordered[k] = tmp;
              tmp_prologue++;
              something_has_been_done = true;
            }
        }
      prologue = tmp_prologue;
    }

  something_has_been_done = true;
  epilogue = 0;
  // Find the epilogue equations
  while (something_has_been_done)
    {
      int tmp_epilogue = epilogue;
      something_has_been_done = false;
625
      for (int i = prologue; i < n - (int) epilogue; i++)
626
627
        {
          int nze = 0;
628
629
          for (int j = prologue; j < n - tmp_epilogue; j++)
            if (IM[j * n + i])
630
              {
631
                nze++;
632
633
                k = j;
              }
634
          if (nze == 1)
635
            {
636
              for (int j = 0; j < n; j++)
637
638
639
640
641
642
643
644
                {
                  bool tmp_bool = IM[(n - 1 - tmp_epilogue) * n + j];
                  IM[(n - 1 - tmp_epilogue) * n + j] = IM[k * n + j];
                  IM[k * n + j] = tmp_bool;
                }
              int tmp = equation_reordered[n - 1 - tmp_epilogue];
              equation_reordered[n - 1 - tmp_epilogue] = equation_reordered[k];
              equation_reordered[k] = tmp;
645
              for (int j = 0; j < n; j++)
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
                {
                  bool tmp_bool = IM[j * n + n - 1 - tmp_epilogue];
                  IM[j * n + n - 1 - tmp_epilogue] = IM[j * n + i];
                  IM[j * n + i] = tmp_bool;
                }
              tmp = variable_reordered[n - 1 - tmp_epilogue];
              variable_reordered[n - 1 - tmp_epilogue] = variable_reordered[i];
              variable_reordered[i] = tmp;
              tmp_epilogue++;
              something_has_been_done = true;
            }
        }
      epilogue = tmp_epilogue;
    }
  free(IM);
}

663
equation_type_and_normalized_equation_t
664
ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
665
{
666
  expr_t lhs;
667
668
  BinaryOpNode *eq_node;
  EquationType Equation_Simulation_Type;
669
  equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
670
671
672
673
674
675
676
  for (unsigned int i = 0; i < equations.size(); i++)
    {
      int eq = Index_Equ_IM[i];
      int var = Index_Var_IM[i];
      eq_node = equations[eq];
      lhs = eq_node->get_arg1();
      Equation_Simulation_Type = E_SOLVE;
677
      auto derivative = first_order_endo_derivatives.find({ eq, { var, 0 } });
678
      pair<bool, expr_t> res;
679
      if (derivative != first_order_endo_derivatives.end())
680
        {
681
          set<pair<int, int>> result;
682
          derivative->second->collectEndogenous(result);
683
          auto d_endo_variable = result.find({ var, 0 });
684
          //Determine whether the equation could be evaluated rather than to be solved
685
          if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
686
687
688
689
690
            {
              Equation_Simulation_Type = E_EVALUATE;
            }
          else
            {
691
              vector<pair<int, pair<expr_t, expr_t>>> List_of_Op_RHS;
692
              res =  equations[eq]->normalizeEquation(var, List_of_Op_RHS);
693
              if (mfs == 2)
694
                {
695
                  if (d_endo_variable == result.end() && res.second)
696
697
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
698
              else if (mfs == 3)
699
                {
700
                  if (res.second) // The equation could be solved analytically
701
702
703
704
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
            }
        }
705
      V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
706
707
708
709
710
    }
  return (V_Equation_Simulation_Type);
}

void
711
ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const
712
713
{
  int nb_endo = symbol_table.endo_nbr();
714
715
  variable_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
  equation_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
716
717
718
  vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
  for (int i = 0; i < nb_endo; i++)
    {
719
      if (i < (int) prologue)
720
721
722
723
        {
          variable_blck[variable_reordered[i]] = i;
          equation_blck[equation_reordered[i]] = i;
        }
724
      else if (i < (int) (components_set.size() + prologue))
725
726
727
728
729
730
731
732
733
734
        {
          variable_blck[variable_reordered[i]] = components_set[i-prologue] + prologue;
          equation_blck[equation_reordered[i]] = components_set[i-prologue] + prologue;
        }
      else
        {
          variable_blck[variable_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
          equation_blck[equation_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
        }
    }
735
  for (const auto & it : dynamic_jacobian)
736
    {
737
738
739
      int lag = it.first.first;
      int j_1 = it.first.second.first;
      int i_1 = it.first.second.second;
740
741
742
      if (variable_blck[i_1] == equation_blck[j_1])
        {
          if (lag > variable_lead_lag[i_1].second)
743
            variable_lead_lag[i_1] = { variable_lead_lag[i_1].first, lag };
744
          if (lag < -variable_lead_lag[i_1].first)
745
            variable_lead_lag[i_1] = { -lag, variable_lead_lag[i_1].second };
746
          if (lag > equation_lead_lag[j_1].second)
747
            equation_lead_lag[j_1] = { equation_lead_lag[j_1].first, lag };
748
          if (lag < -equation_lead_lag[j_1].first)
749
            equation_lead_lag[j_1] = { -lag, equation_lead_lag[j_1].second };
750
751
752
753
754
        }
    }
}

void
755
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const
756
757
758
759
{
  int nb_var = variable_reordered.size();
  int n = nb_var - prologue - epilogue;

760
761
762
  AdjacencyList_t G2(n);

  // It is necessary to manually initialize vertex_index property since this graph uses listS and not vecS as underlying vertex container
763
  auto v_index = get(boost::vertex_index, G2);
764
765
  for (int i = 0; i < n; i++)
    put(v_index, vertex(i, G2), i);
766
767
768

  vector<int> reverse_equation_reordered(nb_var), reverse_variable_reordered(nb_var);

769
  for (int i = 0; i < nb_var; i++)
770
771
772
    {
      reverse_equation_reordered[equation_reordered[i]] = i;
      reverse_variable_reordered[variable_reordered[i]] = i;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
773
774
775
776
    }
  jacob_map_t tmp_normalized_contemporaneous_jacobian;
  if (cutoff == 0)
    {
777
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
778
      for (int i = 0; i < nb_var; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
779
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
780
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
781
          equations[i]->collectEndogenous(endo);
782
          for (const auto & it : endo)
783
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
784
785
786
787
788

        }
    }
  else
    tmp_normalized_contemporaneous_jacobian = static_jacobian;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
789
  for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
790
791
    if (reverse_equation_reordered[it->first.first] >= (int) prologue && reverse_equation_reordered[it->first.first] < (int) (nb_var - epilogue)
        && reverse_variable_reordered[it->first.second] >= (int) prologue && reverse_variable_reordered[it->first.second] < (int) (nb_var - epilogue)
792
        && it->first.first != endo2eq[it->first.second])
793
794
795
      add_edge(vertex(reverse_equation_reordered[endo2eq[it->first.second]]-prologue, G2),
               vertex(reverse_equation_reordered[it->first.first]-prologue, G2),
               G2);
796
797

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
798
  boost::iterator_property_map<int *, boost::property_map<AdjacencyList_t, boost::vertex_index_t>::type, int, int &> endo2block_map(&endo2block[0], get(boost::vertex_index, G2));
799
800

  // Compute strongly connected components
801
  int num = strong_components(G2, endo2block_map);
802

803
  blocks = vector<pair<int, int>>(num, { 0, 0 });
804
805

  // Create directed acyclic graph associated to the strongly connected components
806
  using DirectedGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS>;
807
808
  DirectedGraph dag(num);

809
  for (unsigned int i = 0; i < num_vertices(G2); i++)
810
    {
811
812
      AdjacencyList_t::out_edge_iterator it_out, out_end;
      AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
813
814
      for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
        {
815
816
          int t_b = endo2block_map[target(*it_out, G2)];
          int s_b = endo2block_map[source(*it_out, G2)];
817
818
819
820
821
822
823
824
825
826
827
          if (s_b != t_b)
            add_edge(s_b, t_b, dag);
        }
    }

  // Compute topological sort of DAG (ordered list of unordered SCC)
  deque<int> ordered2unordered;
  topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order

  // Construct mapping from unordered SCC to ordered SCC
  vector<int> unordered2ordered(num);
828
  for (int i = 0; i < num; i++)
829
830
831
832
833
834
    unordered2ordered[ordered2unordered[i]] = i;

  //This vector contains for each block:
  //   - first set = equations belonging to the block,
  //   - second set = the feeback variables,
  //   - third vector = the reordered non-feedback variables.
835
  vector<pair<set<int>, pair<set<int>, vector<int>>>> components_set(num);
836
837
838
839
840
841
842
  for (unsigned int i = 0; i < endo2block.size(); i++)
    {
      endo2block[i] = unordered2ordered[endo2block[i]];
      blocks[endo2block[i]].first++;
      components_set[endo2block[i]].first.insert(i);
    }

843
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
844
845
846
847

  vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
  int order = prologue;
  //Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
848
  if (select_feedback_variable)
849
850
851
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
852
853
854
855
856
            || variable_lag_lead[variable_reordered[i+prologue]].second > 0
            || variable_lag_lead[variable_reordered[i+prologue]].first > 0
            || equation_lag_lead[equation_reordered[i+prologue]].second > 0
            || equation_lag_lead[equation_reordered[i+prologue]].first > 0
            || mfs == 0)
857
          add_edge(vertex(i, G2), vertex(i, G2), G2);
858
859
860
861
862
    }
  else
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
863
          add_edge(vertex(i, G2), vertex(i, G2), G2);
864
    }
865
866
867
868
869
870
  //Determines the dynamic structure of each equation
  n_static = vector<unsigned int>(prologue+num+epilogue, 0);
  n_forward = vector<unsigned int>(prologue+num+epilogue, 0);
  n_backward = vector<unsigned int>(prologue+num+epilogue, 0);
  n_mixed = vector<unsigned int>(prologue+num+epilogue, 0);

871
  for (int i = 0; i < (int) prologue; i++)
872
873
874
875
876
877
878
879
880
881
    {
      if      (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
        n_mixed[i]++;
      else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
        n_forward[i]++;
      else if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
        n_backward[i]++;
      else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
        n_static[i]++;
    }
882
883
884
885
886
887
  //For each block, the minimum set of feedback variable is computed
  // and the non-feedback variables are reordered to get
  // a sub-recursive block without feedback variables

  for (int i = 0; i < num; i++)
    {
888
      AdjacencyList_t G = extract_subgraph(G2, components_set[i].first);
889
890
      set<int> feed_back_vertices;
      //Print(G);
891
      AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
892
      auto v_index = get(boost::vertex_index, G);
893
894
895
896
897
898
      components_set[i].second.first = feed_back_vertices;
      blocks[i].second = feed_back_vertices.size();
      vector<int> Reordered_Vertice;
      Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);

      //First we have the recursive equations conditional on feedback variables
899
      for (int j = 0; j < 4; j++)
900
        {
901
          for (int & its : Reordered_Vertice)
902
903
            {
              bool something_done = false;
904
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
905
906
907
908
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
909
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
910
911
912
913
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
914
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
915
916
917
918
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
919
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
920
921
922
923
924
925
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
926
927
                  equation_reordered[order] = tmp_equation_reordered[its+prologue];
                  variable_reordered[order] = tmp_variable_reordered[its+prologue];
928
929
930
                  order++;
                }
            }
931
932
933
        }
      components_set[i].second.second = Reordered_Vertice;
      //Second we have the equations related to the feedback variables
934
      for (int j = 0; j < 4; j++)
935
        {
936
          for (int feed_back_vertice : feed_back_vertices)
937
938
            {
              bool something_done = false;
939
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
940
941
942
943
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
944
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
945
946
947
948
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
949
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
950
951
952
953
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
954
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
955
956
957
958
959
960
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
961
962
                  equation_reordered[order] = tmp_equation_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
                  variable_reordered[order] = tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
963
964
965
                  order++;
                }
            }
966
967
        }
    }
968

969
  for (int i = 0; i < (int) epilogue; i++)
970
    {
971
      if      (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
972
        n_mixed[prologue+num+i]++;
973
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
974
        n_forward[prologue+num+i]++;
975
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
976
        n_backward[prologue+num+i]++;
977
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
978
979
980
        n_static[prologue+num+i]++;
    }

981
982
  inv_equation_reordered = vector<int>(nb_var);
  inv_variable_reordered = vector<int>(nb_var);
983
  for (int i = 0; i < nb_var; i++)
984
985
986
987
988
989
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
}

990
void
991
ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
{
  int largest_block = 0;
  int Nb_SimulBlocks = 0;
  int Nb_feedback_variable = 0;
  unsigned int Nb_TotalBlocks = getNbBlocks();
  for (unsigned int block = 0; block < Nb_TotalBlocks; block++)
    {
      BlockSimulationType simulation_type = getBlockSimulationType(block);
      if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
        {
          Nb_SimulBlocks++;
          int size = getBlockSize(block);
          if (size > largest_block)
            {
              largest_block = size;
1007
              Nb_feedback_variable = getBlockMfs(block);
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
            }
        }
    }

  int Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
  cout << Nb_TotalBlocks << " block(s) found:" << endl
       << "  " << Nb_RecursBlocks << " recursive block(s) and " << Nb_SimulBlocks << " simultaneous block(s)." << endl
       << "  the largest simultaneous block has " << largest_block << " equation(s)" << endl
       << "                                 and " << Nb_feedback_variable << " feedback variable(s)." << endl;
}

1019
block_type_firstequation_size_mfs_t
1020
ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int>>> &block_col_type, bool linear_decomposition)
1021
1022
1023
1024
1025
{
  int i = 0;
  int count_equ = 0, blck_count_simult = 0;
  int Blck_Size, MFS_Size;