ModelTree.cc 83.4 KB
Newer Older
1
/*
Houtan Bastani's avatar
Houtan Bastani committed
2
 * Copyright (C) 2003-2018 Dynare Team
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

20
#include <cstdlib>
21
#include <cassert>
sebastien's avatar
sebastien committed
22
#include <cmath>
23
#include <iostream>
24
#include <fstream>
25
26

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

using namespace MFS;

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

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

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

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

  // Fill in the graph
51
  set<pair<int, int>> endo;
52

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

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

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

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

67
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
68
69
70
71
72
73
74
75
76
77
78
79
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

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

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

  edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
  bool not_maximum_yet = true;
80
  while (not_maximum_yet)
81
82
83
84
85
86
87
88
89
90
91
    {
      not_maximum_yet = augmentor.augment_matching();
    }
  augmentor.get_current_matching(&mate_map[0]);

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

  assert(check);

#ifdef DEBUG
92
  for (int i = 0; i < n; i++)
93
94
95
96
97
    cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
         << " matched with equation " << (mate_map[i]-n+1) << endl;
#endif

  // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
98
  endo2eq.resize(equations.size());
99
100
101
102
103
104
105
106
  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));

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

  int n1 = 0, n2 = 0;

107
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

      n1++;

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

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

  // Check if all variables are normalized
126
  auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
127
  if (it != mate_map.begin() + n)
sebastien's avatar
sebastien committed
128
129
130
    {
      if (verbose)
        cerr << "ERROR: Could not normalize the model. Variable "
131
             << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
sebastien's avatar
sebastien committed
132
133
134
135
             << " is not in the maximum cardinality matching." << endl;
      check = false;
    }
  return check;
136
137
138
}

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

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

145
  int n = equations.size();
146

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

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

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

  int suppressed = 0;
  while (!check && current_cutoff > 1e-19)
163
    {
164
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
sebastien's avatar
sebastien committed
165
      int suppress = 0;
166
167
      for (auto & iter : normalized_contemporaneous_jacobian)
        if (fabs(iter.second) > max(current_cutoff, cutoff))
168
          tmp_normalized_contemporaneous_jacobian[{ iter.first.first, iter.first.second }] = iter.second;
sebastien's avatar
sebastien committed
169
170
171
172
173
174
175
        else
          suppress++;

      if (suppress != suppressed)
        check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
      suppressed = suppress;
      if (!check)
176
        {
sebastien's avatar
sebastien committed
177
178
179
180
          current_cutoff /= 2;
          // In this last case try to normalize with the complete jacobian
          if (current_cutoff <= 1e-19)
            check = computeNormalization(normalized_contemporaneous_jacobian, false);
181
182
183
        }
    }

sebastien's avatar
sebastien committed
184
  if (!check)
185
    {
sebastien's avatar
sebastien committed
186
187
      cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
      //if no non-singular normalization can be found, try to find a normalization even with a potential singularity
188
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
189
      set<pair<int, int>> endo;
sebastien's avatar
sebastien committed
190
      for (int i = 0; i < n; i++)
191
192
193
        {
          endo.clear();
          equations[i]->collectEndogenous(endo);
194
          for (const auto & it : endo)
195
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
196
        }
sebastien's avatar
sebastien committed
197
198
      check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
      if (check)
199
        {
sebastien's avatar
sebastien committed
200
          // Update the jacobian matrix
201
          for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
sebastien's avatar
sebastien committed
202
            {
203
204
205
206
207
208
              if (static_jacobian.find({ it->first.first, it->first.second }) == static_jacobian.end())
                static_jacobian[{ it->first.first, it->first.second }] = 0;
              if (dynamic_jacobian.find({ 0, { it->first.first, it->first.second } }) == dynamic_jacobian.end())
                dynamic_jacobian[{ 0, { it->first.first, it->first.second } }] = nullptr;
              if (contemporaneous_jacobian.find({ it->first.first, it->first.second }) == contemporaneous_jacobian.end())
                contemporaneous_jacobian[{ it->first.first, it->first.second }] = 0;
209
210
              try
                {
211
212
                  if (first_derivatives.find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == first_derivatives.end())
                    first_derivatives[{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
213
                }
Stéphane Adjemian's avatar
Stéphane Adjemian committed
214
              catch (DataTree::UnknownDerivIDException &e)
215
                {
216
                  cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it->first.second))
217
218
219
                       << " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
                  exit(EXIT_FAILURE);
                }
sebastien's avatar
sebastien committed
220
            }
221
222
        }
    }
sebastien's avatar
sebastien committed
223
224
225
226
227
228

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

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

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

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

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

void
255
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose)
256
257
{
  int nb_elements_contemparenous_Jacobian = 0;
258
  set<pair<int, int>> jacobian_elements_to_delete;
259
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
260
261
262
       it != first_derivatives.end(); it++)
    {
      int deriv_id = it->first.second;
263
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
264
        {
265
          expr_t Id = it->second;
266
267
268
269
270
271
272
273
274
          int eq = it->first.first;
          int symb = getSymbIDByDerivID(deriv_id);
          int var = symbol_table.getTypeSpecificID(symb);
          int lag = getLagByDerivID(deriv_id);
          double val = 0;
          try
            {
              val = Id->eval(eval_context);
            }
275
276
277
278
          catch (ExprNode::EvalExternalFunctionException &e)
            {
              val = 1;
            }
279
280
          catch (ExprNode::EvalException &e)
            {
281
              cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
282
              Id->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, temporary_terms, {});
283
284
285
286
287
288
289
              cerr << endl;
              exit(EXIT_FAILURE);
            }
          if (fabs(val) < cutoff)
            {
              if (verbose)
                cout << "the coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")" << endl;
290
              jacobian_elements_to_delete.emplace(eq, deriv_id);
291
292
293
294
295
296
            }
          else
            {
              if (lag == 0)
                {
                  nb_elements_contemparenous_Jacobian++;
297
                  contemporaneous_jacobian[{ eq, var }] = val;
298
                }
299
300
              if (static_jacobian.find({ eq, var }) != static_jacobian.end())
                static_jacobian[{ eq, var }] += val;
301
              else
302
303
                static_jacobian[{ eq, var }] = val;
              dynamic_jacobian[{ lag, { eq, var } }] = Id;
304
305
306
307
308
            }
        }
    }

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

312
  if (jacobian_elements_to_delete.size() > 0)
313
314
315
316
317
318
    {
      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;
    }
}

319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
vector<pair<int, int> >
ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered,
                                                     vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered,
                                                     lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead,
                                                     vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed)
{
  vector<int> eq2endo(equations.size(), 0);
  /*equation_reordered.resize(equations.size());
  variable_reordered.resize(equations.size());*/
  unsigned int num = 0;
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++)
    if (!is_equation_linear[*it])
      num++;
  vector<int> endo2block = vector<int>(endo2eq.size(), 1);
  vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
  int i = 0, j = 0;
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, j++)
    {
      if (!is_equation_linear[*it])
        {
          equation_reordered[i] = *it;
          variable_reordered[i] = j;
          endo2block[j] = 0;
          components_set[endo2block[j]].first.insert(i);
          /*cout << " -----------------------------------------" << endl;
          cout.flush();
          cout << "equation_reordered[" << *it << "] = " << i << endl;
          cout.flush();
          cout << "variable_reordered[" << j << "] = " << i << endl;
          cout.flush();
          cout << "components_set[" << endo2block[j] << "].first[" << i << "] = " << i << endl;
          cout << "endo2block[" << j << "] = " << 0 << endl;
          cout.flush();
          */
          i++;
        }
    }
/*  for (unsigned int j = 0; j < is_equation_linear.size() ; j++)
     cout << "endo2block[" << j << "] = " << endo2block[j] << endl;*/
  /*cout << "before getVariableLeadLAgByBlock\n";
  cout.flush();*/
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, endo2block.size(), equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
  n_static = vector<unsigned int>(endo2eq.size(), 0);
  n_forward = vector<unsigned int>(endo2eq.size(), 0);
  n_backward = vector<unsigned int>(endo2eq.size(), 0);
  n_mixed = vector<unsigned int>(endo2eq.size(), 0);
  for (unsigned int i = 0; i < endo2eq.size(); i++)
    {
      if      (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second != 0)
        n_mixed[i]++;
      else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second != 0)
        n_forward[i]++;
      else if (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second == 0)
        n_backward[i]++;
      else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second == 0)
        n_static[i]++;
    }
  cout.flush();
  int nb_endo = is_equation_linear.size();
  vector<pair<int, int> > blocks = vector<pair<int, int> >(1, make_pair(i, i));
  inv_equation_reordered = vector<int>(nb_endo);
  inv_variable_reordered = vector<int>(nb_endo);
  for (int i = 0; i < nb_endo; i++)
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
  return blocks;
}

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

418
void
419
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
420
{
421
422
423
  vector<int> eq2endo(equations.size(), 0);
  equation_reordered.resize(equations.size());
  variable_reordered.resize(equations.size());
424
  bool *IM;
425
  int n = equations.size();
426
  IM = (bool *) calloc(n*n, sizeof(bool));
427
  int i = 0;
428
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
429
430
431
432
    {
      eq2endo[*it] = i;
      equation_reordered[i] = i;
      variable_reordered[*it] = i;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
433
434
435
    }
  if (cutoff == 0)
    {
436
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
437
      for (int i = 0; i < n; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
438
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
439
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
440
          equations[i]->collectEndogenous(endo);
441
442
          for (const auto & it : endo)
            IM[i * n + endo2eq[it.first]] = true;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
443
444
445
        }
    }
  else
446
447
    for (const auto & it : static_jacobian_arg)
      IM[it.first.first * n + endo2eq[it.first.second]] = true;
448
449
450
451
452
453
454
455
  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;
456
      for (int i = prologue; i < n; i++)
457
458
        {
          int nze = 0;
459
460
          for (int j = tmp_prologue; j < n; j++)
            if (IM[i * n + j])
461
              {
462
                nze++;
463
464
                k = j;
              }
465
          if (nze == 1)
466
            {
467
              for (int j = 0; j < n; j++)
468
469
470
471
472
473
474
475
                {
                  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;
476
              for (int j = 0; j < n; j++)
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
                {
                  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;
499
      for (int i = prologue; i < n - (int) epilogue; i++)
500
501
        {
          int nze = 0;
502
503
          for (int j = prologue; j < n - tmp_epilogue; j++)
            if (IM[j * n + i])
504
              {
505
                nze++;
506
507
                k = j;
              }
508
          if (nze == 1)
509
            {
510
              for (int j = 0; j < n; j++)
511
512
513
514
515
516
517
518
                {
                  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;
519
              for (int j = 0; j < n; j++)
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
                {
                  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);
}

537
equation_type_and_normalized_equation_t
538
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
539
{
540
  expr_t lhs;
541
542
  BinaryOpNode *eq_node;
  EquationType Equation_Simulation_Type;
543
  equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
544
545
546
547
548
549
550
  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;
551
      auto derivative = first_order_endo_derivatives.find({ eq, { var, 0 } });
552
      pair<bool, expr_t> res;
553
      if (derivative != first_order_endo_derivatives.end())
554
        {
555
          set<pair<int, int>> result;
556
          derivative->second->collectEndogenous(result);
557
          auto d_endo_variable = result.find({ var, 0 });
558
          //Determine whether the equation could be evaluated rather than to be solved
559
          if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
560
561
562
563
564
            {
              Equation_Simulation_Type = E_EVALUATE;
            }
          else
            {
565
              vector<pair<int, pair<expr_t, expr_t>>> List_of_Op_RHS;
566
              res =  equations[eq]->normalizeEquation(var, List_of_Op_RHS);
567
              if (mfs == 2)
568
                {
569
                  if (d_endo_variable == result.end() && res.second)
570
571
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
572
              else if (mfs == 3)
573
                {
574
                  if (res.second) // The equation could be solved analytically
575
576
577
578
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
            }
        }
579
      V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
580
581
582
583
584
    }
  return (V_Equation_Simulation_Type);
}

void
585
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
586
587
{
  int nb_endo = symbol_table.endo_nbr();
588
589
  variable_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
  equation_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
590
591
592
  vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
  for (int i = 0; i < nb_endo; i++)
    {
593
      if (i < (int) prologue)
594
595
596
597
        {
          variable_blck[variable_reordered[i]] = i;
          equation_blck[equation_reordered[i]] = i;
        }
598
      else if (i < (int) (components_set.size() + prologue))
599
600
601
602
603
604
605
606
607
608
        {
          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);
        }
    }
609
  for (const auto & it : dynamic_jacobian)
610
    {
611
612
613
      int lag = it.first.first;
      int j_1 = it.first.second.first;
      int i_1 = it.first.second.second;
614
615
616
      if (variable_blck[i_1] == equation_blck[j_1])
        {
          if (lag > variable_lead_lag[i_1].second)
617
            variable_lead_lag[i_1] = { variable_lead_lag[i_1].first, lag };
618
          if (lag < -variable_lead_lag[i_1].first)
619
            variable_lead_lag[i_1] = { -lag, variable_lead_lag[i_1].second };
620
          if (lag > equation_lead_lag[j_1].second)
621
            equation_lead_lag[j_1] = { equation_lead_lag[j_1].first, lag };
622
          if (lag < -equation_lead_lag[j_1].first)
623
            equation_lead_lag[j_1] = { -lag, equation_lead_lag[j_1].second };
624
625
626
627
628
        }
    }
}

void
629
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
630
631
632
633
{
  int nb_var = variable_reordered.size();
  int n = nb_var - prologue - epilogue;

634
635
636
  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
637
  auto v_index = get(boost::vertex_index, G2);
638
639
  for (int i = 0; i < n; i++)
    put(v_index, vertex(i, G2), i);
640
641
642

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

643
  for (int i = 0; i < nb_var; i++)
644
645
646
    {
      reverse_equation_reordered[equation_reordered[i]] = i;
      reverse_variable_reordered[variable_reordered[i]] = i;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
647
648
649
650
    }
  jacob_map_t tmp_normalized_contemporaneous_jacobian;
  if (cutoff == 0)
    {
651
      set<pair<int, int>> endo;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
652
      for (int i = 0; i < nb_var; i++)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
653
        {
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
654
          endo.clear();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
655
          equations[i]->collectEndogenous(endo);
656
          for (const auto & it : endo)
657
            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
658
659
660
661
662

        }
    }
  else
    tmp_normalized_contemporaneous_jacobian = static_jacobian;
FerhatMihoubi's avatar
Fix :    
FerhatMihoubi committed
663
  for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
664
665
    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)
666
        && it->first.first != endo2eq[it->first.second])
667
668
669
      add_edge(vertex(reverse_equation_reordered[endo2eq[it->first.second]]-prologue, G2),
               vertex(reverse_equation_reordered[it->first.first]-prologue, G2),
               G2);
670
671

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
672
  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));
673
674

  // Compute strongly connected components
675
  int num = strong_components(G2, endo2block_map);
676

677
  blocks = vector<pair<int, int>>(num, { 0, 0 });
678
679

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

683
  for (unsigned int i = 0; i < num_vertices(G2); i++)
684
    {
685
686
      AdjacencyList_t::out_edge_iterator it_out, out_end;
      AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
687
688
      for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
        {
689
690
          int t_b = endo2block_map[target(*it_out, G2)];
          int s_b = endo2block_map[source(*it_out, G2)];
691
692
693
694
695
696
697
698
699
700
701
          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);
702
  for (int i = 0; i < num; i++)
703
704
705
706
707
708
    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.
709
  vector<pair<set<int>, pair<set<int>, vector<int>>>> components_set(num);
710
711
712
713
714
715
716
  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);
    }

717
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
718
719
720
721

  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
722
  if (select_feedback_variable)
723
724
725
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
726
727
728
729
730
            || 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)
731
          add_edge(vertex(i, G2), vertex(i, G2), G2);
732
733
734
735
736
    }
  else
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
737
          add_edge(vertex(i, G2), vertex(i, G2), G2);
738
    }
739
740
741
742
743
744
  //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);

745
  for (int i = 0; i < (int) prologue; i++)
746
747
748
749
750
751
752
753
754
755
    {
      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]++;
    }
756
757
758
759
760
761
  //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++)
    {
762
      AdjacencyList_t G = extract_subgraph(G2, components_set[i].first);
763
764
      set<int> feed_back_vertices;
      //Print(G);
765
      AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
766
      auto v_index = get(boost::vertex_index, G);
767
768
769
770
771
772
      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
773
      for (int j = 0; j < 4; j++)
774
        {
775
          for (int & its : Reordered_Vertice)
776
777
            {
              bool something_done = false;
778
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
779
780
781
782
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
783
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
784
785
786
787
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
788
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
789
790
791
792
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
793
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
794
795
796
797
798
799
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
800
801
                  equation_reordered[order] = tmp_equation_reordered[its+prologue];
                  variable_reordered[order] = tmp_variable_reordered[its+prologue];
802
803
804
                  order++;
                }
            }
805
806
807
        }
      components_set[i].second.second = Reordered_Vertice;
      //Second we have the equations related to the feedback variables
808
      for (int j = 0; j < 4; j++)
809
        {
810
          for (int feed_back_vertice : feed_back_vertices)
811
812
            {
              bool something_done = false;
813
              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)
814
815
816
817
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
818
              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)
819
820
821
822
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
823
              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)
824
825
826
827
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
828
              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)
829
830
831
832
833
834
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
835
836
                  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];
837
838
839
                  order++;
                }
            }
840
841
        }
    }
842

843
  for (int i = 0; i < (int) epilogue; i++)
844
    {
845
      if      (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
846
        n_mixed[prologue+num+i]++;
847
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
848
        n_forward[prologue+num+i]++;
849
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
850
        n_backward[prologue+num+i]++;
851
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
852
853
854
        n_static[prologue+num+i]++;
    }

855
856
  inv_equation_reordered = vector<int>(nb_var);
  inv_variable_reordered = vector<int>(nb_var);
857
  for (int i = 0; i < nb_var; i++)
858
859
860
861
862
863
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
}

864
void
865
ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
{
  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;
881
              Nb_feedback_variable = getBlockMfs(block);
882
883
884
885
886
887
888
889
890
891
892
            }
        }
    }

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

893
block_type_firstequation_size_mfs_t
894
ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int>>> &block_col_type, bool linear_decomposition)
895
896
897
898
899
{
  int i = 0;
  int count_equ = 0, blck_count_simult = 0;
  int Blck_Size, MFS_Size;
  int Lead, Lag;
900
  block_type_firstequation_size_mfs_t block_type_size_mfs;
901
902
  BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
  int eq = 0;
903
904
905
906
  unsigned int l_n_static = 0;
  unsigned int l_n_forward = 0;
  unsigned int l_n_backward = 0;
  unsigned int l_n_mixed = 0;
907
  for (i = 0; i < (int) (prologue+blocks.size()+epilogue); i++)
908
909
    {
      int first_count_equ = count_equ;
910
      if (i < (int) prologue)
911
912
913
914
        {
          Blck_Size = 1;
          MFS_Size = 1;
        }
915
      else if (i < (int) (prologue+blocks.size()))
916
917
918
919
920
        {
          Blck_Size = blocks[blck_count_simult].first;
          MFS_Size = blocks[blck_count_simult].second;
          blck_count_simult++;
        }
921
      else if (i < (int) (prologue+blocks.size()+epilogue))
922
923
924
925
926
927
        {
          Blck_Size = 1;
          MFS_Size = 1;
        }

      Lag = Lead = 0;
928
      set<pair<int, int>> endo;
929
      for (count_equ  = first_count_equ; count_equ  < Blck_Size+first_count_equ; count_equ++)
930
931
932
        {
          endo.clear();
          equations[equation_reordered[count_equ]]->collectEndogenous(endo);
933
          for (const auto & it : endo)
934
            {
935
936
              int curr_variable = it.first;
              int curr_lag = it.second;
937
              auto it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
              if (linear_decomposition)
                {
                  if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end())
                    {
                      if (curr_lag > Lead)
                        Lead = curr_lag;
                      else if (-curr_lag > Lag)
                        Lag = -curr_lag;
                    }
                }
              else
                {
                  if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
                    if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end())
                      {
                        if (curr_lag > Lead)
                          Lead = curr_lag;
                        else if (-curr_lag > Lag)
                          Lag = -curr_lag;
                      }
                }
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
            }
        }
      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;
        }
982
983
984
985
      l_n_static = n_static[i];
      l_n_forward = n_forward[i];
      l_n_backward = n_backward[i];
      l_n_mixed = n_mixed[i];
986
987
      if (Blck_Size == 1)
        {
988
          if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
989
990
991
992
993
994
995
996
            {
              if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
                Simulation_Type = EVALUATE_BACKWARD;
              else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
                Simulation_Type = EVALUATE_FORWARD;
            }
          if (i > 0)
            {
997
998
999
1000
              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
1001
                                 || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))