StaticModel.cc 54.1 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
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
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
43
44
  global_temporary_terms(true),
  cutoff(1e-15),
  mfs(0)
45
46
{
}
47

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

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

74
75
76
77
78
79
80
81
82
83
void
StaticModel::initializeVariablesAndEquations()
{
  for(int j = 0; j < equation_number(); j++)
    {
      equation_reordered.push_back(j);
      variable_reordered.push_back(j);
    }
}

84
85
86
void
StaticModel::computeTemporaryTermsOrdered()
{
87
88
  map<expr_t, pair<int, int> > first_occurence;
  map<expr_t, int> reference_count;
89
  BinaryOpNode *eq_node;
90
91
  first_derivatives_t::const_iterator it;
  first_chain_rule_derivatives_t::const_iterator it_chr;
92
93
94
95
96
  ostringstream tmp_s;
  v_temporary_terms.clear();
  map_idx.clear();

  unsigned int nb_blocks = getNbBlocks();
97
  v_temporary_terms = vector< vector<temporary_terms_t> >(nb_blocks);
98

99
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
100
101

  temporary_terms.clear();
102
  if (!global_temporary_terms)
103
104
105
106
107
108
109
110
111
    {
      for (unsigned int block = 0; block < nb_blocks; block++)
        {

          reference_count.clear();
          temporary_terms.clear();
          unsigned int block_size = getBlockSize(block);
          unsigned int block_nb_mfs = getBlockMfs(block);
          unsigned int block_nb_recursives = block_size - block_nb_mfs;
112
          v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
113
114
          for (unsigned int i = 0; i < block_size; i++)
            {
115
              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
116
                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
117
118
              else
                {
119
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
120
121
122
                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
                }
            }
123
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
124
            {
125
              expr_t id = it->second.second;
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
            }
          set<int> temporary_terms_in_use;
          temporary_terms_in_use.clear();
          v_temporary_terms_inuse[block] = temporary_terms_in_use;
        }
    }
  else
    {
      for (unsigned int block = 0; block < nb_blocks; block++)
        {
          // 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;
141
          v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
142
143
          for (unsigned int i = 0; i < block_size; i++)
            {
144
              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
145
                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
146
147
              else
                {
148
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
149
150
151
                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
                }
            }
152
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
153
            {
154
              expr_t id = it->second.second;
155
156
157
158
159
160
161
162
163
164
165
166
167
              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++)
            {
168
              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
169
                getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
170
171
              else
                {
172
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
173
174
175
                  eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
                }
            }
176
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
177
            {
178
              expr_t id = it->second.second;
179
180
              id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
            }
181
          for (int i = 0; i < (int) getBlockSize(block); i++)
182
            for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
183
184
                 it != v_temporary_terms[block][i].end(); it++)
              (*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
185
186
          v_temporary_terms_inuse[block] = temporary_terms_in_use;
        }
187
      computeTemporaryTermsMapping();
188
    }
189
190
191
192
193
}

void
StaticModel::computeTemporaryTermsMapping()
{
194
  // Add a mapping form node ID to temporary terms order
195
  int j = 0;
196
  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
197
      it != temporary_terms.end(); it++)
198
    map_idx[(*it)->idx] = j++;
199
200
201
202
}

void
StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
203
204
205
{
  string tmp_s, sps;
  ostringstream tmp_output, tmp1_output, global_output;
206
  expr_t lhs = NULL, rhs = NULL;
207
  BinaryOpNode *eq_node;
208
  map<expr_t, int> reference_count;
209
  temporary_terms_t local_temporary_terms;
210
211
212
213
  ofstream  output;
  int nze;
  vector<int> feedback_variables;
  ExprNodeOutputType local_output_type;
214

Sébastien Villemot's avatar
Sébastien Villemot committed
215
  local_output_type = oMatlabStaticModelSparse;
216
  if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
217
    local_temporary_terms = temporary_terms;
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
262
263
264
265
266
267
268
269
270
271
272

  //----------------------------------------------------------------------
  //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
      nze = derivative_endo[block].size();
      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)
        output << "  g1 = zeros(" << block_mfs << ", " << block_mfs << ");" << endl;

      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
273
          for (temporary_terms_inuse_t::const_iterator it = v_temporary_terms_inuse[block].begin();
274
275
276
277
278
279
280
281
282
283
284
285
286
               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];
287
          temporary_terms_t tt2;
288
289
290
291
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
292
              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
                   it != v_temporary_terms[block][i].end(); it++)
                {
                  output << "  " <<  sps;
                  (*it)->writeOutput(output, local_output_type, local_temporary_terms);
                  output << " = ";
                  (*it)->writeOutput(output, local_output_type, tt2);
                  // 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));
309
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
          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("");
337
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
                      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:
381
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
382
383
384
385
386
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
387
              expr_t id = it->second.second;
388
389
390
391
392
393
394
395
396
397
398
399
400
401
              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;
        }
      output.close();
    }
}
402
403

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

  ostringstream tmp_output;
  ofstream code_file;
409
  unsigned int instruction_number = 0;
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
  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
  FDIMT_ fdimt(temporary_terms.size());
428
  fdimt.write(code_file, instruction_number);
429
430
431
432
433
434
435
436
437
438
  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,
439
440
                           u_count_int,
                           symbol_table.endo_nbr()
441
                           );
442
  fbeginblock.write(code_file, instruction_number);
443
444
445

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

451
  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
452
453

  FENDEQU_ fendequ;
454
  fendequ.write(code_file, instruction_number);
455
456
457
458

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

475
          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
476
477

          FSTPSU_ fstpsu(count_u);
478
          fstpsu.write(code_file, instruction_number);
479
480
481
482
483
484
          count_u++;
        }
    }
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
    {
      FLDR_ fldr(i);
485
      fldr.write(code_file, instruction_number);
486
487
488
489
      for(vector<pair<int, int> >::const_iterator it = derivatives[i].begin();
          it != derivatives[i].end(); it++)
        {
          FLDSU_ fldsu(it->second);
490
          fldsu.write(code_file, instruction_number);
491
          FLDSV_ fldsv(eEndogenous, it->first);
492
          fldsv.write(code_file, instruction_number);
493
          FBINARY_ fbinary(oTimes);
494
          fbinary.write(code_file, instruction_number);
495
496
497
          if (it != derivatives[i].begin())
            {
              FBINARY_ fbinary(oPlus);
498
              fbinary.write(code_file, instruction_number);
499
500
501
            }
        }
      FBINARY_ fbinary(oMinus);
502
      fbinary.write(code_file, instruction_number);
503
      FSTPSU_ fstpsu(i);
504
      fstpsu.write(code_file, instruction_number);
505
506
    }
  FENDBLOCK_ fendblock;
507
  fendblock.write(code_file, instruction_number);
508
  FEND_ fend;
509
  fend.write(code_file, instruction_number);
510
511
512
513
  code_file.close();
}

void
514
StaticModel::writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx) const
515
516
{
  struct Uff_l
517
  {
518
519
520
    int u, var, lag;
    Uff_l *pNext;
  };
521

522
523
524
525
526
527
528
529
530
  struct Uff
  {
    Uff_l *Ufl, *Ufl_First;
  };

  int i, v;
  string tmp_s;
  ostringstream tmp_output;
  ofstream code_file;
531
  unsigned int instruction_number = 0;
532
  expr_t lhs = NULL, rhs = NULL;
533
534
  BinaryOpNode *eq_node;
  Uff Uf[symbol_table.endo_nbr()];
535
  map<expr_t, int> reference_count;
536
537
538
539
540
541
542
543
544
545
546
547
548
549
  vector<int> feedback_variables;
  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

  FDIMT_ fdimt(temporary_terms.size());
550
  fdimt.write(code_file, instruction_number);
551
552
553
554
555
556
557

  for (unsigned int block = 0; block < getNbBlocks(); block++)
    {
      feedback_variables.clear();
      if (block > 0)
        {
          FENDBLOCK_ fendblock;
558
          fendblock.write(code_file, instruction_number);
559
560
561
562
563
564
565
566
567
568
569
        }
      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)
        {
570
          Write_Inf_To_Bin_File_Block(file_name, bin_basename, block, u_count_int, file_open);
571
572
573
574
575
576
577
578
579
580
581
582
583
          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,
584
585
                               u_count_int,
                               symbol_table.endo_nbr()
586
                               );
587
      fbeginblock.write(code_file, instruction_number);
588
589
590
591
592

      // The equations
      for (i = 0; i < (int) block_size; i++)
        {
          //The Temporary terms
593
          temporary_terms_t tt2;
594
595
596
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
597
              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
598
599
                   it != v_temporary_terms[block][i].end(); it++)
                {
600
                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
601
602
                  fnumexpr.write(code_file, instruction_number);
                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false);
603
                  FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
604
                  fstpst.write(code_file, instruction_number);
605
606
607
608
609
610
611
612
613
614
615
616
617
                  // Insert current node into tt2
                  tt2.insert(*it);
                }
            }

          int variable_ID, equation_ID;
          EquationType equ_type;
          switch (simulation_type)
            {
            evaluation:
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
              equ_type = getBlockEquationType(block, i);
618
619
              {
                FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
620
                fnumexpr.write(code_file, instruction_number);
621
              }
622
623
              if (equ_type == E_EVALUATE)
                {
624
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
625
626
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
627
628
                  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);
629
630
631
                }
              else if (equ_type == E_EVALUATE_S)
                {
632
                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
633
634
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
635
636
                  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);
637
638
639
640
641
642
643
644
645
646
647
648
649
                }
              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:
650
              FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
651
              fnumexpr.write(code_file, instruction_number);
652
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
653
654
              lhs = eq_node->get_arg1();
              rhs = eq_node->get_arg2();
655
656
              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);
657
658

              FBINARY_ fbinary(oMinus);
659
              fbinary.write(code_file, instruction_number);
660
661

              FSTPR_ fstpr(i - block_recursive);
662
              fstpr.write(code_file, instruction_number);
663
664
665
            }
        }
      FENDEQU_ fendequ;
666
      fendequ.write(code_file, instruction_number);
667
668
669
670
671
672
673
674
      // 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:
675
676
              {
                FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
677
                fnumexpr.write(code_file, instruction_number);
678
              }
679
              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
680
              {
681
                FSTPG_ fstpg(0);
682
                fstpg.write(code_file, instruction_number);
683
              }
684
              break;
685

686
687
688
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              count_u = feedback_variables.size();
689
              for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
690
691
692
693
694
                {
                  unsigned int eq = it->first.first;
                  unsigned int var = it->first.second;
                  unsigned int eqr = getBlockEquationID(block, eq);
                  unsigned int varr = getBlockVariableID(block, var);
695
                  if (eq >= block_recursive && var >= block_recursive)
696
697
698
699
700
701
702
703
704
705
706
707
708
709
                    {
                      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;
710
                      FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
711
712
                      fnumexpr.write(code_file, instruction_number);
                      compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx);
713
                      FSTPSU_ fstpsu(count_u);
714
                      fstpsu.write(code_file, instruction_number);
715
716
717
718
719
720
721
722
                      count_u++;
                    }
                }
              for (i = 0; i < (int) block_size; i++)
                {
                  if (i >= (int) block_recursive)
                    {
                      FLDR_ fldr(i-block_recursive);
723
                      fldr.write(code_file, instruction_number);
724
725

                      FLDZ_ fldz;
726
                      fldz.write(code_file, instruction_number);
727
728
729
730
731

                      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);
732
                          fldsu.write(code_file, instruction_number);
733
                          FLDSV_ fldsv(eEndogenous, Uf[v].Ufl->var);
734
                          fldsv.write(code_file, instruction_number);
735
736

                          FBINARY_ fbinary(oTimes);
737
                          fbinary.write(code_file, instruction_number);
738
739

                          FCUML_ fcuml;
740
                          fcuml.write(code_file, instruction_number);
741
742
743
744
745
746
747
748
749
                        }
                      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);
750
                      fbinary.write(code_file, instruction_number);
751
752

                      FSTPSU_ fstpsu(i - block_recursive);
753
                      fstpsu.write(code_file, instruction_number);
754
755
756
757
758
759
760
761
762
763

                    }
                }
              break;
            default:
              break;
            }
        }
    }
  FENDBLOCK_ fendblock;
764
  fendblock.write(code_file, instruction_number);
765
  FEND_ fend;
766
  fend.write(code_file, instruction_number);
767
768
  code_file.close();
}
769
770

void
771
StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num,
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
                                   int &u_count_int, bool &file_open) const
{
  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;
789
  for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[num].begin(); it != (blocks_derivatives[num]).end(); it++)
790
791
792
793
    {
      unsigned int eq = it->first.first;
      unsigned int var = it->first.second;
      int lag = 0;
794
      if (eq >= block_recursive && var >= block_recursive)
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
        {
          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();
}
819

820
map<pair<int, pair<int, int > >, expr_t>
821
StaticModel::collect_first_order_derivatives_endogenous()
sebastien's avatar
sebastien committed
822
{
823
  map<pair<int, pair<int, int > >, expr_t> endo_derivatives;
824
  for (first_derivatives_t::iterator it2 = first_derivatives.begin();
825
826
       it2 != first_derivatives.end(); it2++)
    {
827
      if (getTypeByDerivID(it2->first.second) == eEndogenous)
828
829
        {
          int eq = it2->first.first;
830
          int var = symbol_table.getTypeSpecificID(it2->first.second);
831
832
833
834
          int lag = 0;
          endo_derivatives[make_pair(eq, make_pair(var, lag))] = it2->second;
        }
    }
835
  return endo_derivatives;
836
837
838
}

void
839
StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode)
840
{
841
  // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
842
843
  set<int> vars;

844
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
845
846
847
848
    vars.insert(symbol_table.getID(eEndogenous, i));

  // Launch computations
  cout << "Computing static model derivatives:" << endl
849
       << " - order 1" << endl;
850
851
852
853
854
855
856
857
858
859
860
861
  first_derivatives.clear();

  computeJacobian(vars);

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

  if (block)
    {
862
      jacob_map_t contemporaneous_jacobian, static_jacobian;
863
      vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
864
865
866
867
868
869

      // 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
870
      computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
871

872
      computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
873

874
      map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
875

876
      equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
877
878
879

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

880
881
882
      lag_lead_vector_t equation_lag_lead, variable_lag_lead;

      computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
883

884
      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type);
885
886
887
888
889
890
891
892
893
894
895
896
897
898

      printBlockDecomposition(blocks);

      computeChainRuleJacobian(blocks_derivatives);

      blocks_linear = BlockLinear(blocks_derivatives, variable_reordered);

      collect_block_first_order_derivatives();

      global_temporary_terms = true;
      if (!no_tmp_terms)
        computeTemporaryTermsOrdered();
    }
  else
899
900
901
902
903
904
905
906
    {
      if (!no_tmp_terms)
        {
          computeTemporaryTerms(true);
          if (bytecode)
            computeTemporaryTermsMapping();
        }
    }
sebastien's avatar
sebastien committed
907
908
909
}

void
910
StaticModel::writeStaticMFile(const string &func_name) const
sebastien's avatar
sebastien committed
911
912
{
  // Writing comments and function definition command
913
914
915
916
917
918
919
920
921
922
923
  string filename = func_name + "_static.m";

  ofstream output;
  output.open(filename.c_str(), ios::out | ios::binary);
  if (!output.is_open())
    {
      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
      exit(EXIT_FAILURE);
    }

  output << "function [residual, g1, g2] = " << func_name + "_static(y, x, params)" << endl
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
         << "%" << endl
         << "% Status : Computes static model for Dynare" << endl
         << "%" << endl
         << "% Warning : this file is generated automatically by Dynare" << endl
         << "%           from model file (.mod)" << endl
         << endl
         << "residual = zeros( " << equations.size() << ", 1);" << endl
         << endl
         << "%" << endl
         << "% Model equations" << endl
         << "%" << endl
         << endl;

  writeModelLocalVariables(output, oMatlabStaticModel);

  writeTemporaryTerms(temporary_terms, output, oMatlabStaticModel);

  writeModelEquations(output, oMatlabStaticModel);

  output << "if ~isreal(residual)" << endl
         << "  residual = real(residual)+imag(residual).^2;" << endl
         << "end" << endl
         << "if nargout >= 2," << endl
         << "  g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl
         << endl
         << "%" << endl
         << "% Jacobian matrix" << endl
         << "%" << endl
         << endl;
sebastien's avatar
sebastien committed
953
954

  // Write Jacobian w.r. to endogenous only
955
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
sebastien's avatar
sebastien committed
956
957
958
       it != first_derivatives.end(); it++)
    {
      int eq = it->first.first;
959
      int symb_id = it->first.second;
960
      expr_t d1 = it->second;
sebastien's avatar
sebastien committed
961

962
      output << "  g1(" << eq+1 << "," << symbol_table.getTypeSpecificID(symb_id)+1 << ")=";
963
964
      d1->writeOutput(output, oMatlabStaticModel, temporary_terms);
      output << ";" << endl;
sebastien's avatar
sebastien committed
965
966
    }

967
968
969
970
971
972
973
974
975
976
977
978
  output << "  if ~isreal(g1)" << endl
         << "    g1 = real(g1)+2*imag(g1);" << endl
         << "  end" << endl
         << "end" << endl
         << "if nargout >= 3," << endl
         << "%" << endl
         << "% Hessian matrix" << endl
         << "%" << endl
         << endl;

  int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
  if (second_derivatives.size())
979
    {
980
      output << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl;
981

982
983
      // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
      int k = 0; // Keep the line of a 2nd derivative in v2
984
      for (second_derivatives_t::const_iterator it = second_derivatives.begin();
985
           it != second_derivatives.end(); it++)
986
        {
987
988
989
          int eq = it->first.first;
          int symb_id1 = it->first.second.first;
          int symb_id2 = it->first.second.second;
990
          expr_t d2 = it->second;
991
992
993

          int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
          int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
sebastien's avatar
sebastien committed
994

995
996
          int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
          int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
sebastien's avatar
sebastien committed
997

998
999
1000
          output << "v2(" << k+1 << ",1)=" << eq + 1 << ";" << endl
                 << "v2(" << k+1 << ",2)=" << col_nb + 1 << ";" << endl
                 << "v2(" << k+1 << ",3)=";