ModelTree.cc 78 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
35
#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 boost;
using namespace MFS;

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

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

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

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

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

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

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

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

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

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

146
  int n = equations.size();
147

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

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

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

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

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

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

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

      int symb_id = lhs->get_symb_id();
      if (symbol_table.getType(symb_id) != eEndogenous)
        continue;

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

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

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

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

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

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

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

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

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

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

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

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

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
574
  iterator_property_map<int *, property_map<AdjacencyList_t, vertex_index_t>::type, int, int &> endo2block_map(&endo2block[0], get(vertex_index, G2));
575
576

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

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

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

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

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

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

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

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

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

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

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

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

      Lag = Lead = 0;
830
      set<pair<int, int>> endo;
831
      for (count_equ  = first_count_equ; count_equ  < Blck_Size+first_count_equ; count_equ++)
832
833
834
        {
          endo.clear();
          equations[equation_reordered[count_equ]]->collectEndogenous(endo);
835
          for (const auto & it : endo)
836
            {
837
838
              int curr_variable = it.first;
              int curr_lag = it.second;
839
              auto it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
840
              if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
841
                if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end())
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
870
                  {
                    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;
        }
871
872
873
874
      l_n_static = n_static[i];
      l_n_forward = n_forward[i];
      l_n_backward = n_backward[i];
      l_n_mixed = n_mixed[i];
875
876
      if (Blck_Size == 1)
        {
877
          if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
878
879
880
881
882
883
884
885
            {
              if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
                Simulation_Type = EVALUATE_BACKWARD;
              else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
                Simulation_Type = EVALUATE_FORWARD;
            }
          if (i > 0)
            {
886
887
888
889
              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
890
                                 || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))
891
892
893
                {
                  for (int j = first_equation; j < first_equation+c_Size; j++)
                    {
894
                      auto it = dynamic_jacobian.find({ -1, { equation_reordered[eq], variable_reordered[j] } });
895
896
                      if (it != dynamic_jacobian.end())
                        is_lag = true;
897
                      it = dynamic_jacobian.find({ +1, { equation_reordered[eq], variable_reordered[j] } });
898
899
900
901
902
903
                      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))
904
905
906
                {
                  //merge the current block with the previous one
                  BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
907
                  c_Size++;
908
                  block_type_size_mfs[block_type_size_mfs.size()-1] = { { c_Type, first_equation }, { c_Size, c_Size } };
909
                  if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
910
                    Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
911
                  if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
912
                    Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
913
                  block_lag_lead[block_type_size_mfs.size()-1] = { Lag, Lead };
914
                  pair< pair< unsigned int, unsigned int>, pair<unsigned int, unsigned int>> tmp = block_col_type[block_col_type.size()-1];
915
                  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 } };
916
917
918
                }
              else
                {
919
920
                  block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
                  block_lag_lead.emplace_back(Lag, Lead);
921
                  block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
922
923
924
925
                }
            }
          else
            {
926
927
              block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
              block_lag_lead.emplace_back(Lag, Lead);
928
              block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
929
930
931
932
            }
        }
      else
        {
933
934
          block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
          block_lag_lead.emplace_back(Lag, Lead);
935
          block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
936
937
938
939
940
941
942
943
        }
      prev_Type = Simulation_Type;
      eq += Blck_Size;
    }
  return (block_type_size_mfs);
}

vector<bool>