ModelTree.cc 85.7 KB
Newer Older
1
/*
2
 * Copyright © 2003-2019 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
#include "MinimumFeedbackSet.hh"
28
29
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wold-style-cast"
30
31
32
33
#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>
34
#pragma GCC diagnostic pop
35
36
37

using namespace MFS;

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

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

49
50
51
52
53
54
55
56
  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;
    };

57
  // Derivatives
58
59
60
61
62
63
64
65
66
67
68
69
  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;
    };
70
71
72
73
74
75

  // 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);
76
77
  for (const auto &it : m.temporary_terms_derivatives)
    temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
78
79
80
  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)
81
    params_derivs_temporary_terms[it.first] = convert_temporary_terms_t(it.second);
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)
89
    nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)};
90
91
}

92
93
94
95
96
97
98
99
100
101
102
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)
{
}

103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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;

134
135
  derivatives.clear();
  params_derivatives.clear();
136
137
138

  temporary_terms.clear();
  temporary_terms_mlv.clear();
139
  temporary_terms_derivatives.clear();
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
  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
164
bool
165
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
166
{
167
  const int n = equations.size();
168
169
170

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

171
  using BipartiteGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
172
173
174
175
176
177
178
179

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

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

  // 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
191
  fill(mate_map.begin(), mate_map.end(), boost::graph_traits<BipartiteGraph>::null_vertex());
192
193
194
195

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

196
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
197
198
199
200
201
202
203
204
205
206
    {
      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);
    }

207
208
  boost::edmonds_augmenting_path_finder<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type> augmentor(g, &mate_map[0], get(boost::vertex_index, g));
  while (augmentor.augment_matching())
209
    {
210
211
    };

212
213
  augmentor.get_current_matching(&mate_map[0]);

214
  bool check = boost::maximum_cardinality_matching_verifier<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(boost::vertex_index, g));
215
216
217
218
219
#endif

  assert(check);

#ifdef DEBUG
220
  for (int i = 0; i < n; i++)
221
222
223
224
225
    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
226
  endo2eq.resize(equations.size());
227
  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](int i) { return i-n; });
228
229
230
231
232
233
234

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

  int n1 = 0, n2 = 0;

235
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
236
237
238
239
240
241
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

      n1++;

242
243
244
      auto x = natural_endo2eqs.equal_range(i);
      if (find_if(x.first, x.second, [=](auto y) { return y.second == endo2eq[i]; }) == x.second)
        cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, i))
245
246
247
248
249
250
251
252
253
             << " 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
254
  auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
255
  if (it != mate_map.begin() + n)
sebastien's avatar
sebastien committed
256
257
258
    {
      if (verbose)
        cerr << "ERROR: Could not normalize the model. Variable "
259
             << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
sebastien's avatar
sebastien committed
260
261
262
263
             << " is not in the maximum cardinality matching." << endl;
      check = false;
    }
  return check;
264
265
266
}

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

271
272
  cout << "Normalizing the model..." << endl;

273
  int n = equations.size();
274

sebastien's avatar
sebastien committed
275
276
  // compute the maximum value of each row of the contemporaneous Jacobian matrix
  //jacob_map normalized_contemporaneous_jacobian;
277
  jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
sebastien's avatar
sebastien committed
278
  vector<double> max_val(n, 0.0);
279
280
281
  for (const auto &it : contemporaneous_jacobian)
    if (fabs(it.second) > max_val[it.first.first])
      max_val[it.first.first] = fabs(it.second);
282

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

  //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)
291
    {
292
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
sebastien's avatar
sebastien committed
293
      int suppress = 0;
294
295
      for (auto & iter : normalized_contemporaneous_jacobian)
        if (fabs(iter.second) > max(current_cutoff, cutoff))
296
          tmp_normalized_contemporaneous_jacobian[{ iter.first.first, iter.first.second }] = iter.second;
sebastien's avatar
sebastien committed
297
298
299
300
301
302
303
        else
          suppress++;

      if (suppress != suppressed)
        check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
      suppressed = suppress;
      if (!check)
304
        {
sebastien's avatar
sebastien committed
305
306
307
308
          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);
309
310
311
        }
    }

sebastien's avatar
sebastien committed
312
  if (!check)
313
    {
sebastien's avatar
sebastien committed
314
315
      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
316
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
317
      set<pair<int, int>> endo;
sebastien's avatar
sebastien committed
318
      for (int i = 0; i < n; i++)
319
320
321
        {
          endo.clear();
          equations[i]->collectEndogenous(endo);
322
          for (const auto & it : endo)
323
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
324
        }
sebastien's avatar
sebastien committed
325
326
      check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
      if (check)
327
        {
sebastien's avatar
sebastien committed
328
          // Update the jacobian matrix
329
          for (const auto &it : tmp_normalized_contemporaneous_jacobian)
sebastien's avatar
sebastien committed
330
            {
331
332
333
334
335
336
              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;
337
338
              try
                {
339
340
                  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;
341
                }
Stéphane Adjemian's avatar
Stéphane Adjemian committed
342
              catch (DataTree::UnknownDerivIDException &e)
343
                {
344
                  cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it.first.second))
345
346
347
                       << " 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
348
            }
349
350
        }
    }
sebastien's avatar
sebastien committed
351
352
353
354
355
356

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

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

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

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

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

void
383
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)
384
385
{
  int nb_elements_contemparenous_Jacobian = 0;
386
  set<vector<int>> jacobian_elements_to_delete;
387
  for (const auto &it : derivatives[1])
388
    {
389
      int deriv_id = it.first[1];
390
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
391
        {
392
393
          expr_t Id = it.second;
          int eq = it.first[0];
394
395
396
397
398
399
400
401
          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);
            }
402
403
404
405
          catch (ExprNode::EvalExternalFunctionException &e)
            {
              val = 1;
            }
406
407
          catch (ExprNode::EvalException &e)
            {
408
              cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
409
              Id->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, temporary_terms, {});
410
411
412
413
414
415
416
              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;
417
              jacobian_elements_to_delete.insert({ eq, deriv_id });
418
419
420
421
422
423
            }
          else
            {
              if (lag == 0)
                {
                  nb_elements_contemparenous_Jacobian++;
424
                  contemporaneous_jacobian[{ eq, var }] = val;
425
                }
426
427
              if (static_jacobian.find({ eq, var }) != static_jacobian.end())
                static_jacobian[{ eq, var }] += val;
428
              else
429
                static_jacobian[{ eq, var }] = val;
430
              dynamic_jacobian[{ lag, eq, var }] = Id;
431
432
433
434
435
            }
        }
    }

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

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

446
447
448
449
450
451
452
453
454
455
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;
456
457
  for (auto it : endo2eq)
    if (!is_equation_linear[it])
458
459
460
461
      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;
462
  for (auto it : endo2eq)
463
    {
464
      if (!is_equation_linear[it])
465
        {
466
          equation_reordered[i] = it;
467
468
469
470
471
472
473
474
475
476
477
478
479
480
          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++;
481
          j++;
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
        }
    }
/*  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();
506
  vector<pair<int, int>> blocks(1, {i, i});
507
508
  inv_equation_reordered.resize(nb_endo);
  inv_variable_reordered.resize(nb_endo);
509
510
511
512
513
514
515
516
517
518
519
520
521
522
  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());
523
  for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
524
525
526
    if (!is_equation_linear[eq])
      {
        BinaryOpNode *eq_node = equations[eq];
527
        expr_t lhs = eq_node->arg1;
528
529
530
531
532
        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 !
533
             auto it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
534
535
536
537
538
539
540
541
542
543
544
545
             if (it == endo2eq.end())
               endo2eq[result.begin()->first] = eq;
             else
              {
                bool_result = false;
                break;
              }
          }
      }
  return bool_result;
}

546
void
547
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
548
{
549
550
551
552
  vector<int> eq2endo(equations.size(), 0);
  equation_reordered.resize(equations.size());
  variable_reordered.resize(equations.size());
  int n = equations.size();
553
  vector<bool> IM(n*n);
554
  int i = 0;
555
  for (auto it : endo2eq)
556
    {
557
      eq2endo[it] = i;
558
      equation_reordered[i] = i;
559
560
      variable_reordered[it] = i;
      i++;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
561
562
563
    }
  if (cutoff == 0)
    {
564
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
565
      for (int i = 0; i < n; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
566
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
567
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
568
          equations[i]->collectEndogenous(endo);
569
570
          for (const auto & it : endo)
            IM[i * n + endo2eq[it.first]] = true;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
571
572
573
        }
    }
  else
574
575
    for (const auto & it : static_jacobian_arg)
      IM[it.first.first * n + endo2eq[it.first.second]] = true;
576
577
578
579
580
581
582
583
  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;
584
      for (int i = prologue; i < n; i++)
585
586
        {
          int nze = 0;
587
588
          for (int j = tmp_prologue; j < n; j++)
            if (IM[i * n + j])
589
              {
590
                nze++;
591
592
                k = j;
              }
593
          if (nze == 1)
594
            {
595
              for (int j = 0; j < n; j++)
596
597
598
599
600
601
602
603
                {
                  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;
604
              for (int j = 0; j < n; j++)
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
                {
                  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;
627
      for (int i = prologue; i < n - static_cast<int>(epilogue); i++)
628
629
        {
          int nze = 0;
630
631
          for (int j = prologue; j < n - tmp_epilogue; j++)
            if (IM[j * n + i])
632
              {
633
                nze++;
634
635
                k = j;
              }
636
          if (nze == 1)
637
            {
638
              for (int j = 0; j < n; j++)
639
640
641
642
643
644
645
646
                {
                  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;
647
              for (int j = 0; j < n; j++)
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
                {
                  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;
    }
}

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

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

void
753
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
754
755
756
757
{
  int nb_var = variable_reordered.size();
  int n = nb_var - prologue - epilogue;

758
759
760
  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
761
  auto v_index = get(boost::vertex_index, G2);
762
763
  for (int i = 0; i < n; i++)
    put(v_index, vertex(i, G2), i);
764
765
766

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

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

        }
    }
  else
    tmp_normalized_contemporaneous_jacobian = static_jacobian;
787
  for (const auto &it : tmp_normalized_contemporaneous_jacobian)
788
789
    if (reverse_equation_reordered[it.first.first] >= static_cast<int>(prologue) && reverse_equation_reordered[it.first.first] < static_cast<int>(nb_var - epilogue)
        && reverse_variable_reordered[it.first.second] >= static_cast<int>(prologue) && reverse_variable_reordered[it.first.second] < static_cast<int>(nb_var - epilogue)
790
791
792
        && it.first.first != endo2eq[it.first.second])
      add_edge(vertex(reverse_equation_reordered[endo2eq[it.first.second]]-prologue, G2),
               vertex(reverse_equation_reordered[it.first.first]-prologue, G2),
793
               G2);
794
795

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
796
  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));
797
798

  // Compute strongly connected components
799
  int num = strong_components(G2, endo2block_map);
800

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

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

807
  for (unsigned int i = 0; i < num_vertices(G2); i++)
808
    {
809
810
      AdjacencyList_t::out_edge_iterator it_out, out_end;
      AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
811
812
      for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
        {
813
814
          int t_b = endo2block_map[target(*it_out, G2)];
          int s_b = endo2block_map[source(*it_out, G2)];
815
816
817
818
819
820
821
822
823
824
825
          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);
826
  for (int i = 0; i < num; i++)
827
828
829
830
831
832
    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.
833
  vector<tuple<set<int>, set<int>, vector<int>>> components_set(num);
834
835
836
837
  for (unsigned int i = 0; i < endo2block.size(); i++)
    {
      endo2block[i] = unordered2ordered[endo2block[i]];
      blocks[endo2block[i]].first++;
838
      get<0>(components_set[endo2block[i]]).insert(i);
839
840
    }

841
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
842
843
844
845

  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
846
  if (select_feedback_variable)
847
848
849
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
850
851
852
853
854
            || 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)
855
          add_edge(vertex(i, G2), vertex(i, G2), G2);
856
857
858
859
860
    }
  else
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
861
          add_edge(vertex(i, G2), vertex(i, G2), G2);
862
    }
863
864
865
866
867
868
  //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);

869
  for (int i = 0; i < static_cast<int>(prologue); i++)
870
871
872
873
874
875
876
877
878
879
    {
      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]++;
    }
880
881
882
883
884
885
  //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++)
    {
886
      AdjacencyList_t G = extract_subgraph(G2, get<0>(components_set[i]));
887
888
      set<int> feed_back_vertices;
      //Print(G);
889
      AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
890
      auto v_index = get(boost::vertex_index, G);
891
      get<1>(components_set[i]) = feed_back_vertices;
892
893
894
895
896
      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
897
      for (int j = 0; j < 4; j++)
898
        {
899
          for (int its : Reordered_Vertice)
900
901
            {
              bool something_done = false;
902
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
903
904
905
906
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
907
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
908
909
910
911
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
912
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
913
914
915
916
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
917
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
918
919
920
921
922
923
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
924
925
                  equation_reordered[order] = tmp_equation_reordered[its+prologue];
                  variable_reordered[order] = tmp_variable_reordered[its+prologue];
926
927
928
                  order++;
                }
            }
929
        }
930
      get<2>(components_set[i]) = Reordered_Vertice;
931
      //Second we have the equations related to the feedback variables
932
      for (int j = 0; j < 4; j++)
933
        {
934
          for (int feed_back_vertice : feed_back_vertices)
935
936
            {
              bool something_done = false;
937
              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)
938
939
940
941
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
942
              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)
943
944
945
946
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
947
              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)
948
949
950
951
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
952
              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)
953
954
955
956
957
958
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {