ModelTree.cc 55.4 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
sebastien's avatar
sebastien committed
2
 * Copyright (C) 2003-2010 Dynare Team
sebastien's avatar
sebastien committed
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>
sebastien's avatar
sebastien committed
25

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

using namespace boost;
using namespace MFS;

sebastien's avatar
sebastien committed
36
bool
37
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
{
  const int n = equation_number();

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

  typedef adjacency_list<vecS, vecS, undirectedS> BipartiteGraph;

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

54
  for (jacob_map_t::const_iterator it = contemporaneous_jacobian.begin(); it != contemporaneous_jacobian.end(); it++)
sebastien's avatar
sebastien committed
55
    add_edge(it->first.first + n, it->first.second, g);
56 57 58 59 60 61 62 63 64 65 66 67

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

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

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

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

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

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

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

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

  assert(check);

#ifdef DEBUG
93
  for (int i = 0; i < n; i++)
94 95 96 97 98 99 100 101 102 103 104 105 106 107
    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
  endo2eq.resize(equation_number());
  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));

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

  int n1 = 0, n2 = 0;

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

      n1++;

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

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

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

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

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

sebastien's avatar
sebastien committed
146
  int n = equation_number();
147

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

156
  for (jacob_map_t::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
sebastien's avatar
sebastien committed
157 158 159 160 161 162 163
    iter->second /= max_val[iter->first.first];

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

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

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

sebastien's avatar
sebastien committed
185
  if (!check)
186
    {
sebastien's avatar
sebastien committed
187 188
      cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
      //if no non-singular normalization can be found, try to find a normalization even with a potential singularity
189
      jacob_map_t tmp_normalized_contemporaneous_jacobian;
190
      set<pair<int, int> > endo;
sebastien's avatar
sebastien committed
191
      for (int i = 0; i < n; i++)
192 193 194
        {
          endo.clear();
          equations[i]->collectEndogenous(endo);
195
          for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
sebastien's avatar
sebastien committed
196
            tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1;
197
        }
sebastien's avatar
sebastien committed
198 199
      check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
      if (check)
200
        {
sebastien's avatar
sebastien committed
201
          // Update the jacobian matrix
202
          for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
sebastien's avatar
sebastien committed
203 204 205 206 207 208 209
            {
              if (static_jacobian.find(make_pair(it->first.first, it->first.second)) == static_jacobian.end())
                static_jacobian[make_pair(it->first.first, it->first.second)] = 0;
              if (dynamic_jacobian.find(make_pair(0, make_pair(it->first.first, it->first.second))) == dynamic_jacobian.end())
                dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
              if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
                contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
210 211
              if (first_derivatives.find(make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))) == first_derivatives.end())
                first_derivatives[make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))] = Zero;
sebastien's avatar
sebastien committed
212
            }
213 214
        }
    }
sebastien's avatar
sebastien committed
215 216 217 218 219 220

  if (!check)
    {
      cerr << "No normalization could be computed. Aborting." << endl;
      exit(EXIT_FAILURE);
    }
221 222 223 224 225
}

void
ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
{
226
  for (int i = 0; i < equation_number(); i++)
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
    {
      VariableNode *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
      if (lhs == NULL)
        continue;

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

      set<pair<int, int> > endo;
      equations[i]->get_arg2()->collectEndogenous(endo);
      if (endo.find(make_pair(symbol_table.getTypeSpecificID(symb_id), 0)) != endo.end())
        continue;

      endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
      cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
    }
}

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

  // Get rid of the elements of the Jacobian matrix below the cutoff
301
  for (set<pair<int, int> >::const_iterator it = jacobian_elements_to_delete.begin(); it != jacobian_elements_to_delete.end(); it++)
302 303
    first_derivatives.erase(*it);

304
  if (jacobian_elements_to_delete.size() > 0)
305 306 307 308 309 310 311
    {
      cout << jacobian_elements_to_delete.size() << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
           << "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
    }
}

void
312
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
313
{
314
  vector<int> eq2endo(equation_number(), 0);
315 316 317 318
  equation_reordered.resize(equation_number());
  variable_reordered.resize(equation_number());
  bool *IM;
  int n = equation_number();
319
  IM = (bool *) calloc(n*n, sizeof(bool));
320
  int i = 0;
321
  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
322 323 324 325 326
    {
      eq2endo[*it] = i;
      equation_reordered[i] = i;
      variable_reordered[*it] = i;
    }
327
  for (jacob_map_t::const_iterator it = static_jacobian_arg.begin(); it != static_jacobian_arg.end(); it++)
328 329 330 331 332 333 334 335 336
    IM[it->first.first * n + endo2eq[it->first.second]] = true;
  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;
337
      for (int i = prologue; i < n; i++)
338 339
        {
          int nze = 0;
340 341
          for (int j = tmp_prologue; j < n; j++)
            if (IM[i * n + j])
342
              {
343
                nze++;
344 345
                k = j;
              }
346
          if (nze == 1)
347
            {
348
              for (int j = 0; j < n; j++)
349 350 351 352 353 354 355 356
                {
                  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;
357
              for (int j = 0; j < n; j++)
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379
                {
                  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;
380
      for (int i = prologue; i < n - (int) epilogue; i++)
381 382
        {
          int nze = 0;
383 384
          for (int j = prologue; j < n - tmp_epilogue; j++)
            if (IM[j * n + i])
385
              {
386
                nze++;
387 388
                k = j;
              }
389
          if (nze == 1)
390
            {
391
              for (int j = 0; j < n; j++)
392 393 394 395 396 397 398 399
                {
                  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;
400
              for (int j = 0; j < n; j++)
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417
                {
                  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);
}

418
equation_type_and_normalized_equation_t
419
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
420
{
421
  expr_t lhs, rhs;
422 423
  BinaryOpNode *eq_node;
  EquationType Equation_Simulation_Type;
424
  equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
425 426 427 428 429 430 431 432
  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();
      rhs = eq_node->get_arg2();
      Equation_Simulation_Type = E_SOLVE;
433 434
      map<pair<int, pair<int, int> >, expr_t>::const_iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
      pair<bool, expr_t> res;
435
      if (derivative != first_order_endo_derivatives.end())
436 437 438 439 440
        {
          set<pair<int, int> > result;
          derivative->second->collectEndogenous(result);
          set<pair<int, int> >::const_iterator d_endo_variable = result.find(make_pair(var, 0));
          //Determine whether the equation could be evaluated rather than to be solved
441
          if (lhs->isVariableNodeEqualTo(eEndogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
442 443 444 445 446
            {
              Equation_Simulation_Type = E_EVALUATE;
            }
          else
            {
447
              vector<pair<int, pair<expr_t, expr_t> > > List_of_Op_RHS;
448
              res =  equations[eq]->normalizeEquation(var, List_of_Op_RHS);
449
              if (mfs == 2)
450
                {
451
                  if (d_endo_variable == result.end() && res.second)
452 453
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
454
              else if (mfs == 3)
455
                {
456
                  if (res.second) // The equation could be solved analytically
457 458 459 460 461 462 463 464 465 466
                    Equation_Simulation_Type = E_EVALUATE_S;
                }
            }
        }
      V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second));
    }
  return (V_Equation_Simulation_Type);
}

void
467
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
468 469
{
  int nb_endo = symbol_table.endo_nbr();
470 471
  variable_lead_lag = lag_lead_vector_t(nb_endo, make_pair(0, 0));
  equation_lead_lag = lag_lead_vector_t(nb_endo, make_pair(0, 0));
472 473 474 475 476 477 478 479
  vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
  for (int i = 0; i < nb_endo; i++)
    {
      if (i < prologue)
        {
          variable_blck[variable_reordered[i]] = i;
          equation_blck[equation_reordered[i]] = i;
        }
480
      else if (i < (int) components_set.size() + prologue)
481 482 483 484 485 486 487 488 489 490
        {
          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);
        }
    }
491
  for (dynamic_jacob_map_t::const_iterator it = dynamic_jacobian.begin(); it != dynamic_jacobian.end(); it++)
492 493
    {
      int lag = it->first.first;
494
      int j_1 = it->first.second.first;
495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
      int i_1 = it->first.second.second;
      if (variable_blck[i_1] == equation_blck[j_1])
        {
          if (lag > variable_lead_lag[i_1].second)
            variable_lead_lag[i_1] = make_pair(variable_lead_lag[i_1].first, lag);
          if (lag < -variable_lead_lag[i_1].first)
            variable_lead_lag[i_1] = make_pair(-lag, variable_lead_lag[i_1].second);
          if (lag > equation_lead_lag[j_1].second)
            equation_lead_lag[j_1] = make_pair(equation_lead_lag[j_1].first, lag);
          if (lag < -equation_lead_lag[j_1].first)
            equation_lead_lag[j_1] = make_pair(-lag, equation_lead_lag[j_1].second);
        }
    }
}

void
511
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
512 513 514 515 516 517 518 519 520
{
  int nb_var = variable_reordered.size();
  int n = nb_var - prologue - epilogue;
  typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;

  GraphvizDigraph G2(n);

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

521
  for (int i = 0; i < nb_var; i++)
522 523 524 525 526
    {
      reverse_equation_reordered[equation_reordered[i]] = i;
      reverse_variable_reordered[variable_reordered[i]] = i;
    }

527
  for (jacob_map_t::const_iterator it = static_jacobian.begin(); it != static_jacobian.end(); it++)
528 529 530
    if (reverse_equation_reordered[it->first.first] >= prologue && reverse_equation_reordered[it->first.first] < nb_var - epilogue
        && reverse_variable_reordered[it->first.second] >= prologue && reverse_variable_reordered[it->first.second] < nb_var - epilogue
        && it->first.first != endo2eq[it->first.second])
531
      add_edge(reverse_equation_reordered[endo2eq[it->first.second]]-prologue, reverse_equation_reordered[it->first.first]-prologue, G2);
532 533 534 535 536 537 538 539 540 541 542 543

  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));

  // Compute strongly connected components
  int num = strong_components(G2, &endo2block[0]);

  blocks = vector<pair<int, int> >(num, make_pair(0, 0));

  // Create directed acyclic graph associated to the strongly connected components
  typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
  DirectedGraph dag(num);

544
  for (unsigned int i = 0; i < num_vertices(G2); i++)
545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562
    {
      GraphvizDigraph::out_edge_iterator it_out, out_end;
      GraphvizDigraph::vertex_descriptor vi = vertex(i, G2);
      for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
        {
          int t_b = endo2block[target(*it_out, G2)];
          int s_b = endo2block[source(*it_out, G2)];
          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);
563
  for (int i = 0; i < num; i++)
564 565 566 567 568 569 570 571 572 573 574 575 576 577
    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.
  vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
  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);
    }

578
  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
579 580 581 582

  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
583
  if (select_feedback_variable)
584 585 586
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
587 588 589 590 591
            || 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)
592 593 594 595 596 597 598 599
          add_edge(i, i, G2);
    }
  else
    {
      for (int i = 0; i < n; i++)
        if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
          add_edge(i, i, G2);
    }
600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616
  //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);

  for (int i = 0; i < prologue; i++)
    {
      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]++;
    }
617 618 619 620 621 622
  //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++)
    {
623
      AdjacencyList_t G = GraphvizDigraph_2_AdjacencyList(G2, components_set[i].first);
624 625
      set<int> feed_back_vertices;
      //Print(G);
626 627
      AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
      property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
628 629 630 631 632 633
      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
634
      for (int j = 0; j < 4; j++)
635
        {
636 637 638
          for (vector<int>::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
            {
              bool something_done = false;
639
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
640 641 642 643
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
644
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
645 646 647 648
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
649
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
650 651 652 653
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
654
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
655 656 657 658 659 660 661 662 663 664 665
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
                  equation_reordered[order] = tmp_equation_reordered[*its+prologue];
                  variable_reordered[order] = tmp_variable_reordered[*its+prologue];
                  order++;
                }
            }
666 667 668
        }
      components_set[i].second.second = Reordered_Vertice;
      //Second we have the equations related to the feedback variables
669
      for (int j = 0; j < 4; j++)
670
        {
671 672 673
          for (set<int>::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
            {
              bool something_done = false;
674
              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
675 676 677 678
                {
                  n_mixed[prologue+i]++;
                  something_done = true;
                }
679
              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
680 681 682 683
                {
                  n_forward[prologue+i]++;
                  something_done = true;
                }
684
              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
685 686 687 688
                {
                  n_backward[prologue+i]++;
                  something_done = true;
                }
689
              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
690 691 692 693 694 695 696 697 698 699 700
                {
                  n_static[prologue+i]++;
                  something_done = true;
                }
              if (something_done)
                {
                  equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
                  variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
                  order++;
                }
            }
701 702
        }
    }
703 704 705

  for (int i = 0; i < epilogue; i++)
    {
706
      if      (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
707
        n_mixed[prologue+num+i]++;
708
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
709
        n_forward[prologue+num+i]++;
710
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
711
        n_backward[prologue+num+i]++;
712
      else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
713 714 715
        n_static[prologue+num+i]++;
    }

716 717
  inv_equation_reordered = vector<int>(nb_var);
  inv_variable_reordered = vector<int>(nb_var);
718
  for (int i = 0; i < nb_var; i++)
719 720 721 722 723 724
    {
      inv_variable_reordered[variable_reordered[i]] = i;
      inv_equation_reordered[equation_reordered[i]] = i;
    }
}

725
void
726
ModelTree::printBlockDecomposition(const vector<pair<int, int> > &blocks) const
727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753
{
  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;
              Nb_feedback_variable = blocks[Nb_SimulBlocks-1].second;
            }
        }
    }

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

754
block_type_firstequation_size_mfs_t
755
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)
756 757 758 759 760
{
  int i = 0;
  int count_equ = 0, blck_count_simult = 0;
  int Blck_Size, MFS_Size;
  int Lead, Lag;
761
  block_type_firstequation_size_mfs_t block_type_size_mfs;
762 763
  BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
  int eq = 0;
764 765 766 767
  unsigned int l_n_static = 0;
  unsigned int l_n_forward = 0;
  unsigned int l_n_backward = 0;
  unsigned int l_n_mixed = 0;
768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789
  for (i = 0; i < prologue+(int) blocks.size()+epilogue; i++)
    {
      int first_count_equ = count_equ;
      if (i < prologue)
        {
          Blck_Size = 1;
          MFS_Size = 1;
        }
      else if (i < prologue+(int) blocks.size())
        {
          Blck_Size = blocks[blck_count_simult].first;
          MFS_Size = blocks[blck_count_simult].second;
          blck_count_simult++;
        }
      else if (i < prologue+(int) blocks.size()+epilogue)
        {
          Blck_Size = 1;
          MFS_Size = 1;
        }

      Lag = Lead = 0;
      set<pair<int, int> > endo;
790
      for (count_equ  = first_count_equ; count_equ  < Blck_Size+first_count_equ; count_equ++)
791 792 793 794 795 796 797 798
        {
          endo.clear();
          equations[equation_reordered[count_equ]]->collectEndogenous(endo);
          for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
            {
              int curr_variable = it->first;
              int curr_lag = it->second;
              vector<int>::const_iterator it = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
799
              if (it != variable_reordered.begin()+(first_count_equ+Blck_Size))
800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829
                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;
                  }
            }
        }
      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;
        }
830 831 832 833
      l_n_static = n_static[i];
      l_n_forward = n_forward[i];
      l_n_backward = n_backward[i];
      l_n_mixed = n_mixed[i];
834 835
      if (Blck_Size == 1)
        {
836
          if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
837 838 839 840 841 842 843 844
            {
              if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
                Simulation_Type = EVALUATE_BACKWARD;
              else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
                Simulation_Type = EVALUATE_FORWARD;
            }
          if (i > 0)
            {
845 846
              if ((prev_Type ==  EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD)
                  || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD))
847 848 849 850 851
                {
                  //merge the current block with the previous one
                  BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
                  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;
852 853
                  c_Size++;
                  block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(c_Size, c_Size));
854
                  if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
855
                    Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
856
                  if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
857 858
                    Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
                  block_lag_lead[block_type_size_mfs.size()-1] = make_pair(Lag, Lead);
859 860
                  pair< pair< unsigned int, unsigned int>, pair<unsigned int, unsigned int> > tmp = block_col_type[block_col_type.size()-1];
                  block_col_type[block_col_type.size()-1] = make_pair( make_pair(tmp.first.first+l_n_static, tmp.first.second+l_n_forward), make_pair(tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed) );
861 862 863 864 865
                }
              else
                {
                  block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
                  block_lag_lead.push_back(make_pair(Lag, Lead));
866
                  block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
867 868 869 870 871 872
                }
            }
          else
            {
              block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
              block_lag_lead.push_back(make_pair(Lag, Lead));
873
              block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
874 875 876 877 878 879
            }
        }
      else
        {
          block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
          block_lag_lead.push_back(make_pair(Lag, Lead));
880
          block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
881 882 883 884 885 886 887 888
        }
      prev_Type = Simulation_Type;
      eq += Blck_Size;
    }
  return (block_type_size_mfs);
}

vector<bool>
889
ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector<int> &variable_reordered) const
890 891 892
{
  unsigned int nb_blocks = getNbBlocks();
  vector<bool> blocks_linear(nb_blocks, true);
893
  for (unsigned int block = 0; block < nb_blocks; block++)
894 895 896
    {
      BlockSimulationType simulation_type = getBlockSimulationType(block);
      int block_size = getBlockSize(block);
897
      block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block];
898
      int first_variable_position = getBlockFirstEquation(block);
899
      if (simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
900
        {
901
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
902 903 904 905
            {
              int lag = it->second.first;
              if (lag == 0)
                {
906
                  expr_t Id = it->second.second;
907 908 909 910
                  set<pair<int, int> > endogenous;
                  Id->collectEndogenous(endogenous);
                  if (endogenous.size() > 0)
                    {
911
                      for (int l = 0; l < block_size; l++)
912 913 914 915 916 917 918 919 920 921 922
                        {
                          if (endogenous.find(make_pair(variable_reordered[first_variable_position+l], 0)) != endogenous.end())
                            {
                              blocks_linear[block] = false;
                              goto the_end;
                            }
                        }
                    }
                }
            }
        }
923
      else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
924
        {
925
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
926 927
            {
              int lag = it->second.first;
928
              expr_t Id = it->second.second; //
929 930 931 932
              set<pair<int, int> > endogenous;
              Id->collectEndogenous(endogenous);
              if (endogenous.size() > 0)
                {
933
                  for (int l = 0; l < block_size; l++)
934 935 936 937 938 939 940 941 942 943
                    {
                      if (endogenous.find(make_pair(variable_reordered[first_variable_position+l], lag)) != endogenous.end())
                        {
                          blocks_linear[block] = false;
                          goto the_end;
                        }
                    }
                }
            }
        }
944 945
    the_end:
      ;
946
    }
947
  return (blocks_linear);
948 949
}

sebastien's avatar
sebastien committed
950
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
951 952 953
                     NumericalConstants &num_constants_arg,
                     ExternalFunctionsTable &external_functions_table_arg) :
  DataTree(symbol_table_arg, num_constants_arg, external_functions_table_arg)
sebastien's avatar
sebastien committed
954
{
955
  for (int i = 0; i < 3; i++)
956
    NNZDerivatives[i] = 0;
sebastien's avatar
sebastien committed
957 958 959 960
}

int
ModelTree::equation_number() const
961
{
962
  return (equations.size());
963
}
sebastien's avatar
sebastien committed
964 965

void
966 967
ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
                           ExprNodeOutputType output_type,
968
                           const temporary_terms_t &temporary_terms) const
969
{
970
  first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
971 972 973 974 975
  if (it != first_derivatives.end())
    (it->second)->writeOutput(output, output_type, temporary_terms);
  else
    output << 0;
}
976

sebastien's avatar
sebastien committed
977
void
978
ModelTree::computeJacobian(const set<int> &vars)
sebastien's avatar
sebastien committed
979
{
980 981
  for (set<int>::const_iterator it = vars.begin();
       it != vars.end(); it++)
982
    for (int eq = 0; eq < (int) equations.size(); eq++)
sebastien's avatar
sebastien committed
983
      {
984
        expr_t d1 = equations[eq]->getDerivative(*it);
sebastien's avatar
sebastien committed
985 986
        if (d1 == Zero)
          continue;
987
        first_derivatives[make_pair(eq, *it)] = d1;
988
        ++NNZDerivatives[0];
sebastien's avatar
sebastien committed
989
      }
990
}
sebastien's avatar
sebastien committed
991

992 993 994
void
ModelTree::computeHessian(const set<int> &vars)
{
995
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
996
       it != first_derivatives.end(); it++)
sebastien's avatar
sebastien committed
997
    {
998 999
      int eq = it->first.first;
      int var1 = it->first.second;
1000
      expr_t d1 = it->second;
1001 1002

      // Store only second derivatives with var2 <= var1
1003 1004
      for (set<int>::const_iterator it2 = vars.begin();
           it2 != vars.end(); it2++)
sebastien's avatar
sebastien committed
1005
        {
1006 1007 1008 1009
          int var2 = *it2;
          if (var2 > var1)
            continue;

1010
          expr_t d2 = d1->getDerivative(var2);
1011 1012 1013
          if (d2 == Zero)
            continue;
          second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
1014 1015 1016 1017
          if (var2 == var1)
            ++NNZDerivatives[1];
          else
            NNZDerivatives[1] += 2;
sebastien's avatar
sebastien committed
1018 1019
        }
    }
1020
}
sebastien's avatar
sebastien committed
1021

1022 1023 1024
void
ModelTree::computeThirdDerivatives(const set<int> &vars)
{