StaticModel.cc 69.8 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2011 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
21
#include <iostream>
#include <cmath>
22
#include <cstdlib>
23
#include <cassert>
24
25
#include <cstdio>
#include <cerrno>
sebastien's avatar
sebastien committed
26
#include <algorithm>
27
28
29
30
31
32
33
34
35
#include "StaticModel.hh"

// For mkdir() and chdir()
#ifdef _WIN32
# include <direct.h>
#else
# include <unistd.h>
# include <sys/stat.h>
# include <sys/types.h>
36
#endif
sebastien's avatar
sebastien committed
37

38
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
39
40
41
                         NumericalConstants &num_constants_arg,
                         ExternalFunctionsTable &external_functions_table_arg) :
  ModelTree(symbol_table_arg, num_constants_arg, external_functions_table_arg),
42
  global_temporary_terms(true)
43
44
{
}
45

46
void
47
StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
48
{
49
  first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id)));
50
  if (it != first_derivatives.end())
51
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
52
53
54
  else
    {
      FLDZ_ fldz;
55
      fldz.write(code_file, instruction_number);
56
57
    }
}
sebastien's avatar
sebastien committed
58

59
void
60
StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
61
{
62
  map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
63
  if (it != first_chain_rule_derivatives.end())
64
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
65
66
67
  else
    {
      FLDZ_ fldz;
68
      fldz.write(code_file, instruction_number);
69
70
71
72
73
74
    }
}

void
StaticModel::computeTemporaryTermsOrdered()
{
75
76
  map<expr_t, pair<int, int> > first_occurence;
  map<expr_t, int> reference_count;
77
  BinaryOpNode *eq_node;
78
79
  first_derivatives_t::const_iterator it;
  first_chain_rule_derivatives_t::const_iterator it_chr;
80
81
82
83
84
  ostringstream tmp_s;
  v_temporary_terms.clear();
  map_idx.clear();

  unsigned int nb_blocks = getNbBlocks();
85
  v_temporary_terms = vector< vector<temporary_terms_t> >(nb_blocks);
86
  v_temporary_terms_local = vector< vector<temporary_terms_t> >(nb_blocks);
87

88
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
89

90
91
  map_idx2 = vector<map_idx_t>(nb_blocks);

92
  temporary_terms.clear();
93
94
95

  //local temporay terms
  for (unsigned int block = 0; block < nb_blocks; block++)
96
    {
97
98
99
100
101
102
      map<expr_t, int> reference_count_local;
      reference_count_local.clear();
      map<expr_t, pair<int, int> > first_occurence_local;
      first_occurence_local.clear();
      temporary_terms_t temporary_terms_l;
      temporary_terms_l.clear();
103

104
105
106
107
108
109
110
111
112
113
      unsigned int block_size = getBlockSize(block);
      unsigned int block_nb_mfs = getBlockMfs(block);
      unsigned int block_nb_recursives = block_size - block_nb_mfs;
      v_temporary_terms_local[block] = vector<temporary_terms_t>(block_size);

      for (unsigned int i = 0; i < block_size; i++)
        {
          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  i);
          else
114
            {
115
116
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  i);
117
118
            }
        }
119
120
121
122
123
124
125
126
127
      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
        {
          expr_t id = it->second.second;
          id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  block_size-1);
        }
      set<int> temporary_terms_in_use;
      temporary_terms_in_use.clear();
      v_temporary_terms_inuse[block] = temporary_terms_in_use;
      computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
128
    }
129
130
131

  // global temporay terms
  for (unsigned int block = 0; block < nb_blocks; block++)
132
    {
133
134
135
136
137
138
      // Compute the temporary terms reordered
      unsigned int block_size = getBlockSize(block);
      unsigned int block_nb_mfs = getBlockMfs(block);
      unsigned int block_nb_recursives = block_size - block_nb_mfs;
      v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
      for (unsigned int i = 0; i < block_size; i++)
139
        {
140
141
142
          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
          else
143
            {
144
145
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
146
147
            }
        }
148
      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
149
        {
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
          expr_t id = it->second.second;
          id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
        }
    }

  for (unsigned int block = 0; block < nb_blocks; block++)
    {
      // Collecte the temporary terms reordered
      unsigned int block_size = getBlockSize(block);
      unsigned int block_nb_mfs = getBlockMfs(block);
      unsigned int block_nb_recursives = block_size - block_nb_mfs;
      set<int> temporary_terms_in_use;
      for (unsigned int i = 0; i < block_size; i++)
        {
          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
            getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          else
167
            {
168
169
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
170
171
            }
        }
172
173
174
175
176
177
178
179
180
181
      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
        {
          expr_t id = it->second.second;
          id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
        }
      for (int i = 0; i < (int) getBlockSize(block); i++)
        for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
             it != v_temporary_terms[block][i].end(); it++)
          (*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
      v_temporary_terms_inuse[block] = temporary_terms_in_use;
182
    }
183
  computeTemporaryTermsMapping(temporary_terms, map_idx);
184
185
186
}

void
187
StaticModel::computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx)
188
{
189
  // Add a mapping form node ID to temporary terms order
190
  int j = 0;
191
  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
192
       it != temporary_terms.end(); it++)
193
    map_idx[(*it)->idx] = j++;
194
195
196
197
}

void
StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
198
199
200
{
  string tmp_s, sps;
  ostringstream tmp_output, tmp1_output, global_output;
201
  expr_t lhs = NULL, rhs = NULL;
202
  BinaryOpNode *eq_node;
203
  map<expr_t, int> reference_count;
204
  temporary_terms_t local_temporary_terms;
205
206
  ofstream  output;
  vector<int> feedback_variables;
207
  deriv_node_temp_terms_t tef_terms;
208
  ExprNodeOutputType local_output_type;
209

Sébastien Villemot's avatar
Sébastien Villemot committed
210
  local_output_type = oMatlabStaticModelSparse;
211
  if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
212
    local_temporary_terms = temporary_terms;
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261

  //----------------------------------------------------------------------
  //For each block
  for (unsigned int block = 0; block < getNbBlocks(); block++)
    {
      //recursive_variables.clear();
      feedback_variables.clear();
      //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
      BlockSimulationType simulation_type = getBlockSimulationType(block);
      unsigned int block_size = getBlockSize(block);
      unsigned int block_mfs = getBlockMfs(block);
      unsigned int block_recursive = block_size - block_mfs;

      tmp1_output.str("");
      tmp1_output << static_basename << "_" << block+1 << ".m";
      output.open(tmp1_output.str().c_str(), ios::out | ios::binary);
      output << "%\n";
      output << "% " << tmp1_output.str() << " : Computes static model for Dynare\n";
      output << "%\n";
      output << "% Warning : this file is generated automatically by Dynare\n";
      output << "%           from model file (.mod)\n\n";
      output << "%/\n";
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        output << "function y = " << static_basename << "_" << block+1 << "(y, x, params)\n";
      else
        output << "function [residual, y, g1] = " << static_basename << "_" << block+1 << "(y, x, params)\n";

      BlockType block_type;
      if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
        block_type = SIMULTANS;
      else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
                || simulation_type == EVALUATE_BACKWARD    || simulation_type == EVALUATE_FORWARD)
               && getBlockFirstEquation(block) < prologue)
        block_type = PROLOGUE;
      else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
                || simulation_type == EVALUATE_BACKWARD    || simulation_type == EVALUATE_FORWARD)
               && getBlockFirstEquation(block) >= equations.size() - epilogue)
        block_type = EPILOGUE;
      else
        block_type = SIMULTANS;
      output << "  % ////////////////////////////////////////////////////////////////////////" << endl
             << "  % //" << string("                     Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
             << "          //" << endl
             << "  % //                     Simulation type "
             << BlockSim(simulation_type) << "  //" << endl
             << "  % ////////////////////////////////////////////////////////////////////////" << endl;
      output << "  global options_;" << endl;
      //The Temporary terms
      if (simulation_type != EVALUATE_BACKWARD  && simulation_type != EVALUATE_FORWARD)
262
        output << " g1 = spalloc("  << block_mfs << ", " << block_mfs << ", " << derivative_endo[block].size() << ");" << endl;
263
264
265
266

      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
267
          for (temporary_terms_inuse_t::const_iterator it = v_temporary_terms_inuse[block].begin();
268
269
270
271
272
273
274
275
276
277
278
279
280
               it != v_temporary_terms_inuse[block].end(); it++)
            tmp_output << " T" << *it;
          output << "  global" << tmp_output.str() << ";\n";
        }

      if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
        output << "  residual=zeros(" << block_mfs << ",1);\n";

      // The equations
      for (unsigned int i = 0; i < block_size; i++)
        {
          if (!global_temporary_terms)
            local_temporary_terms = v_temporary_terms[block][i];
281
          temporary_terms_t tt2;
282
283
284
285
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
286
              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
287
288
                   it != v_temporary_terms[block][i].end(); it++)
                {
289
290
291
                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);

292
                  output << "  " <<  sps;
293
                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
294
                  output << " = ";
295
                  (*it)->writeOutput(output, local_output_type, tt2, tef_terms);
296
297
298
299
300
301
302
303
304
305
                  // Insert current node into tt2
                  tt2.insert(*it);
                  output << ";" << endl;
                }
            }

          int variable_ID = getBlockVariableID(block, i);
          int equation_ID = getBlockEquationID(block, i);
          EquationType equ_type = getBlockEquationType(block, i);
          string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, variable_ID));
306
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
          lhs = eq_node->get_arg1();
          rhs = eq_node->get_arg2();
          tmp_output.str("");
          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms);
          switch (simulation_type)
            {
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
            evaluation:
              output << "  % equation " << getBlockEquationID(block, i)+1 << " variable : " << sModel
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
              output << "  ";
              if (equ_type == E_EVALUATE)
                {
                  output << tmp_output.str();
                  output << " = ";
                  rhs->writeOutput(output, local_output_type, local_temporary_terms);
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  output << "%" << tmp_output.str();
                  output << " = ";
                  if (isBlockEquationRenormalized(block, i))
                    {
                      rhs->writeOutput(output, local_output_type, local_temporary_terms);
                      output << "\n  ";
                      tmp_output.str("");
334
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
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
                      lhs = eq_node->get_arg1();
                      rhs = eq_node->get_arg2();
                      lhs->writeOutput(output, local_output_type, local_temporary_terms);
                      output << " = ";
                      rhs->writeOutput(output, local_output_type, local_temporary_terms);
                    }
                }
              else
                {
                  cerr << "Type missmatch for equation " << equation_ID+1  << "\n";
                  exit(EXIT_FAILURE);
                }
              output << ";\n";
              break;
            case SOLVE_BACKWARD_SIMPLE:
            case SOLVE_FORWARD_SIMPLE:
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              if (i < block_recursive)
                goto evaluation;
              feedback_variables.push_back(variable_ID);
              output << "  % equation " << equation_ID+1 << " variable : " << sModel
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
              output << "  " << "residual(" << i+1-block_recursive << ") = (";
              goto end;
            default:
            end:
              output << tmp_output.str();
              output << ") - (";
              rhs->writeOutput(output, local_output_type, local_temporary_terms);
              output << ");\n";
            }
        }
      // The Jacobian if we have to solve the block
      if (simulation_type == SOLVE_BACKWARD_SIMPLE   || simulation_type == SOLVE_FORWARD_SIMPLE
          || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
        output << "  " << sps << "% Jacobian  " << endl;
      switch (simulation_type)
        {
        case SOLVE_BACKWARD_SIMPLE:
        case SOLVE_FORWARD_SIMPLE:
        case SOLVE_BACKWARD_COMPLETE:
        case SOLVE_FORWARD_COMPLETE:
378
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
379
380
381
382
383
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
384
              expr_t id = it->second.second;
385
386
387
388
389
390
391
392
393
394
395
              output << "    g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
              id->writeOutput(output, local_output_type, local_temporary_terms);
              output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
                     << "(" << 0
                     << ") " << varr+1
                     << ", equation=" << eqr+1 << endl;
            }
          break;
        default:
          break;
        }
396
      output << "end" << endl;
397
398
399
      output.close();
    }
}
400
401

void
402
StaticModel::writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_t map_idx) const
403
404
405
406
{

  ostringstream tmp_output;
  ofstream code_file;
407
  unsigned int instruction_number = 0;
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
  bool file_open = false;

  string main_name = file_name;
  main_name += ".cod";
  code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate);
  if (!code_file.is_open())
    {
      cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
      exit(EXIT_FAILURE);
    }
  int count_u;
  int u_count_int = 0;

  Write_Inf_To_Bin_File(file_name, u_count_int, file_open, false, symbol_table.endo_nbr());
  file_open = true;

  //Temporary variables declaration
425
426
  FDIMST_ fdimst(temporary_terms.size());
  fdimst.write(code_file, instruction_number);
427
428
429
430
431
432
433
434
435
436
  FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
                           SOLVE_FORWARD_COMPLETE,
                           0,
                           symbol_table.endo_nbr(),
                           variable_reordered,
                           equation_reordered,
                           false,
                           symbol_table.endo_nbr(),
                           0,
                           0,
437
438
                           u_count_int,
                           symbol_table.endo_nbr()
439
                           );
440
  fbeginblock.write(code_file, instruction_number);
441
442
443

  // Add a mapping form node ID to temporary terms order
  int j = 0;
444
  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
445
446
       it != temporary_terms.end(); it++)
    map_idx[(*it)->idx] = j++;
447
  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
448

449
  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
450
451

  FENDEQU_ fendequ;
452
  fendequ.write(code_file, instruction_number);
453

454
455
456
457
458
459
  // Get the current code_file position and jump if eval = true
  streampos pos1 = code_file.tellp();
  FJMPIFEVAL_ fjmp_if_eval(0);
  fjmp_if_eval.write(code_file, instruction_number);
  int prev_instruction_number = instruction_number;

460
461
462
  vector<vector<pair<int, int> > > derivatives;
  derivatives.resize(symbol_table.endo_nbr());
  count_u = symbol_table.endo_nbr();
463
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
464
465
466
467
468
       it != first_derivatives.end(); it++)
    {
      int deriv_id = it->first.second;
      if (getTypeByDerivID(deriv_id) == eEndogenous)
        {
469
          expr_t d1 = it->second;
470
471
472
          unsigned int eq = it->first.first;
          int symb = getSymbIDByDerivID(deriv_id);
          unsigned int var = symbol_table.getTypeSpecificID(symb);
473
          FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
474
          fnumexpr.write(code_file, instruction_number);
475
476
477
478
          if (!derivatives[eq].size())
            derivatives[eq].clear();
          derivatives[eq].push_back(make_pair(var, count_u));

479
          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
480
481

          FSTPSU_ fstpsu(count_u);
482
          fstpsu.write(code_file, instruction_number);
483
484
485
486
487
488
          count_u++;
        }
    }
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
    {
      FLDR_ fldr(i);
489
      fldr.write(code_file, instruction_number);
490
      if (derivatives[i].size())
491
        {
492
493
          for (vector<pair<int, int> >::const_iterator it = derivatives[i].begin();
               it != derivatives[i].end(); it++)
494
            {
495
496
497
498
499
              FLDSU_ fldsu(it->second);
              fldsu.write(code_file, instruction_number);
              FLDSV_ fldsv(eEndogenous, it->first);
              fldsv.write(code_file, instruction_number);
              FBINARY_ fbinary(oTimes);
500
              fbinary.write(code_file, instruction_number);
501
502
503
504
505
              if (it != derivatives[i].begin())
                {
                  FBINARY_ fbinary(oPlus);
                  fbinary.write(code_file, instruction_number);
                }
506
            }
507
508
          FBINARY_ fbinary(oMinus);
          fbinary.write(code_file, instruction_number);
509
510
        }
      FSTPSU_ fstpsu(i);
511
      fstpsu.write(code_file, instruction_number);
512
    }
513
514
515
516
517
518
519
520
521
522
  // Get the current code_file position and jump = true
  streampos pos2 = code_file.tellp();
  FJMP_ fjmp(0);
  fjmp.write(code_file, instruction_number);
  // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
  streampos pos3 = code_file.tellp();
  code_file.seekp(pos1);
  FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
  fjmp_if_eval1.write(code_file, instruction_number);
  code_file.seekp(pos3);
523
  prev_instruction_number = instruction_number;
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547

  temporary_terms_t tt2;
  tt2.clear();
  temporary_terms_t tt3;
  tt3.clear();

  // The Jacobian if we have to solve the block determinsitic bloc
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
    {
      int deriv_id = it->first.second;
      if (getTypeByDerivID(deriv_id) == eEndogenous)
        {
          expr_t d1 = it->second;
          unsigned int eq = it->first.first;
          int symb = getSymbIDByDerivID(deriv_id);
          unsigned int var = symbol_table.getTypeSpecificID(symb);
          FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
          fnumexpr.write(code_file, instruction_number);
          if (!derivatives[eq].size())
            derivatives[eq].clear();
          derivatives[eq].push_back(make_pair(var, count_u));

          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
548
          FSTPG2_ fstpg2(eq, var);
549
550
551
552
553
554
555
556
557
558
559
          fstpg2.write(code_file, instruction_number);
        }
    }

  // Set codefile position to previous JMP_ and set the number of instructions to jump
  pos1 = code_file.tellp();
  code_file.seekp(pos2);
  FJMP_ fjmp1(instruction_number - prev_instruction_number);
  fjmp1.write(code_file, instruction_number);
  code_file.seekp(pos1);

560
  FENDBLOCK_ fendblock;
561
  fendblock.write(code_file, instruction_number);
562
  FEND_ fend;
563
  fend.write(code_file, instruction_number);
564
565
566
567
  code_file.close();
}

void
568
StaticModel::writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const
569
570
{
  struct Uff_l
571
  {
572
573
574
    int u, var, lag;
    Uff_l *pNext;
  };
575

576
577
578
579
580
581
582
583
584
  struct Uff
  {
    Uff_l *Ufl, *Ufl_First;
  };

  int i, v;
  string tmp_s;
  ostringstream tmp_output;
  ofstream code_file;
585
  unsigned int instruction_number = 0;
586
  expr_t lhs = NULL, rhs = NULL;
587
588
  BinaryOpNode *eq_node;
  Uff Uf[symbol_table.endo_nbr()];
589
  map<expr_t, int> reference_count;
590
  vector<int> feedback_variables;
591
  deriv_node_temp_terms_t tef_terms;
592
593
594
595
596
597
598
599
600
601
602
603
  bool file_open = false;

  string main_name = file_name;
  main_name += ".cod";
  code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate);
  if (!code_file.is_open())
    {
      cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
      exit(EXIT_FAILURE);
    }
  //Temporary variables declaration

604
605
  FDIMST_ fdimst(temporary_terms.size());
  fdimst.write(code_file, instruction_number);
606
607
608
609
610
611
612

  for (unsigned int block = 0; block < getNbBlocks(); block++)
    {
      feedback_variables.clear();
      if (block > 0)
        {
          FENDBLOCK_ fendblock;
613
          fendblock.write(code_file, instruction_number);
614
615
616
617
618
619
620
621
622
623
624
        }
      int count_u;
      int u_count_int = 0;
      BlockSimulationType simulation_type = getBlockSimulationType(block);
      unsigned int block_size = getBlockSize(block);
      unsigned int block_mfs = getBlockMfs(block);
      unsigned int block_recursive = block_size - block_mfs;

      if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE
          || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
        {
625
          Write_Inf_To_Bin_File_Block(file_name, bin_basename, block, u_count_int, file_open);
626
627
628
629
630
631
632
633
634
635
636
637
638
          file_open = true;
        }

      FBEGINBLOCK_ fbeginblock(block_mfs,
                               simulation_type,
                               getBlockFirstEquation(block),
                               block_size,
                               variable_reordered,
                               equation_reordered,
                               blocks_linear[block],
                               symbol_table.endo_nbr(),
                               0,
                               0,
639
                               u_count_int,
640
                               /*symbol_table.endo_nbr()*/ block_size
641
                               );
642

643
      fbeginblock.write(code_file, instruction_number);
644

645
646
647
648
649
650
      // Get the current code_file position and jump if eval = true
      streampos pos1 = code_file.tellp();
      FJMPIFEVAL_ fjmp_if_eval(0);
      fjmp_if_eval.write(code_file, instruction_number);
      int prev_instruction_number = instruction_number;

651
652
653
      for (i = 0; i < (int) block_size; i++)
        {
          //The Temporary terms
654
          temporary_terms_t tt2;
655
656
657
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
658
              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
659
660
                   it != v_temporary_terms[block][i].end(); it++)
                {
661
662
663
                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);

664
                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int) (map_idx.find((*it)->idx)->second));
665
                  fnumexpr.write(code_file, instruction_number);
666
                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
667
                  FSTPST_ fstpst((int) (map_idx.find((*it)->idx)->second));
668
                  fstpst.write(code_file, instruction_number);
669
670
671
672
673
                  // Insert current node into tt2
                  tt2.insert(*it);
                }
            }

674
          // The equations
675
676
677
678
679
680
681
682
          int variable_ID, equation_ID;
          EquationType equ_type;
          switch (simulation_type)
            {
            evaluation:
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
              equ_type = getBlockEquationType(block, i);
683
684
              {
                FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
685
                fnumexpr.write(code_file, instruction_number);
686
              }
687
688
              if (equ_type == E_EVALUATE)
                {
689
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
690
691
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
692
693
                  rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
                  lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
694
695
696
                }
              else if (equ_type == E_EVALUATE_S)
                {
697
                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
698
699
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
700
701
                  rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
                  lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
702
703
704
705
706
707
708
709
710
711
712
713
714
                }
              break;
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              if (i < (int) block_recursive)
                goto evaluation;
              variable_ID = getBlockVariableID(block, i);
              equation_ID = getBlockEquationID(block, i);
              feedback_variables.push_back(variable_ID);
              Uf[equation_ID].Ufl = NULL;
              goto end;
            default:
            end:
715
              FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
716
              fnumexpr.write(code_file, instruction_number);
717
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
718
719
              lhs = eq_node->get_arg1();
              rhs = eq_node->get_arg2();
720
721
              lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
              rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
722
723

              FBINARY_ fbinary(oMinus);
724
              fbinary.write(code_file, instruction_number);
725
726

              FSTPR_ fstpr(i - block_recursive);
727
              fstpr.write(code_file, instruction_number);
728
729
730
            }
        }
      FENDEQU_ fendequ;
731
      fendequ.write(code_file, instruction_number);
732

733
734
735
736
737
738
739
740
      // The Jacobian if we have to solve the block
      if    (simulation_type != EVALUATE_BACKWARD
             && simulation_type != EVALUATE_FORWARD)
        {
          switch (simulation_type)
            {
            case SOLVE_BACKWARD_SIMPLE:
            case SOLVE_FORWARD_SIMPLE:
741
742
              {
                FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
743
                fnumexpr.write(code_file, instruction_number);
744
              }
745
              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx, temporary_terms);
746
              {
747
                FSTPG_ fstpg(0);
748
                fstpg.write(code_file, instruction_number);
749
              }
750
              break;
751

752
753
754
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              count_u = feedback_variables.size();
755
              for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
756
757
758
759
760
                {
                  unsigned int eq = it->first.first;
                  unsigned int var = it->first.second;
                  unsigned int eqr = getBlockEquationID(block, eq);
                  unsigned int varr = getBlockVariableID(block, var);
761
                  if (eq >= block_recursive && var >= block_recursive)
762
763
764
765
766
767
768
769
770
771
772
773
774
775
                    {
                      if (!Uf[eqr].Ufl)
                        {
                          Uf[eqr].Ufl = (Uff_l *) malloc(sizeof(Uff_l));
                          Uf[eqr].Ufl_First = Uf[eqr].Ufl;
                        }
                      else
                        {
                          Uf[eqr].Ufl->pNext = (Uff_l *) malloc(sizeof(Uff_l));
                          Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
                        }
                      Uf[eqr].Ufl->pNext = NULL;
                      Uf[eqr].Ufl->u = count_u;
                      Uf[eqr].Ufl->var = varr;
776
                      FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
777
                      fnumexpr.write(code_file, instruction_number);
778
                      compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx, temporary_terms);
779
                      FSTPSU_ fstpsu(count_u);
780
                      fstpsu.write(code_file, instruction_number);
781
782
783
784
785
786
787
788
                      count_u++;
                    }
                }
              for (i = 0; i < (int) block_size; i++)
                {
                  if (i >= (int) block_recursive)
                    {
                      FLDR_ fldr(i-block_recursive);
789
                      fldr.write(code_file, instruction_number);
790
791

                      FLDZ_ fldz;
792
                      fldz.write(code_file, instruction_number);
793
794
795
796
797

                      v = getBlockEquationID(block, i);
                      for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
                        {
                          FLDSU_ fldsu(Uf[v].Ufl->u);
798
                          fldsu.write(code_file, instruction_number);
799
                          FLDSV_ fldsv(eEndogenous, Uf[v].Ufl->var);
800
                          fldsv.write(code_file, instruction_number);
801
802

                          FBINARY_ fbinary(oTimes);
803
                          fbinary.write(code_file, instruction_number);
804
805

                          FCUML_ fcuml;
806
                          fcuml.write(code_file, instruction_number);
807
808
809
810
811
812
813
814
815
                        }
                      Uf[v].Ufl = Uf[v].Ufl_First;
                      while (Uf[v].Ufl)
                        {
                          Uf[v].Ufl_First = Uf[v].Ufl->pNext;
                          free(Uf[v].Ufl);
                          Uf[v].Ufl = Uf[v].Ufl_First;
                        }
                      FBINARY_ fbinary(oMinus);
816
                      fbinary.write(code_file, instruction_number);
817
818

                      FSTPSU_ fstpsu(i - block_recursive);
819
                      fstpsu.write(code_file, instruction_number);
820
821
822
823
824
825
826
827

                    }
                }
              break;
            default:
              break;
            }
        }
828
829
830
831
832
833
834
835
836
837
838

      // Get the current code_file position and jump = true
      streampos pos2 = code_file.tellp();
      FJMP_ fjmp(0);
      fjmp.write(code_file, instruction_number);
      // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
      streampos pos3 = code_file.tellp();
      code_file.seekp(pos1);
      FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
      fjmp_if_eval1.write(code_file, instruction_number);
      code_file.seekp(pos3);
839
      prev_instruction_number = instruction_number;
840
841
842
843
844

      temporary_terms_t tt2;
      tt2.clear();
      temporary_terms_t tt3;
      tt3.clear();
845
      deriv_node_temp_terms_t tef_terms2;
846
847
848
849
850
851
852
853

      for (i = 0; i < (int) block_size; i++)
        {
          if (v_temporary_terms_local[block].size())
            {
              for (temporary_terms_t::const_iterator it = v_temporary_terms_local[block][i].begin();
                   it != v_temporary_terms_local[block][i].end(); it++)
                {
854
                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
855
                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
856

857
                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int) (map_idx2[block].find((*it)->idx)->second));
858
                  fnumexpr.write(code_file, instruction_number);
859
860
861

                  (*it)->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);

862
                  FSTPST_ fstpst((int) (map_idx2[block].find((*it)->idx)->second));
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
                  fstpst.write(code_file, instruction_number);
                  // Insert current node into tt2
                  tt3.insert(*it);
                  tt2.insert(*it);
                }
            }

          // The equations
          int variable_ID, equation_ID;
          EquationType equ_type;
          switch (simulation_type)
            {
            evaluation_l:
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
              equ_type = getBlockEquationType(block, i);
              {
                FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
                fnumexpr.write(code_file, instruction_number);
              }
              if (equ_type == E_EVALUATE)
                {
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
888
889
                  rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
                  lhs->compile(code_file, instruction_number, true, tt2, map_idx2[block], false, false);
890
891
892
893
894
895
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
896
897
                  rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
                  lhs->compile(code_file, instruction_number, true, tt2, map_idx2[block], false, false);
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
                }
              break;
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              if (i < (int) block_recursive)
                goto evaluation_l;
              variable_ID = getBlockVariableID(block, i);
              equation_ID = getBlockEquationID(block, i);
              feedback_variables.push_back(variable_ID);
              Uf[equation_ID].Ufl = NULL;
              goto end_l;
            default:
            end_l:
              FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
              fnumexpr.write(code_file, instruction_number);
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              lhs = eq_node->get_arg1();
              rhs = eq_node->get_arg2();
916
917
              lhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
              rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937

              FBINARY_ fbinary(oMinus);
              fbinary.write(code_file, instruction_number);

              FSTPR_ fstpr(i - block_recursive);
              fstpr.write(code_file, instruction_number);
            }
        }
      FENDEQU_ fendequ_l;
      fendequ_l.write(code_file, instruction_number);

      // The Jacobian if we have to solve the block determinsitic bloc
      switch (simulation_type)
        {
        case SOLVE_BACKWARD_SIMPLE:
        case SOLVE_FORWARD_SIMPLE:
          {
            FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
            fnumexpr.write(code_file, instruction_number);
          }
938
          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx2[block], tt2 /*temporary_terms*/);
939
          {
940
            FSTPG2_ fstpg2(0, 0);
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
            fstpg2.write(code_file, instruction_number);
          }
          break;
        case EVALUATE_BACKWARD:
        case EVALUATE_FORWARD:
        case SOLVE_BACKWARD_COMPLETE:
        case SOLVE_FORWARD_COMPLETE:
          count_u = feedback_variables.size();
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
              FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
              fnumexpr.write(code_file, instruction_number);

958
              compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2 /*temporary_terms*/);
959

960
              FSTPG2_ fstpg2(eq, var);
961
962
963
964
965
966
967
968
969
970
971
972
              fstpg2.write(code_file, instruction_number);
            }
          break;
        default:
          break;
        }
      // Set codefile position to previous JMP_ and set the number of instructions to jump
      pos1 = code_file.tellp();
      code_file.seekp(pos2);
      FJMP_ fjmp1(instruction_number - prev_instruction_number);
      fjmp1.write(code_file, instruction_number);
      code_file.seekp(pos1);
973
974
    }
  FENDBLOCK_ fendblock;
975
  fendblock.write(code_file, instruction_number);
976
  FEND_ fend;
977
  fend.write(code_file, instruction_number);
978
979
  code_file.close();
}
980
981

void
982
StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num,
983
                                         int &u_count_int, bool &file_open) const
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
{
  int j;
  std::ofstream SaveCode;
  if (file_open)
    SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::in | ios::binary | ios::ate);
  else
    SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::binary);
  if (!SaveCode.is_open())
    {
      cout << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing\n";
      exit(EXIT_FAILURE);
    }
  u_count_int = 0;
  unsigned int block_size = getBlockSize(num);
  unsigned int block_mfs = getBlockMfs(num);
  unsigned int block_recursive = block_size - block_mfs;
1000
  for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[num].begin(); it != (blocks_derivatives[num]).end(); it++)
1001
1002
1003
1004
    {
      unsigned int eq = it->first.first;
      unsigned int var = it->first.second;
      int lag = 0;
1005
      if (eq >= block_recursive && var >= block_recursive)
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
        {
          int v = eq - block_recursive;
          SaveCode.write(reinterpret_cast<char *>(&v), sizeof(v));
          int varr = var - block_recursive;
          SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
          SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
          int u = u_count_int + block_mfs;
          SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
          u_count_int++;
        }
    }

  for (j = block_recursive; j < (int) block_size; j++)
    {
      unsigned int varr = getBlockVariableID(num, j);
      SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
    }
  for (j = block_recursive; j < (int) block_size; j++)
    {
      unsigned int eqr = getBlockEquationID(num, j);
      SaveCode.write(reinterpret_cast<char *>(&eqr), sizeof(eqr));
    }
  SaveCode.close();
}
1030

1031
map<pair<int, pair<int, int > >, expr_t>
1032
StaticModel::collect_first_order_derivatives_endogenous()
sebastien's avatar
sebastien committed
1033
{
1034
  map<pair<int, pair<int, int > >, expr_t> endo_derivatives;
1035
  for (first_derivatives_t::iterator it2 = first_derivatives.begin();
1036
1037
       it2 != first_derivatives.end(); it2++)
    {
1038
      if (getTypeByDerivID(it2->first.second) == eEndogenous)
1039
1040
        {
          int eq = it2->first.first;
1041
          int var = symbol_table.getTypeSpecificID(it2->first.second);
1042
1043
1044
1045
          int lag = 0;
          endo_derivatives[make_pair(eq, make_pair(var, lag))] = it2->second;
        }
    }
1046
  return endo_derivatives;
1047
1048
1049
}

void
1050
StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode)
1051
{
1052
1053
  initializeVariablesAndEquations();

1054
  // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
1055
1056
  set<int> vars;

1057
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
1058
1059
1060
1061
    vars.insert(symbol_table.getID(eEndogenous, i));

  // Launch computations
  cout << "Computing static model derivatives:" << endl
1062
       << " - order 1" << endl;
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
  first_derivatives.clear();

  computeJacobian(vars);

  if (hessian)
    {
      cout << " - order 2" << endl;
      computeHessian(vars);
    }

  if (block)
    {
1075
      jacob_map_t contemporaneous_jacobian, static_jacobian;
1076
      vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
1077
1078
1079
1080
1081
1082

      // for each block contains pair<Size, Feddback_variable>
      vector<pair<int, int> > blocks;

      evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);

sebastien's avatar
sebastien committed
1083
      computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
1084

1085
      computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
1086

1087
      map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
1088

1089
      equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
1090
1091
1092

      cout << "Finding the optimal block decomposition of the model ...\n";

Ferhat Mihoubi's avatar