ModelTree.cc 85.3 KB
Newer Older
1
/*
2
 * Copyright (C) 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
28
29
30
31
32
33
34
#include "MinimumFeedbackSet.hh"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>

using namespace MFS;

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

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

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

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

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

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

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

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

100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
ModelTree::ModelTree(const ModelTree &m) :
  DataTree {m},
  equations_lineno {m.equations_lineno},
  equation_tags {m.equation_tags},
  NNZDerivatives {m.NNZDerivatives},
  equation_reordered {m.equation_reordered},
  variable_reordered {m.variable_reordered},
  inv_equation_reordered {m.inv_equation_reordered},
  inv_variable_reordered {m.inv_variable_reordered},
  is_equation_linear {m.is_equation_linear},
  endo2eq {m.endo2eq},
  epilogue {m.epilogue},
  prologue {m.prologue},
  block_lag_lead {m.block_lag_lead},
  cutoff {m.cutoff},
  mfs {m.mfs}
{
  copyHelper(m);
}

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

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

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

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

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

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

  copyHelper(m);

  return *this;
}


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

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

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

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

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

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

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

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

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

204
205
  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())
206
    {
207
208
    };

209
210
  augmentor.get_current_matching(&mate_map[0]);

211
  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));
212
213
214
215
216
#endif

  assert(check);

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

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

  int n1 = 0, n2 = 0;

232
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
233
234
235
236
237
238
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

      n1++;

239
240
241
      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))
242
243
244
245
246
247
248
249
250
             << " 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
251
  auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
252
  if (it != mate_map.begin() + n)
sebastien's avatar
sebastien committed
253
254
255
    {
      if (verbose)
        cerr << "ERROR: Could not normalize the model. Variable "
256
             << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
sebastien's avatar
sebastien committed
257
258
259
260
             << " is not in the maximum cardinality matching." << endl;
      check = false;
    }
  return check;
261
262
263
}

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

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

270
  int n = equations.size();
271

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

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

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

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

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

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

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

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

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

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

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

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

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

443
444
445
446
447
448
449
450
451
452
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;
453
454
  for (auto it : endo2eq)
    if (!is_equation_linear[it])
455
456
457
458
      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;
459
  for (auto it : endo2eq)
460
    {
461
      if (!is_equation_linear[it])
462
        {
463
          equation_reordered[i] = it;
464
465
466
467
468
469
470
471
472
473
474
475
476
477
          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++;
478
          j++;
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
        }
    }
/*  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();
503
504
505
  vector<pair<int, int>> blocks(1, make_pair(i, i));
  inv_equation_reordered.resize(nb_endo);
  inv_variable_reordered.resize(nb_endo);
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
  for (int i = 0; i < nb_endo; i++)
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
  return blocks;
}

bool
ModelTree::computeNaturalNormalization()
{
  bool bool_result = true;
  set<pair<int, int> > result;
  endo2eq.resize(equations.size());
  for (int eq = 0; eq < (int) equations.size(); eq++)
    if (!is_equation_linear[eq])
      {
        BinaryOpNode *eq_node = equations[eq];
524
        expr_t lhs = eq_node->arg1;
525
526
527
528
529
        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 !
530
             auto it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
531
532
533
534
535
536
537
538
539
540
541
542
             if (it == endo2eq.end())
               endo2eq[result.begin()->first] = eq;
             else
              {
                bool_result = false;
                break;
              }
          }
      }
  return bool_result;
}

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

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

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

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

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

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

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

        }
    }
  else
    tmp_normalized_contemporaneous_jacobian = static_jacobian;
784
785
786
787
788
789
  for (const auto &it : tmp_normalized_contemporaneous_jacobian)
    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)
        && 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),
790
               G2);
791
792

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
793
  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));
794
795

  // Compute strongly connected components
796
  int num = strong_components(G2, endo2block_map);
797

798
  blocks = vector<pair<int, int>>(num, { 0, 0 });
799
800

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

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

838
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
839
840
841
842

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

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

964
  for (int i = 0; i < (int) epilogue; i++)
965
    {
966
      if      (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
967
        n_mixed[prologue+num+i]++;
968
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_le