ModelTree.cc 77.8 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;

sebastien's avatar
sebastien committed
35
bool
36
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
37
{
38
  const int n = equations.size();
39
40
41

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

42
  using BipartiteGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
43
44
45
46
47
48
49
50

  /*
    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
51
  set<pair<int, int>> endo;
52

53
54
  for (const auto & it : contemporaneous_jacobian)
    add_edge(it.first.first + n, it.first.second, g);
55
56
57
58
59
60
61
62
63
64
65
66

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

67
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
68
69
70
71
72
73
74
75
76
77
78
79
    {
      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;
80
  while (not_maximum_yet)
81
82
83
84
85
86
87
88
89
90
91
    {
      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
92
  for (int i = 0; i < n; i++)
93
94
95
96
97
    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
98
  endo2eq.resize(equations.size());
99
100
101
102
103
104
105
106
  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;

107
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
    {
      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
126
  auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
127
  if (it != mate_map.begin() + n)
sebastien's avatar
sebastien committed
128
129
130
    {
      if (verbose)
        cerr << "ERROR: Could not normalize the model. Variable "
131
             << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
sebastien's avatar
sebastien committed
132
133
134
135
             << " is not in the maximum cardinality matching." << endl;
      check = false;
    }
  return check;
136
137
138
}

void
139
ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian)
140
{
sebastien's avatar
sebastien committed
141
142
  bool check = false;

143
144
  cout << "Normalizing the model..." << endl;

145
  int n = equations.size();
146

sebastien's avatar
sebastien committed
147
148
  // compute the maximum value of each row of the contemporaneous Jacobian matrix
  //jacob_map normalized_contemporaneous_jacobian;
149
  jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
sebastien's avatar
sebastien committed
150
  vector<double> max_val(n, 0.0);
151
  for (jacob_map_t::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++)
sebastien's avatar
sebastien committed
152
153
    if (fabs(iter->second) > max_val[iter->first.first])
      max_val[iter->first.first] = fabs(iter->second);
154

155
156
  for (auto & iter : normalized_contemporaneous_jacobian)
    iter.second /= max_val[iter.first.first];
sebastien's avatar
sebastien committed
157
158
159
160
161
162

  //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)
163
    {
164
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
sebastien's avatar
sebastien committed
165
      int suppress = 0;
166
167
      for (auto & iter : normalized_contemporaneous_jacobian)
        if (fabs(iter.second) > max(current_cutoff, cutoff))
168
          tmp_normalized_contemporaneous_jacobian[{ iter.first.first, iter.first.second }] = iter.second;
sebastien's avatar
sebastien committed
169
170
171
172
173
174
175
        else
          suppress++;

      if (suppress != suppressed)
        check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
      suppressed = suppress;
      if (!check)
176
        {
sebastien's avatar
sebastien committed
177
178
179
180
          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);
181
182
183
        }
    }

sebastien's avatar
sebastien committed
184
  if (!check)
185
    {
sebastien's avatar
sebastien committed
186
187
      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
188
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
189
      set<pair<int, int>> endo;
sebastien's avatar
sebastien committed
190
      for (int i = 0; i < n; i++)
191
192
193
        {
          endo.clear();
          equations[i]->collectEndogenous(endo);
194
          for (const auto & it : endo)
195
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
196
        }
sebastien's avatar
sebastien committed
197
198
      check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
      if (check)
199
        {
sebastien's avatar
sebastien committed
200
          // Update the jacobian matrix
201
          for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
sebastien's avatar
sebastien committed
202
            {
203
204
205
206
207
208
              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;
209
210
              try
                {
211
212
                  if (first_derivatives.find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == first_derivatives.end())
                    first_derivatives[{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
213
                }
Stéphane Adjemian's avatar
Stéphane Adjemian committed
214
              catch (DataTree::UnknownDerivIDException &e)
215
                {
216
                  cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it->first.second))
217
218
219
                       << " 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
220
            }
221
222
        }
    }
sebastien's avatar
sebastien committed
223
224
225
226
227
228

  if (!check)
    {
      cerr << "No normalization could be computed. Aborting." << endl;
      exit(EXIT_FAILURE);
    }
229
230
231
232
233
}

void
ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
{
234
  for (size_t i = 0; i < equations.size(); i++)
235
    {
236
      auto *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
237
      if (lhs == nullptr)
238
239
240
        continue;

      int symb_id = lhs->get_symb_id();
241
      if (symbol_table.getType(symb_id) != SymbolType::endogenous)
242
243
        continue;

244
      set<pair<int, int>> endo;
245
      equations[i]->get_arg2()->collectEndogenous(endo);
246
      if (endo.find({ symbol_table.getTypeSpecificID(symb_id), 0 }) != endo.end())
247
248
        continue;

249
      endo2eqs.emplace(symbol_table.getTypeSpecificID(symb_id), (int) i);
250
251
252
253
254
      cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
    }
}

void
255
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)
256
257
{
  int nb_elements_contemparenous_Jacobian = 0;
258
  set<pair<int, int>> jacobian_elements_to_delete;
259
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
260
261
262
       it != first_derivatives.end(); it++)
    {
      int deriv_id = it->first.second;
263
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
264
        {
265
          expr_t Id = it->second;
266
267
268
269
270
271
272
273
274
          int eq = it->first.first;
          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);
            }
275
276
277
278
          catch (ExprNode::EvalExternalFunctionException &e)
            {
              val = 1;
            }
279
280
          catch (ExprNode::EvalException &e)
            {
281
              cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
282
              Id->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, temporary_terms, {});
283
284
285
286
287
288
289
              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;
290
              jacobian_elements_to_delete.emplace(eq, deriv_id);
291
292
293
294
295
296
            }
          else
            {
              if (lag == 0)
                {
                  nb_elements_contemparenous_Jacobian++;
297
                  contemporaneous_jacobian[{ eq, var }] = val;
298
                }
299
300
              if (static_jacobian.find({ eq, var }) != static_jacobian.end())
                static_jacobian[{ eq, var }] += val;
301
              else
302
303
                static_jacobian[{ eq, var }] = val;
              dynamic_jacobian[{ lag, { eq, var } }] = Id;
304
305
306
307
308
            }
        }
    }

  // Get rid of the elements of the Jacobian matrix below the cutoff
309
310
  for (const auto & it : jacobian_elements_to_delete)
    first_derivatives.erase(it);
311

312
  if (jacobian_elements_to_delete.size() > 0)
313
314
315
316
317
318
319
    {
      cout << jacobian_elements_to_delete.size() << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
           << "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
    }
}

void
320
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
321
{
322
323
324
  vector<int> eq2endo(equations.size(), 0);
  equation_reordered.resize(equations.size());
  variable_reordered.resize(equations.size());
325
  bool *IM;
326
  int n = equations.size();
327
  IM = (bool *) calloc(n*n, sizeof(bool));
328
  int i = 0;
329
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
330
331
332
333
    {
      eq2endo[*it] = i;
      equation_reordered[i] = i;
      variable_reordered[*it] = i;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
334
335
336
    }
  if (cutoff == 0)
    {
337
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
338
      for (int i = 0; i < n; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
339
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
340
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
341
          equations[i]->collectEndogenous(endo);
342
343
          for (const auto & it : endo)
            IM[i * n + endo2eq[it.first]] = true;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
344
345
346
        }
    }
  else
347
348
    for (const auto & it : static_jacobian_arg)
      IM[it.first.first * n + endo2eq[it.first.second]] = true;
349
350
351
352
353
354
355
356
  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;
357
      for (int i = prologue; i < n; i++)
358
359
        {
          int nze = 0;
360
361
          for (int j = tmp_prologue; j < n; j++)
            if (IM[i * n + j])
362
              {
363
                nze++;
364
365
                k = j;
              }
366
          if (nze == 1)
367
            {
368
              for (int j = 0; j < n; j++)
369
370
371
372
373
374
375
376
                {
                  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;
377
              for (int j = 0; j < n; j++)
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
                {
                  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;
400
      for (int i = prologue; i < n - (int) epilogue; i++)
401
402
        {
          int nze = 0;
403
404
          for (int j = prologue; j < n - tmp_epilogue; j++)
            if (IM[j * n + i])
405
              {
406
                nze++;
407
408
                k = j;
              }
409
          if (nze == 1)
410
            {
411
              for (int j = 0; j < n; j++)
412
413
414
415
416
417
418
419
                {
                  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;
420
              for (int j = 0; j < n; j++)
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
                {
                  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);
}

438
equation_type_and_normalized_equation_t
439
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
440
{
441
  expr_t lhs;
442
443
  BinaryOpNode *eq_node;
  EquationType Equation_Simulation_Type;
444
  equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
445
446
447
448
449
450
451
  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;
452
      auto derivative = first_order_endo_derivatives.find({ eq, { var, 0 } });
453
      pair<bool, expr_t> res;
454
      if (derivative != first_order_endo_derivatives.end())
455
        {
456
          set<pair<int, int>> result;
457
          derivative->second->collectEndogenous(result);
458
          auto d_endo_variable = result.find({ var, 0 });
459
          //Determine whether the equation could be evaluated rather than to be solved
460
          if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
461
462
463
464
465
            {
              Equation_Simulation_Type = E_EVALUATE;
            }
          else
            {
466
              vector<pair<int, pair<expr_t, expr_t>>> List_of_Op_RHS;
467
              res =  equations[eq]->normalizeEquation(var, List_of_Op_RHS);
468
              if (mfs == 2)
469
                {
470
                  if (d_endo_variable == result.end() && res.second)
471
472
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
473
              else if (mfs == 3)
474
                {
475
                  if (res.second) // The equation could be solved analytically
476
477
478
479
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
            }
        }
480
      V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
481
482
483
484
485
    }
  return (V_Equation_Simulation_Type);
}

void
486
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
487
488
{
  int nb_endo = symbol_table.endo_nbr();
489
490
  variable_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
  equation_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
491
492
493
  vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
  for (int i = 0; i < nb_endo; i++)
    {
494
      if (i < (int) prologue)
495
496
497
498
        {
          variable_blck[variable_reordered[i]] = i;
          equation_blck[equation_reordered[i]] = i;
        }
499
      else if (i < (int) (components_set.size() + prologue))
500
501
502
503
504
505
506
507
508
509
        {
          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);
        }
    }
510
  for (const auto & it : dynamic_jacobian)
511
    {
512
513
514
      int lag = it.first.first;
      int j_1 = it.first.second.first;
      int i_1 = it.first.second.second;
515
516
517
      if (variable_blck[i_1] == equation_blck[j_1])
        {
          if (lag > variable_lead_lag[i_1].second)
518
            variable_lead_lag[i_1] = { variable_lead_lag[i_1].first, lag };
519
          if (lag < -variable_lead_lag[i_1].first)
520
            variable_lead_lag[i_1] = { -lag, variable_lead_lag[i_1].second };
521
          if (lag > equation_lead_lag[j_1].second)
522
            equation_lead_lag[j_1] = { equation_lead_lag[j_1].first, lag };
523
          if (lag < -equation_lead_lag[j_1].first)
524
            equation_lead_lag[j_1] = { -lag, equation_lead_lag[j_1].second };
525
526
527
528
529
        }
    }
}

void
530
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
531
532
533
534
{
  int nb_var = variable_reordered.size();
  int n = nb_var - prologue - epilogue;

535
536
537
  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
538
  auto v_index = get(boost::vertex_index, G2);
539
540
  for (int i = 0; i < n; i++)
    put(v_index, vertex(i, G2), i);
541
542
543

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

544
  for (int i = 0; i < nb_var; i++)
545
546
547
    {
      reverse_equation_reordered[equation_reordered[i]] = i;
      reverse_variable_reordered[variable_reordered[i]] = i;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
548
549
550
551
    }
  jacob_map_t tmp_normalized_contemporaneous_jacobian;
  if (cutoff == 0)
    {
552
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
553
      for (int i = 0; i < nb_var; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
554
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
555
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
556
          equations[i]->collectEndogenous(endo);
557
          for (const auto & it : endo)
558
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
559
560
561
562
563

        }
    }
  else
    tmp_normalized_contemporaneous_jacobian = static_jacobian;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
564
  for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
565
566
    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)
567
        && it->first.first != endo2eq[it->first.second])
568
569
570
      add_edge(vertex(reverse_equation_reordered[endo2eq[it->first.second]]-prologue, G2),
               vertex(reverse_equation_reordered[it->first.first]-prologue, G2),
               G2);
571
572

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
573
  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));
574
575

  // Compute strongly connected components
576
  int num = strong_components(G2, endo2block_map);
577

578
  blocks = vector<pair<int, int>>(num, { 0, 0 });
579
580

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

584
  for (unsigned int i = 0; i < num_vertices(G2); i++)
585
    {
586
587
      AdjacencyList_t::out_edge_iterator it_out, out_end;
      AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
588
589
      for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
        {
590
591
          int t_b = endo2block_map[target(*it_out, G2)];
          int s_b = endo2block_map[source(*it_out, G2)];
592
593
594
595
596
597
598
599
600
601
602
          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);
603
  for (int i = 0; i < num; i++)
604
605
606
607
608
609
    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.
610
  vector<pair<set<int>, pair<set<int>, vector<int>>>> components_set(num);
611
612
613
614
615
616
617
  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);
    }

618
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
619
620
621
622

  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
623
  if (select_feedback_variable)
624
625
626
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
627
628
629
630
631
            || 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)
632
          add_edge(vertex(i, G2), vertex(i, G2), G2);
633
634
635
636
637
    }
  else
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
638
          add_edge(vertex(i, G2), vertex(i, G2), G2);
639
    }
640
641
642
643
644
645
  //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);

646
  for (int i = 0; i < (int) prologue; i++)
647
648
649
650
651
652
653
654
655
656
    {
      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]++;
    }
657
658
659
660
661
662
  //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++)
    {
663
      AdjacencyList_t G = extract_subgraph(G2, components_set[i].first);
664
665
      set<int> feed_back_vertices;
      //Print(G);
666
      AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
667
      auto v_index = get(boost::vertex_index, G);
668
669
670
671
672
673
      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
674
      for (int j = 0; j < 4; j++)
675
        {
676
          for (int & its : Reordered_Vertice)
677
678
            {
              bool something_done = false;
679
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
680
681
682
683
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
684
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
685
686
687
688
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
689
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
690
691
692
693
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
694
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
695
696
697
698
699
700
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
701
702
                  equation_reordered[order] = tmp_equation_reordered[its+prologue];
                  variable_reordered[order] = tmp_variable_reordered[its+prologue];
703
704
705
                  order++;
                }
            }
706
707
708
        }
      components_set[i].second.second = Reordered_Vertice;
      //Second we have the equations related to the feedback variables
709
      for (int j = 0; j < 4; j++)
710
        {
711
          for (int feed_back_vertice : feed_back_vertices)
712
713
            {
              bool something_done = false;
714
              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)
715
716
717
718
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
719
              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)
720
721
722
723
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
724
              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)
725
726
727
728
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
729
              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)
730
731
732
733
734
735
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
736
737
                  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];
738
739
740
                  order++;
                }
            }
741
742
        }
    }
743

744
  for (int i = 0; i < (int) epilogue; i++)
745
    {
746
      if      (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
747
        n_mixed[prologue+num+i]++;
748
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
749
        n_forward[prologue+num+i]++;
750
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
751
        n_backward[prologue+num+i]++;
752
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
753
754
755
        n_static[prologue+num+i]++;
    }

756
757
  inv_equation_reordered = vector<int>(nb_var);
  inv_variable_reordered = vector<int>(nb_var);
758
  for (int i = 0; i < nb_var; i++)
759
760
761
762
763
764
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
}

765
void
766
ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
{
  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;
782
              Nb_feedback_variable = getBlockMfs(block);
783
784
785
786
787
788
789
790
791
792
793
            }
        }
    }

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

794
block_type_firstequation_size_mfs_t
795
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)
796
797
798
799
800
{
  int i = 0;
  int count_equ = 0, blck_count_simult = 0;
  int Blck_Size, MFS_Size;
  int Lead, Lag;
801
  block_type_firstequation_size_mfs_t block_type_size_mfs;
802
803
  BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
  int eq = 0;
804
805
806
807
  unsigned int l_n_static = 0;
  unsigned int l_n_forward = 0;
  unsigned int l_n_backward = 0;
  unsigned int l_n_mixed = 0;
808
  for (i = 0; i < (int) (prologue+blocks.size()+epilogue); i++)
809
810
    {
      int first_count_equ = count_equ;
811
      if (i < (int) prologue)
812
813
814
815
        {
          Blck_Size = 1;
          MFS_Size = 1;
        }
816
      else if (i < (int) (prologue+blocks.size()))
817
818
819
820
821
        {
          Blck_Size = blocks[blck_count_simult].first;
          MFS_Size = blocks[blck_count_simult].second;
          blck_count_simult++;
        }
822
      else if (i < (int) (prologue+blocks.size()+epilogue))
823
824
825
826
827
828
        {
          Blck_Size = 1;
          MFS_Size = 1;
        }

      Lag = Lead = 0;
829
      set<pair<int, int>> endo;
830
      for (count_equ  = first_count_equ; count_equ  < Blck_Size+first_count_equ; count_equ++)
831
832
833
        {
          endo.clear();
          equations[equation_reordered[count_equ]]->collectEndogenous(endo);
834
          for (const auto & it : endo)
835
            {
836
837
              int curr_variable = it.first;
              int curr_lag = it.second;
838
              auto it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
839
              if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
840
                if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end())
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
                  {
                    if (curr_lag > Lead)
                      Lead = curr_lag;
                    else if (-curr_lag > Lag)
                      Lag = -curr_lag;
                  }
            }
        }
      if ((Lag > 0) && (Lead > 0))
        {
          if (Blck_Size == 1)
            Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
          else
            Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
        }
      else if (Blck_Size > 1)
        {
          if (Lead > 0)
            Simulation_Type = SOLVE_BACKWARD_COMPLETE;
          else
            Simulation_Type = SOLVE_FORWARD_COMPLETE;
        }
      else
        {
          if (Lead > 0)
            Simulation_Type = SOLVE_BACKWARD_SIMPLE;
          else
            Simulation_Type = SOLVE_FORWARD_SIMPLE;
        }
870
871
872
873
      l_n_static = n_static[i];
      l_n_forward = n_forward[i];
      l_n_backward = n_backward[i];
      l_n_mixed = n_mixed[i];
874
875
      if (Blck_Size == 1)
        {
876
          if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
877
878
879
880
881
882
883
884
            {
              if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
                Simulation_Type = EVALUATE_BACKWARD;
              else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
                Simulation_Type = EVALUATE_FORWARD;
            }
          if (i > 0)
            {
885
886
887
888
              bool is_lead = false, is_lag = false;
              int c_Size = (block_type_size_mfs[block_type_size_mfs.size()-1]).second.first;
              int first_equation = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.second;
              if (c_Size > 0 && ((prev_Type ==  EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
889
                                 || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))
890
891
892
                {
                  for (int j = first_equation; j < first_equation+c_Size; j++)
                    {
893
                      auto it = dynamic_jacobian.find({ -1, { equation_reordered[eq], variable_reordered[j] } });
894
895
                      if (it != dynamic_jacobian.end())
                        is_lag = true;
896
                      it = dynamic_jacobian.find({ +1, { equation_reordered[eq], variable_reordered[j] } });
897
898
899
900
901
902
                      if (it != dynamic_jacobian.end())
                        is_lead = true;
                    }
                }
              if ((prev_Type ==  EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
                  || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag))
903
904
905
                {
                  //merge the current block with the previous one
                  BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
906
                  c_Size++;
907
                  block_type_size_mfs[block_type_size_mfs.size()-1] = { { c_Type, first_equation }, { c_Size, c_Size } };
908
                  if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
909
                    Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
910
                  if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
911
                    Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
912
                  block_lag_lead[block_type_size_mfs.size()-1] = { Lag, Lead };
913
                  pair< pair< unsigned int, unsigned int>, pair<unsigned int, unsigned int>> tmp = block_col_type[block_col_type.size()-1];
914
                  block_col_type[block_col_type.size()-1] = { { tmp.first.first+l_n_static, tmp.first.second+l_n_forward }, { tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed } };
915
916
917
                }
              else
                {
918
919
                  block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
                  block_lag_lead.emplace_back(Lag, Lead);
920
                  block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
921
922
923
924
                }
            }
          else
            {