StaticModel.cc 142 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
Houtan Bastani's avatar
Houtan Bastani committed
2
 * Copyright (C) 2003-2018 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
#include <boost/filesystem.hpp>

#include "StaticModel.hh"
sebastien's avatar
sebastien committed
31

32
33
34
void
StaticModel::copyHelper(const StaticModel &m)
{
35
  auto f = [this](expr_t e) { return e->clone(*this); };
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

  auto convert_vector_tt = [this,f](vector<temporary_terms_t> vtt)
    {
      vector<temporary_terms_t> vtt2;
      for (const auto &tt : vtt)
        {
          temporary_terms_t tt2;
          for (const auto &it : tt)
            tt2.insert(f(it));
          vtt2.push_back(tt2);
        }
      return vtt2;
    };

  for (const auto &it : m.v_temporary_terms)
    v_temporary_terms.push_back(convert_vector_tt(it));
  for (const auto &it : m.v_temporary_terms_local)
    v_temporary_terms_local.push_back(convert_vector_tt(it));

  for (const auto &it : m.first_chain_rule_derivatives)
    first_chain_rule_derivatives[it.first] = f(it.second);

  for (const auto &it : m.equation_type_and_normalized_equation)
    equation_type_and_normalized_equation.push_back(make_pair(it.first, f(it.second)));

  for (const auto &it : m.blocks_derivatives)
    {
      block_derivatives_equation_variable_laglead_nodeid_t v;
      for (const auto &it2 : it)
        v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
      blocks_derivatives.push_back(v);
    }

  for (const auto &it : m.dynamic_jacobian)
    dynamic_jacobian[it.first] = f(it.second);

  auto convert_derivative_t = [this,f](derivative_t dt)
    {
      derivative_t dt2;
      for (const auto &it : dt)
        dt2[it.first] = f(it.second);
      return dt2;
    };
  for (const auto &it : m.derivative_endo)
    derivative_endo.push_back(convert_derivative_t(it));
  for (const auto &it : m.derivative_other_endo)
    derivative_other_endo.push_back(convert_derivative_t(it));
  for (const auto &it : m.derivative_exo)
    derivative_exo.push_back(convert_derivative_t(it));
  for (const auto &it : m.derivative_exo_det)
    derivative_exo_det.push_back(convert_derivative_t(it));
}

89
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
90
                         NumericalConstants &num_constants_arg,
91
                         ExternalFunctionsTable &external_functions_table_arg) :
92
  ModelTree{symbol_table_arg, num_constants_arg, external_functions_table_arg}
93
94
{
}
95

96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
StaticModel::StaticModel(const StaticModel &m) :
  ModelTree {m},
  v_temporary_terms_inuse {m.v_temporary_terms_inuse},
  map_idx {m.map_idx},
  map_idx2 {m.map_idx2},
  global_temporary_terms {m.global_temporary_terms},
  block_type_firstequation_size_mfs {m.block_type_firstequation_size_mfs},
  blocks_linear {m.blocks_linear},
  other_endo_block {m.other_endo_block},
  exo_block {m.exo_block},
  exo_det_block {m.exo_det_block},
  block_col_type {m.block_col_type},
  variable_block_lead_lag {m.variable_block_lead_lag},
  equation_block {m.equation_block},
  endo_max_leadlag_block {m.endo_max_leadlag_block},
  other_endo_max_leadlag_block {m.other_endo_max_leadlag_block},
  exo_max_leadlag_block {m.exo_max_leadlag_block},
  exo_det_max_leadlag_block {m.exo_det_max_leadlag_block},
  max_leadlag_block {m.max_leadlag_block}
{
  copyHelper(m);
}

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

  v_temporary_terms.clear();
  v_temporary_terms_local.clear();

  v_temporary_terms_inuse = m.v_temporary_terms_inuse;

  first_chain_rule_derivatives.clear();

  map_idx = m.map_idx;
  map_idx2 = m.map_idx2;
  global_temporary_terms = m.global_temporary_terms;

  equation_type_and_normalized_equation.clear();

  block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;

  blocks_derivatives.clear();
  dynamic_jacobian.clear();

  blocks_linear = m.blocks_linear;

  derivative_endo.clear();
  derivative_other_endo.clear();
  derivative_exo.clear();
  derivative_exo_det.clear();

  other_endo_block = m.other_endo_block;
  exo_block = m.exo_block;
  exo_det_block = m.exo_det_block;
  block_col_type = m.block_col_type;
  variable_block_lead_lag = m.variable_block_lead_lag;
  equation_block = m.equation_block;
  endo_max_leadlag_block = m.endo_max_leadlag_block;
  other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
  exo_max_leadlag_block = m.exo_max_leadlag_block;
  exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
  max_leadlag_block = m.max_leadlag_block;

  copyHelper(m);

  return *this;
}

166
void
167
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
168
{
169
  auto it = first_derivatives.find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
170
  if (it != first_derivatives.end())
171
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
172
173
174
  else
    {
      FLDZ_ fldz;
175
      fldz.write(code_file, instruction_number);
176
177
    }
}
sebastien's avatar
sebastien committed
178

179
void
180
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
181
{
182
  auto it = first_chain_rule_derivatives.find({ eqr, { varr, lag } });
183
  if (it != first_chain_rule_derivatives.end())
184
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
185
186
187
  else
    {
      FLDZ_ fldz;
188
      fldz.write(code_file, instruction_number);
189
190
191
192
193
194
    }
}

void
StaticModel::computeTemporaryTermsOrdered()
{
195
  map<expr_t, pair<int, int>> first_occurence;
196
  map<expr_t, int> reference_count;
197
  BinaryOpNode *eq_node;
198
199
  first_derivatives_t::const_iterator it;
  first_chain_rule_derivatives_t::const_iterator it_chr;
200
201
202
203
204
  ostringstream tmp_s;
  v_temporary_terms.clear();
  map_idx.clear();

  unsigned int nb_blocks = getNbBlocks();
205
206
  v_temporary_terms = vector< vector<temporary_terms_t>>(nb_blocks);
  v_temporary_terms_local = vector< vector<temporary_terms_t>>(nb_blocks);
207

208
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
209

210
211
  map_idx2 = vector<map_idx_t>(nb_blocks);

212
  temporary_terms.clear();
213
214
215

  //local temporay terms
  for (unsigned int block = 0; block < nb_blocks; block++)
216
    {
217
218
      map<expr_t, int> reference_count_local;
      reference_count_local.clear();
219
      map<expr_t, pair<int, int>> first_occurence_local;
220
221
222
      first_occurence_local.clear();
      temporary_terms_t temporary_terms_l;
      temporary_terms_l.clear();
223

224
225
226
227
228
229
230
231
232
233
      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
234
            {
235
236
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  i);
237
238
            }
        }
239
240
241
242
243
244
245
246
247
      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]);
248
    }
249
250
251

  // global temporay terms
  for (unsigned int block = 0; block < nb_blocks; block++)
252
    {
253
254
255
256
257
258
      // 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++)
259
        {
260
261
262
          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
          else
263
            {
264
265
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
266
267
            }
        }
268
      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
269
        {
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
          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
287
            {
288
289
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
              eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
290
291
            }
        }
292
293
294
295
296
297
      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++)
298
        for (auto it = v_temporary_terms[block][i].begin();
299
300
301
             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;
302
    }
303
  computeTemporaryTermsMapping(temporary_terms, map_idx);
304
305
306
}

void
307
StaticModel::computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx)
308
{
309
  // Add a mapping form node ID to temporary terms order
310
  int j = 0;
311
312
  for (auto temporary_term : temporary_terms)
    map_idx[temporary_term->idx] = j++;
313
314
315
}

void
316
StaticModel::writeModelEquationsOrdered_M(const string &basename) const
317
318
319
{
  string tmp_s, sps;
  ostringstream tmp_output, tmp1_output, global_output;
320
  expr_t lhs = nullptr, rhs = nullptr;
321
  BinaryOpNode *eq_node;
322
  map<expr_t, int> reference_count;
323
  temporary_terms_t local_temporary_terms;
324
325
  ofstream  output;
  vector<int> feedback_variables;
326
  deriv_node_temp_terms_t tef_terms;
327
  ExprNodeOutputType local_output_type;
328

329
  local_output_type = ExprNodeOutputType::matlabStaticModelSparse;
330
  if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
331
    local_temporary_terms = temporary_terms;
332
333
334
335
336
337
338
339
340
341
342
343
344
345

  //----------------------------------------------------------------------
  //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("");
346
      tmp1_output << packageDir(basename + ".block") << "/static_" << block+1 << ".m";
347
      output.open(tmp1_output.str(), ios::out | ios::binary);
348
349
350
351
352
353
354
      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)
355
        output << "function y = static_" << block+1 << "(y, x, params)\n";
356
      else
357
        output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)\n";
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380

      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)
381
        output << " g1 = spalloc("  << block_mfs << ", " << block_mfs << ", " << derivative_endo[block].size() << ");" << endl;
382
383
384
385

      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
386
387
          for (int it : v_temporary_terms_inuse[block])
            tmp_output << " T" << it;
388
389
390
391
392
393
394
          output << "  global" << tmp_output.str() << ";\n";
        }

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

      // The equations
Houtan Bastani's avatar
Houtan Bastani committed
395
      temporary_terms_idxs_t temporary_terms_idxs;
396
397
398
399
      for (unsigned int i = 0; i < block_size; i++)
        {
          if (!global_temporary_terms)
            local_temporary_terms = v_temporary_terms[block][i];
400
          temporary_terms_t tt2;
401
402
403
404
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
405
              for (auto it : v_temporary_terms[block][i])
406
                {
407
                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
408
                    it->writeExternalFunctionOutput(output, local_output_type, tt2, {}, tef_terms);
409

410
                  output << "  " <<  sps;
411
                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
412
                  output << " = ";
413
                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
414
                  // Insert current node into tt2
415
                  tt2.insert(it);
416
417
418
419
420
421
422
                  output << ";" << endl;
                }
            }

          int variable_ID = getBlockVariableID(block, i);
          int equation_ID = getBlockEquationID(block, i);
          EquationType equ_type = getBlockEquationType(block, i);
423
          string sModel = symbol_table.getName(symbol_table.getID(SymbolType::endogenous, variable_ID));
424
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
425
426
427
          lhs = eq_node->get_arg1();
          rhs = eq_node->get_arg2();
          tmp_output.str("");
428
          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
429
430
431
432
433
434
435
436
437
438
439
440
          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 << " = ";
441
                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
442
443
444
445
446
447
448
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  output << "%" << tmp_output.str();
                  output << " = ";
                  if (isBlockEquationRenormalized(block, i))
                    {
449
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
450
451
                      output << "\n  ";
                      tmp_output.str("");
452
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
453
454
                      lhs = eq_node->get_arg1();
                      rhs = eq_node->get_arg2();
455
                      lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
456
                      output << " = ";
457
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
458
459
460
461
                    }
                }
              else
                {
462
                  cerr << "Type mismatch for equation " << equation_ID+1  << "\n";
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
                  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 << ") - (";
482
              rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
483
484
485
486
487
488
489
490
491
492
493
494
495
              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:
496
          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
497
498
499
500
501
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
502
              expr_t id = it->second.second;
503
              output << "    g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
504
              id->writeOutput(output, local_output_type, local_temporary_terms, {});
505
              output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
506
507
508
509
510
511
512
513
                     << "(" << 0
                     << ") " << varr+1
                     << ", equation=" << eqr+1 << endl;
            }
          break;
        default:
          break;
        }
514
      output << "end" << endl;
515
516
517
      output.close();
    }
}
518
519

void
520
StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx) const
521
522
523
524
{

  ostringstream tmp_output;
  ofstream code_file;
525
  unsigned int instruction_number = 0;
526
527
  bool file_open = false;

528
529
530
  boost::filesystem::create_directories(basename + "/model/bytecode");

  string main_name = basename + "/model/bytecode/static.cod";
531
  code_file.open(main_name, ios::out | ios::binary | ios::ate);
532
533
  if (!code_file.is_open())
    {
534
      cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
535
536
537
538
539
      exit(EXIT_FAILURE);
    }
  int count_u;
  int u_count_int = 0;

540
  Write_Inf_To_Bin_File(basename + "/model/bytecode/static.bin", u_count_int, file_open, false, symbol_table.endo_nbr());
541
542
543
  file_open = true;

  //Temporary variables declaration
544
545
  FDIMST_ fdimst(temporary_terms.size());
  fdimst.write(code_file, instruction_number);
546
547
548
549
550
551
552
553
554
555
  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,
556
557
                           u_count_int,
                           symbol_table.endo_nbr()
558
                           );
559
  fbeginblock.write(code_file, instruction_number);
560
561
562

  // Add a mapping form node ID to temporary terms order
  int j = 0;
563
564
  for (auto temporary_term : temporary_terms)
    map_idx[temporary_term->idx] = j++;
565
  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
566

567
  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
568
569

  FENDEQU_ fendequ;
570
  fendequ.write(code_file, instruction_number);
571

572
573
574
575
576
577
  // 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;

578
  vector<vector<pair<int, int>>> derivatives;
579
580
  derivatives.resize(symbol_table.endo_nbr());
  count_u = symbol_table.endo_nbr();
581
  for (const auto & first_derivative : first_derivatives)
582
    {
583
      int deriv_id = first_derivative.first.second;
584
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
585
        {
586
587
          expr_t d1 = first_derivative.second;
          unsigned int eq = first_derivative.first.first;
588
589
          int symb = getSymbIDByDerivID(deriv_id);
          unsigned int var = symbol_table.getTypeSpecificID(symb);
590
          FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
591
          fnumexpr.write(code_file, instruction_number);
592
593
          if (!derivatives[eq].size())
            derivatives[eq].clear();
594
          derivatives[eq].emplace_back(var, count_u);
595

596
          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
597
598

          FSTPSU_ fstpsu(count_u);
599
          fstpsu.write(code_file, instruction_number);
600
601
602
603
604
605
          count_u++;
        }
    }
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
    {
      FLDR_ fldr(i);
606
      fldr.write(code_file, instruction_number);
607
      if (derivatives[i].size())
608
        {
609
          for (vector<pair<int, int>>::const_iterator it = derivatives[i].begin();
610
               it != derivatives[i].end(); it++)
611
            {
612
613
              FLDSU_ fldsu(it->second);
              fldsu.write(code_file, instruction_number);
614
              FLDSV_ fldsv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(it->first)};
615
              fldsv.write(code_file, instruction_number);
616
              FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
617
              fbinary.write(code_file, instruction_number);
618
619
              if (it != derivatives[i].begin())
                {
620
                  FBINARY_ fbinary{static_cast<int>(BinaryOpcode::plus)};
621
622
                  fbinary.write(code_file, instruction_number);
                }
623
            }
624
          FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
625
          fbinary.write(code_file, instruction_number);
626
627
        }
      FSTPSU_ fstpsu(i);
628
      fstpsu.write(code_file, instruction_number);
629
    }
630
631
632
633
634
635
636
637
638
639
  // 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);
640
  prev_instruction_number = instruction_number;
641
642
643
644
645
646
647

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

  // The Jacobian if we have to solve the block determinsitic bloc
648
  for (const auto & first_derivative : first_derivatives)
649
    {
650
      int deriv_id = first_derivative.first.second;
651
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
652
        {
653
654
          expr_t d1 = first_derivative.second;
          unsigned int eq = first_derivative.first.first;
655
656
657
658
659
660
          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();
661
          derivatives[eq].emplace_back(var, count_u);
662
663

          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
664
          FSTPG2_ fstpg2(eq, var);
665
666
667
668
669
670
671
672
673
674
675
          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);

676
  FENDBLOCK_ fendblock;
677
  fendblock.write(code_file, instruction_number);
678
  FEND_ fend;
679
  fend.write(code_file, instruction_number);
680
681
682
683
  code_file.close();
}

void
684
StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const
685
686
{
  struct Uff_l
687
  {
688
689
690
    int u, var, lag;
    Uff_l *pNext;
  };
691

692
693
694
695
696
697
698
699
700
  struct Uff
  {
    Uff_l *Ufl, *Ufl_First;
  };

  int i, v;
  string tmp_s;
  ostringstream tmp_output;
  ofstream code_file;
701
  unsigned int instruction_number = 0;
702
  expr_t lhs = nullptr, rhs = nullptr;
703
704
  BinaryOpNode *eq_node;
  Uff Uf[symbol_table.endo_nbr()];
705
  map<expr_t, int> reference_count;
706
  vector<int> feedback_variables;
707
  deriv_node_temp_terms_t tef_terms;
708
709
  bool file_open = false;

710
711
712
  boost::filesystem::create_directories(basename + "/model/bytecode");

  string main_name = basename + "/model/bytecode/static.cod";
713
  code_file.open(main_name, ios::out | ios::binary | ios::ate);
714
715
  if (!code_file.is_open())
    {
716
      cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
717
718
719
720
      exit(EXIT_FAILURE);
    }
  //Temporary variables declaration

721
722
  FDIMST_ fdimst(temporary_terms.size());
  fdimst.write(code_file, instruction_number);
723
724
725
726
727
728
729

  for (unsigned int block = 0; block < getNbBlocks(); block++)
    {
      feedback_variables.clear();
      if (block > 0)
        {
          FENDBLOCK_ fendblock;
730
          fendblock.write(code_file, instruction_number);
731
732
733
734
735
736
737
738
739
740
741
        }
      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)
        {
742
          Write_Inf_To_Bin_File_Block(basename, block, u_count_int, file_open);
743
744
745
746
747
748
749
750
751
752
753
754
755
          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,
756
                               u_count_int,
757
                               /*symbol_table.endo_nbr()*/ block_size
758
                               );
759

760
      fbeginblock.write(code_file, instruction_number);
761

762
763
764
765
766
767
      // 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;

768
769
770
      for (i = 0; i < (int) block_size; i++)
        {
          //The Temporary terms
771
          temporary_terms_t tt2;
772
773
774
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
775
              for (auto it : v_temporary_terms[block][i])
776
                {
777
                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
778
                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
779

780
                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find(it->idx)->second));
781
                  fnumexpr.write(code_file, instruction_number);
782
783
                  it->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
                  FSTPST_ fstpst((int)(map_idx.find(it->idx)->second));
784
                  fstpst.write(code_file, instruction_number);
785
                  // Insert current node into tt2
786
                  tt2.insert(it);
787
788
789
                }
            }

790
          // The equations
791
792
793
794
795
796
797
798
          int variable_ID, equation_ID;
          EquationType equ_type;
          switch (simulation_type)
            {
            evaluation:
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
              equ_type = getBlockEquationType(block, i);
799
800
              {
                FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
801
                fnumexpr.write(code_file, instruction_number);
802
              }
803
804
              if (equ_type == E_EVALUATE)
                {
805
                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
806
807
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
808
809
                  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);
810
811
812
                }
              else if (equ_type == E_EVALUATE_S)
                {
813
                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
814
815
                  lhs = eq_node->get_arg1();
                  rhs = eq_node->get_arg2();
816
817
                  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);
818
819
820
821
822
823
824
825
826
                }
              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);
827
              Uf[equation_ID].Ufl = nullptr;
828
829
830
              goto end;
            default:
            end:
831
              FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
832
              fnumexpr.write(code_file, instruction_number);
833
              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
834
835
              lhs = eq_node->get_arg1();
              rhs = eq_node->get_arg2();
836
837
              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);
838

839
              FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
840
              fbinary.write(code_file, instruction_number);
841
842

              FSTPR_ fstpr(i - block_recursive);
843
              fstpr.write(code_file, instruction_number);
844
845
846
            }
        }
      FENDEQU_ fendequ;
847
      fendequ.write(code_file, instruction_number);
848

849
850
851
852
853
854
855
856
      // 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:
857
858
              {
                FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
859
                fnumexpr.write(code_file, instruction_number);
860
              }
861
              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx, temporary_terms);
862
              {
863
                FSTPG_ fstpg(0);
864
                fstpg.write(code_file, instruction_number);
865
              }
866
              break;
867

868
869
870
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              count_u = feedback_variables.size();
871
              for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
872
873
874
875
876
                {
                  unsigned int eq = it->first.first;
                  unsigned int var = it->first.second;
                  unsigned int eqr = getBlockEquationID(block, eq);
                  unsigned int varr = getBlockVariableID(block, var);
877
                  if (eq >= block_recursive && var >= block_recursive)
878
879
880
881
882
883
884
885
886
887
888
                    {
                      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;
                        }
889
                      Uf[eqr].Ufl->pNext = nullptr;
890
891
                      Uf[eqr].Ufl->u = count_u;
                      Uf[eqr].Ufl->var = varr;
892
                      FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
893
                      fnumexpr.write(code_file, instruction_number);
894
                      compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx, temporary_terms);
895
                      FSTPSU_ fstpsu(count_u);
896
                      fstpsu.write(code_file, instruction_number);
897
898
899
900
901
902
903
904
                      count_u++;
                    }
                }
              for (i = 0; i < (int) block_size; i++)
                {
                  if (i >= (int) block_recursive)
                    {
                      FLDR_ fldr(i-block_recursive);
905
                      fldr.write(code_file, instruction_number);
906
907

                      FLDZ_ fldz;
908
                      fldz.write(code_file, instruction_number);
909
910
911
912
913

                      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);
914
                          fldsu.write(code_file, instruction_number);
915
                          FLDSV_ fldsv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(Uf[v].Ufl->var)};
916
                          fldsv.write(code_file, instruction_number);
917

918
                          FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
919
                          fbinary.write(code_file, instruction_number);
920
921

                          FCUML_ fcuml;
922
                          fcuml.write(code_file, instruction_number);
923
924
925
926
927
928
929
930
                        }
                      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;
                        }
931
                      FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
932
                      fbinary.write(code_file, instruction_number);
933
934

                      FSTPSU_ fstpsu(i - block_recursive);
935
                      fstpsu.write(code_file, instruction_number);
936
937
938
939
940
941
942
943

                    }
                }
              break;
            default:
              break;
            }
        }
944
945
946
947
948
949
950
951
952
953
954

      // 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);
955
      prev_instruction_number = instruction_number;
956
957
958
959
960

      temporary_terms_t tt2;
      tt2.clear();
      temporary_terms_t tt3;
      tt3.clear();
961
      deriv_node_temp_terms_t tef_terms2;
962
963
964
965
966

      for (i = 0; i < (int) block_size; i++)
        {
          if (v_temporary_terms_local[block].size())
            {
967
              for (auto it = v_temporary_terms_local[block][i].begin();
968
969
                   it != v_temporary_terms_local[block][i].end(); it++)
                {
970
                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != nullptr)
971
                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
972

Stéphane Adjemian's avatar
Stéphane Adjemian committed
973
                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find((*it)->idx)->second));
974
                  fnumexpr.write(code_file, instruction_number);
975
976
977

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

Stéphane Adjemian's avatar
Stéphane Adjemian committed
978
                  FSTPST_ fstpst((int)(map_idx2[block].find((*it)->idx)->second));
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
                  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)
                {