DynamicModel.cc 292 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
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
#include <iostream>
sebastien's avatar
sebastien committed
21
#include <cmath>
22
#include <cstdlib>
23
#include <cassert>
24
25
#include <cstdio>
#include <cerrno>
26
#include <algorithm>
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
27
#include <iterator>
sebastien's avatar
sebastien committed
28
29
30
31
32
33
34
35
36
37
38
39
#include "DynamicModel.hh"

// For mkdir() and chdir()
#ifdef _WIN32
# include <direct.h>
#else
# include <unistd.h>
# include <sys/stat.h>
# include <sys/types.h>
#endif

DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
40
41
42
                           NumericalConstants &num_constants_arg,
                           ExternalFunctionsTable &external_functions_table_arg) :
  ModelTree(symbol_table_arg, num_constants_arg, external_functions_table_arg),
43
44
45
46
  max_lag(0), max_lead(0),
  max_endo_lag(0), max_endo_lead(0),
  max_exo_lag(0), max_exo_lead(0),
  max_exo_det_lag(0), max_exo_det_lead(0),
47
48
49
50
  max_lag_orig(0), max_lead_orig(0),
  max_endo_lag_orig(0), max_endo_lead_orig(0),
  max_exo_lag_orig(0), max_exo_lead_orig(0),
  max_exo_det_lag_orig(0), max_exo_det_lead_orig(0),
51
  dynJacobianColsNbr(0),
52
  global_temporary_terms(true)
sebastien's avatar
sebastien committed
53
54
55
{
}

sebastien's avatar
sebastien committed
56
57
VariableNode *
DynamicModel::AddVariable(int symb_id, int lag)
sebastien's avatar
sebastien committed
58
{
sebastien's avatar
sebastien committed
59
  return AddVariableInternal(symb_id, lag);
sebastien's avatar
sebastien committed
60
61
}

sebastien's avatar
sebastien committed
62
void
63
DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
64
{
65
  auto it = first_derivatives.find({ eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag) });
66
  if (it != first_derivatives.end())
67
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
68
69
70
  else
    {
      FLDZ_ fldz;
71
      fldz.write(code_file, instruction_number);
72
73
    }
}
74
75

void
76
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
77
{
78
  auto it = first_chain_rule_derivatives.find({ eqr, { varr, lag } });
79
  if (it != first_chain_rule_derivatives.end())
80
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
81
  else
82
83
    {
      FLDZ_ fldz;
84
      fldz.write(code_file, instruction_number);
85
    }
86
87
}

sebastien's avatar
sebastien committed
88
void
89
DynamicModel::computeTemporaryTermsOrdered()
sebastien's avatar
sebastien committed
90
{
91
  map<expr_t, pair<int, int>> first_occurence;
92
  map<expr_t, int> reference_count;
sebastien's avatar
sebastien committed
93
  BinaryOpNode *eq_node;
94
95
  first_derivatives_t::const_iterator it;
  first_chain_rule_derivatives_t::const_iterator it_chr;
sebastien's avatar
sebastien committed
96
  ostringstream tmp_s;
97
98
  v_temporary_terms.clear();
  map_idx.clear();
sebastien's avatar
sebastien committed
99

100
  unsigned int nb_blocks = getNbBlocks();
101
  v_temporary_terms = vector<vector<temporary_terms_t>>(nb_blocks);
102
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
sebastien's avatar
sebastien committed
103
  temporary_terms.clear();
104

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

204
205
206
207
208
void
DynamicModel::computeTemporaryTermsMapping()
{
  // Add a mapping form node ID to temporary terms order
  int j = 0;
209
210
  for (auto temporary_term : temporary_terms)
    map_idx[temporary_term->idx] = j++;
211
212
}

sebastien's avatar
sebastien committed
213
void
214
DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
215
216
217
{
  string tmp_s, sps;
  ostringstream tmp_output, tmp1_output, global_output;
218
  expr_t lhs = nullptr, rhs = nullptr;
219
  BinaryOpNode *eq_node;
220
221
  ostringstream Ufoss;
  vector<string> Uf(symbol_table.endo_nbr(), "");
222
  map<expr_t, int> reference_count;
223
  temporary_terms_t local_temporary_terms;
224
  ofstream  output;
225
  int nze, nze_exo, nze_exo_det, nze_other_endo;
226
227
  vector<int> feedback_variables;
  ExprNodeOutputType local_output_type;
228
  Ufoss.str("");
sebastien's avatar
sebastien committed
229

Sébastien Villemot's avatar
Sébastien Villemot committed
230
  local_output_type = oMatlabDynamicModelSparse;
231
  if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
232
    local_temporary_terms = temporary_terms;
233
234
235
236
237
238
239
240
241
242
243
244

  //----------------------------------------------------------------------
  //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();
      nze_other_endo = derivative_other_endo[block].size();
      nze_exo = derivative_exo[block].size();
245
      nze_exo_det = derivative_exo_det[block].size();
246
247
248
249
      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;
250
      deriv_node_temp_terms_t tef_terms;
Sébastien Villemot's avatar
Sébastien Villemot committed
251
      local_output_type = oMatlabDynamicModelSparse;
252
      if (global_temporary_terms)
Sébastien Villemot's avatar
Sébastien Villemot committed
253
        local_temporary_terms = temporary_terms;
254

255
256
      int prev_lag;
      unsigned int prev_var, count_col, count_col_endo, count_col_exo, count_col_exo_det, count_col_other_endo;
257
      map<pair<int, pair<int, int>>, expr_t> tmp_block_endo_derivative;
258
      for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
259
        tmp_block_endo_derivative[{ it->second.first, { it->first.second, it->first.first } }] = it->second.second;
260
261
262
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_endo = 0;
263
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
264
265
266
267
268
269
270
271
272
273
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_endo++;
            }
        }
274
      map<pair<int, pair<int, int>>, expr_t> tmp_block_exo_derivative;
275
      for (auto it = derivative_exo[block].begin(); it != (derivative_exo[block]).end(); it++)
276
        tmp_block_exo_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
277
278
279
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_exo = 0;
280
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
281
282
283
284
285
286
287
288
289
290
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_exo++;
            }
        }
291
      map<pair<int, pair<int, int>>, expr_t> tmp_block_exo_det_derivative;
292
      for (auto it = derivative_exo_det[block].begin(); it != (derivative_exo_det[block]).end(); it++)
293
        tmp_block_exo_det_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
294
295
296
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_exo_det = 0;
297
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
298
299
300
301
302
303
304
305
306
307
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_exo_det++;
            }
        }
308
      map<pair<int, pair<int, int>>, expr_t> tmp_block_other_endo_derivative;
309
      for (auto it = derivative_other_endo[block].begin(); it != (derivative_other_endo[block]).end(); it++)
310
        tmp_block_other_endo_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
311
312
313
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_other_endo = 0;
314
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_other_endo_derivative.begin(); it != tmp_block_other_endo_derivative.end(); it++)
315
316
317
318
319
320
321
322
323
324
325
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          if (var != prev_var || lag != prev_lag)
            {
              prev_var = var;
              prev_lag = lag;
              count_col_other_endo++;
            }
        }

326
327
328
329
330
331
332
333
334
335
336
      tmp1_output.str("");
      tmp1_output << dynamic_basename << "_" << block+1 << ".m";
      output.open(tmp1_output.str().c_str(), ios::out | ios::binary);
      output << "%\n";
      output << "% " << tmp1_output.str() << " : Computes dynamic 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)
        {
337
          output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)\n";
338
339
        }
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
340
        output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
341
      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
342
        output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
343
      else
344
        output << "function [residual, y, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)\n";
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
      BlockType block_type;
      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        block_type = SIMULTAN;
      else 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;
      //The Temporary terms
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        {
          output << "  if(jacobian_eval)\n";
370
371
          output << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
372
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
373
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
374
375
376
377
378
          output << "  end;\n";
        }
      else
        {
          output << "  if(jacobian_eval)\n";
379
380
          output << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
381
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
382
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
383
384
385
          output << "  else\n";
          if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
            {
386
387
388
              output << "    g1 = spalloc(" << block_mfs << "*Periods, "
                     << block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
                     << ", " << nze << "*Periods);\n";
389
            }
ferhat's avatar
ferhat committed
390
          else
391
392
393
394
395
396
            {
              output << "    g1 = spalloc(" << block_mfs
                     << ", " << block_mfs << ", " << nze << ");\n";
            }
          output << "  end;\n";
        }
397

398
399
400
401
      output << "  g2=0;g3=0;\n";
      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
402
403
          for (int it : v_temporary_terms_inuse[block])
            tmp_output << " T" << it;
404
405
406
407
          output << "  global" << tmp_output.str() << ";\n";
        }
      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        {
408
          temporary_terms_t tt2;
409
410
411
412
413
414
415
          tt2.clear();
          for (int i = 0; i < (int) block_size; i++)
            {
              if (v_temporary_terms[block][i].size() && global_temporary_terms)
                {
                  output << "  " << "% //Temporary variables initialization" << endl
                         << "  " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
416
                  for (auto it : v_temporary_terms[block][i])
417
418
                    {
                      output << "  ";
419
                      // In the following, "Static" is used to avoid getting the "(it_)" subscripting
420
                      it->writeOutput(output, oMatlabStaticModelSparse, local_temporary_terms, {});
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
                      output << " = T_zeros;" << endl;
                    }
                }
            }
        }
      if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
        output << "  residual=zeros(" << block_mfs << ",1);\n";
      else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        output << "  residual=zeros(" << block_mfs << ",y_kmin+periods);\n";
      if (simulation_type == EVALUATE_BACKWARD)
        output << "  for it_ = (y_kmin+periods):y_kmin+1\n";
      if (simulation_type == EVALUATE_FORWARD)
        output << "  for it_ = y_kmin+1:(y_kmin+periods)\n";

      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        {
          output << "  b = zeros(periods*y_size,1);" << endl
                 << "  for it_ = y_kmin+1:(periods+y_kmin)" << endl
                 << "    Per_y_=it_*y_size;" << endl
                 << "    Per_J_=(it_-y_kmin-1)*y_size;" << endl
                 << "    Per_K_=(it_-1)*y_size;" << endl;
          sps = "  ";
        }
      else
        if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
          sps = "  ";
        else
          sps = "";
      // The equations
Houtan Bastani's avatar
Houtan Bastani committed
450
      temporary_terms_idxs_t temporary_terms_idxs;
451
452
      for (unsigned int i = 0; i < block_size; i++)
        {
453
          temporary_terms_t tt2;
454
455
456
457
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
458
              for (auto it : v_temporary_terms[block][i])
459
                {
460
                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
461
                    it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
462

463
                  output << "  " <<  sps;
464
                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
465
                  output << " = ";
466
                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
467
                  // Insert current node into tt2
468
                  tt2.insert(it);
469
470
471
472
473
474
475
476
                  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));
477
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
478
479
480
          lhs = eq_node->get_arg1();
          rhs = eq_node->get_arg2();
          tmp_output.str("");
481
          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
482
483
484
485
486
487
488
489
490
491
492
493
          switch (simulation_type)
            {
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
            evaluation:     if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
                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 << " = ";
494
                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
495
496
497
498
499
500
501
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  output << "%" << tmp_output.str();
                  output << " = ";
                  if (isBlockEquationRenormalized(block, i))
                    {
502
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
503
504
                      output << "\n    ";
                      tmp_output.str("");
505
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
506
507
                      lhs = eq_node->get_arg1();
                      rhs = eq_node->get_arg2();
508
                      lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
509
                      output << " = ";
510
                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
511
512
513
514
                    }
                }
              else
                {
515
                  cerr << "Type mismatch for equation " << equation_ID+1  << "\n";
516
517
518
519
520
521
522
523
524
525
526
527
                  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
528
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
529
530
531
532
533
534
535
536
              output << "  " << "residual(" << i+1-block_recursive << ") = (";
              goto end;
            case SOLVE_TWO_BOUNDARIES_COMPLETE:
            case SOLVE_TWO_BOUNDARIES_SIMPLE:
              if (i < block_recursive)
                goto evaluation;
              feedback_variables.push_back(variable_ID);
              output << "    % equation " << equation_ID+1 << " variable : " << sModel
537
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
538
              Ufoss << "    b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
539
              Uf[equation_ID] += Ufoss.str();
540
              Ufoss.str("");
541
542
543
544
545
546
              output << "    residual(" << i+1-block_recursive << ", it_) = (";
              goto end;
            default:
            end:
              output << tmp_output.str();
              output << ") - (";
547
              rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
548
              output << ");\n";
sebastien's avatar
sebastien committed
549
#ifdef CONDITION
550
551
              if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
                output << "  condition(" << i+1 << ")=0;\n";
sebastien's avatar
sebastien committed
552
#endif
553
554
555
556
            }
        }
      // The Jacobian if we have to solve the block
      if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
557
        output << "  " << sps << "% Jacobian  " << endl << "    if jacobian_eval" << endl;
558
559
560
561
      else
        if (simulation_type == SOLVE_BACKWARD_SIMPLE   || simulation_type == SOLVE_FORWARD_SIMPLE
            || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
          output << "  % Jacobian  " << endl << "  if jacobian_eval" << endl;
sebastien's avatar
sebastien committed
562
        else
563
          output << "    % Jacobian  " << endl << "    if jacobian_eval" << endl;
564
565
566
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
567
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
568
        {
569
570
571
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
572
573
          int eqr = getBlockEquationID(block, eq);
          int varr = getBlockVariableID(block, var);
574
          if (var != prev_var || lag != prev_lag)
575
            {
576
577
578
579
              prev_var = var;
              prev_lag = lag;
              count_col++;
            }
580

581
          expr_t id = it->second;
582

583
          output << "      g1(" << eq+1 << ", " << count_col << ") = ";
584
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
585
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
586
                 << "(" << lag
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
587
588
                 << ") " << varr+1 << ", " << var+1
                 << ", equation=" << eqr+1 << ", " << eq+1 << endl;
589
590
591
592
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
593
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
594
595
596
597
598
599
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
600
            {
601
602
603
              prev_var = var;
              prev_lag = lag;
              count_col++;
604
            }
605
606
          expr_t id = it->second;
          output << "      g1_x(" << eqr+1 << ", " << count_col << ") = ";
607
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
608
609
610
611
612
613
614
615
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
616
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_det_derivative.begin(); it != tmp_block_exo_det_derivative.end(); it++)
617
618
619
620
621
622
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
623
            {
624
625
626
              prev_var = var;
              prev_lag = lag;
              count_col++;
627
            }
628
629
          expr_t id = it->second;
          output << "      g1_xd(" << eqr+1 << ", " << count_col << ") = ";
630
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
631
632
633
634
635
636
637
638
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
639
      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_other_endo_derivative.begin(); it != tmp_block_other_endo_derivative.end(); it++)
640
641
642
643
644
645
        {
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          int eqr = getBlockInitialEquationID(block, eq);
          if (var != prev_var || lag != prev_lag)
646
            {
647
648
649
              prev_var = var;
              prev_lag = lag;
              count_col++;
650
            }
651
652
          expr_t id = it->second;

653
          output << "      g1_o(" << eqr+1 << ", " << /*var+1+(lag+block_max_lag)*block_size*/ count_col << ") = ";
654
          id->writeOutput(output, local_output_type, local_temporary_terms, {});
655
656
657
658
659
660
661
662
663
664
665
666
667
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      output << "      varargout{1}=g1_x;\n";
      output << "      varargout{2}=g1_xd;\n";
      output << "      varargout{3}=g1_o;\n";

      switch (simulation_type)
        {
        case EVALUATE_FORWARD:
        case EVALUATE_BACKWARD:
668
669
670
671
672
673
674
675
          output << "    end;" << endl;
          output << "  end;" << endl;
          break;
        case SOLVE_BACKWARD_SIMPLE:
        case SOLVE_FORWARD_SIMPLE:
        case SOLVE_BACKWARD_COMPLETE:
        case SOLVE_FORWARD_COMPLETE:
          output << "  else" << endl;
676
          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
677
678
679
680
681
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
682
              expr_t id = it->second.second;
683
              int lag = it->second.first;
684
685
686
              if (lag == 0)
                {
                  output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
687
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
688
689
690
691
692
693
                  output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
                         << "(" << lag
                         << ") " << varr+1
                         << ", equation=" << eqr+1 << endl;
                }

694
695
696
697
698
            }
          output << "  end;\n";
          break;
        case SOLVE_TWO_BOUNDARIES_SIMPLE:
        case SOLVE_TWO_BOUNDARIES_COMPLETE:
699
          output << "    else" << endl;
700
          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
701
702
703
704
705
706
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
              ostringstream tmp_output;
707
              expr_t id = it->second.second;
708
              int lag = it->second.first;
709
              if (eq >= block_recursive && var >= block_recursive)
710
711
                {
                  if (lag == 0)
712
713
714
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+Per_K_)*y(it_, " << varr+1 << ")";
715
                  else if (lag == 1)
716
717
718
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+Per_y_)*y(it_+1, " << varr+1 << ")";
719
                  else if (lag > 0)
720
721
722
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")";
723
                  else
724
725
726
                    Ufoss << "+g1(" << eq+1-block_recursive
                          << "+Per_J_, " << var+1-block_recursive
                          << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")";
727
                  Uf[eqr] += Ufoss.str();
728
729
                  Ufoss.str("");

730
731
732
733
734
735
736
737
738
739
740
741
742
                  if (lag == 0)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+Per_K_) = ";
                  else if (lag == 1)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+Per_y_) = ";
                  else if (lag > 0)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+y_size*(it_+" << lag-1 << ")) = ";
                  else if (lag < 0)
                    tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                               << var+1-block_recursive << "+y_size*(it_" << lag-1 << ")) = ";
                  output << " " << tmp_output.str();
743
                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
744
745
746
747
748
                  output << ";";
                  output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
                         << "(" << lag << ") " << varr+1
                         << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
                }
749

sebastien's avatar
sebastien committed
750
#ifdef CONDITION
751
752
              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
              output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
sebastien's avatar
sebastien committed
753
#endif
754
755
756
757
            }
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i >= block_recursive)
758
                output << "  " << Uf[getBlockEquationID(block, i)] << ";\n";
sebastien's avatar
sebastien committed
759
#ifdef CONDITION
760
761
              output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
              output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
sebastien's avatar
sebastien committed
762
#endif
763
            }
sebastien's avatar
sebastien committed
764
#ifdef CONDITION
765
766
767
768
769
770
771
772
773
774
775
776
777
778
          for (m = 0; m <= ModelBlock->Block_List[block].Max_Lead+ModelBlock->Block_List[block].Max_Lag; m++)
            {
              k = m-ModelBlock->Block_List[block].Max_Lag;
              for (i = 0; i < ModelBlock->Block_List[block].IM_lead_lag[m].size; i++)
                {
                  unsigned int eq = ModelBlock->Block_List[block].IM_lead_lag[m].Equ_Index[i];
                  unsigned int var = ModelBlock->Block_List[block].IM_lead_lag[m].Var_Index[i];
                  unsigned int u = ModelBlock->Block_List[block].IM_lead_lag[m].u[i];
                  unsigned int eqr = ModelBlock->Block_List[block].IM_lead_lag[m].Equ[i];
                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
                }
            }
          for (i = 0; i < ModelBlock->Block_List[block].Size; i++)
            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
sebastien's avatar
sebastien committed
779
#endif
780
781
          output << "    end;" << endl;
          output << "  end;" << endl;
782
783
784
785
          break;
        default:
          break;
        }
786
      output << "end" << endl;
787
788
789
      output.close();
    }
}
sebastien's avatar
sebastien committed
790
791

void
792
DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basename, const map_idx_t &map_idx) const
793
{
794

795
796
  ostringstream tmp_output;
  ofstream code_file;
797
  unsigned int instruction_number = 0;
798
799
800
801
802
803
804
  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())
    {
805
      cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
806
807
808
809
810
811
812
813
814
815
816
817
818
      exit(EXIT_FAILURE);
    }

  int count_u;
  int u_count_int = 0;
  BlockSimulationType simulation_type;
  if ((max_endo_lag > 0) && (max_endo_lead > 0))
    simulation_type = SOLVE_TWO_BOUNDARIES_COMPLETE;
  else if ((max_endo_lag >= 0) && (max_endo_lead == 0))
    simulation_type = SOLVE_FORWARD_COMPLETE;
  else
    simulation_type = SOLVE_BACKWARD_COMPLETE;

819
  Write_Inf_To_Bin_File(file_name, u_count_int, file_open, simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE, symbol_table.endo_nbr());
820
821
822
823
  file_open = true;

  //Temporary variables declaration
  FDIMT_ fdimt(temporary_terms.size());
824
825
826
  fdimt.write(code_file, instruction_number);

  vector<unsigned int> exo, exo_det, other_endo;
827

828
  for (int i = 0; i < symbol_table.exo_det_nbr(); i++)
829
    exo_det.push_back(i);
830
  for (int i = 0; i < symbol_table.exo_nbr(); i++)
831
    exo.push_back(i);
832

833
834
  map<pair< int, pair<int, int>>, expr_t> first_derivatives_reordered_endo;
  map<pair< pair<int, int>, pair<int, int>>, expr_t>  first_derivatives_reordered_exo;
835
  for (const auto & first_derivative : first_derivatives)
836
    {
837
838
      int deriv_id = first_derivative.first.second;
      unsigned int eq = first_derivative.first.first;
839
840
841
842
      int symb = getSymbIDByDerivID(deriv_id);
      unsigned int var = symbol_table.getTypeSpecificID(symb);
      int lag = getLagByDerivID(deriv_id);
      if (getTypeByDerivID(deriv_id) == eEndogenous)
843
        first_derivatives_reordered_endo[{ lag, make_pair(var, eq) }] = first_derivative.second;
844
      else if (getTypeByDerivID(deriv_id) == eExogenous || getTypeByDerivID(deriv_id) == eExogenousDet)
845
        first_derivatives_reordered_exo[{ { lag, getTypeByDerivID(deriv_id) }, { var, eq } }] = first_derivative.second;
846
847
848
849
    }
  int prev_var = -1;
  int prev_lag = -999999999;
  int count_col_endo = 0;
850
  for (map<pair< int, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
851
852
853
854
       it != first_derivatives_reordered_endo.end(); it++)
    {
      int var = it->first.second.first;
      int lag = it->first.first;
855
      if (prev_var != var || prev_lag != lag)
856
857
858
859
860
861
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
    }
862
863
  prev_var = -1;
  prev_lag = -999999999;
864
  int prev_type = -1;
865
  int count_col_exo = 0;
866
  int count_col_det_exo = 0;
867

868
  for (map<pair< pair<int, int>, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_exo.begin();
869
870
871
       it != first_derivatives_reordered_exo.end(); it++)
    {
      int var = it->first.second.first;
872
873
874
      int lag = it->first.first.first;
      int type = it->first.first.second;
      if (prev_var != var || prev_lag != lag || prev_type != type)
875
876
877
        {
          prev_var = var;
          prev_lag = lag;
878
          prev_type = type;
879
880
881
882
          if (type == eExogenous)
            count_col_exo++;
          else if (type == eExogenousDet)
            count_col_det_exo++;
883
884
        }
    }
885

886
887
888
889
890
891
892
893
  FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
                           simulation_type,
                           0,
                           symbol_table.endo_nbr(),
                           variable_reordered,
                           equation_reordered,
                           false,
                           symbol_table.endo_nbr(),
894
895
                           max_endo_lag,
                           max_endo_lead,
896
                           u_count_int,
897
                           count_col_endo,
898
                           symbol_table.exo_det_nbr(),
899
900
                           count_col_det_exo,
                           symbol_table.exo_nbr(),
901
                           count_col_exo,
902
                           0,
903
904
905
906
                           0,
                           exo_det,
                           exo,
                           other_endo
907
                           );
908
  fbeginblock.write(code_file, instruction_number);
909

910
  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
911

912
  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, true, false);
913
914

  FENDEQU_ fendequ;
915
  fendequ.write(code_file, instruction_number);
916
917
918
919
920
921
922

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

923
  vector<vector<pair<pair<int, int>, int >>> derivatives;
924
925
  derivatives.resize(symbol_table.endo_nbr());
  count_u = symbol_table.endo_nbr();
926
  for (const auto & first_derivative : first_derivatives)
927
    {
928
      int deriv_id = first_derivative.first.second;
929
930
      if (getTypeByDerivID(deriv_id) == eEndogenous)
        {
931
932
          expr_t d1 = first_derivative.second;
          unsigned int eq = first_derivative.first.first;
933
934
935
          int symb = getSymbIDByDerivID(deriv_id);
          unsigned int var = symbol_table.getTypeSpecificID(symb);
          int lag = getLagByDerivID(deriv_id);
936
          FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
937
          fnumexpr.write(code_file, instruction_number);
938
939
          if (!derivatives[eq].size())
            derivatives[eq].clear();
940
          derivatives[eq].emplace_back(make_pair(var, lag), count_u);
941
          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
942
943

          FSTPU_ fstpu(count_u);
944
          fstpu.write(code_file, instruction_number);
945
946
947
948
949
950
          count_u++;
        }
    }
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
    {
      FLDR_ fldr(i);
951
      fldr.write(code_file, instruction_number);
952
      if (derivatives[i].size())
953
        {
954
          for (vector<pair<pair<int, int>, int>>::const_iterator it = derivatives[i].begin();
955
               it != derivatives[i].end(); it++)
956
            {
957
958
959
960
961
              FLDU_ fldu(it->second);
              fldu.write(code_file, instruction_number);
              FLDV_ fldv(eEndogenous, it->first.first, it->first.second);
              fldv.write(code_file, instruction_number);
              FBINARY_ fbinary(oTimes);
962
              fbinary.write(code_file, instruction_number);
963
964
965
966
967
              if (it != derivatives[i].begin())
                {
                  FBINARY_ fbinary(oPlus);
                  fbinary.write(code_file, instruction_number);
                }
968
            }
969
970
          FBINARY_ fbinary(oMinus);
          fbinary.write(code_file, instruction_number);
971
972
        }
      FSTPU_ fstpu(i);
973
      fstpu.write(code_file, instruction_number);
974
    }
975
976
977
978
979
980
981
982
983
984
985

  // 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);
986
  prev_instruction_number = instruction_number;
987
988
989
990
991

  // The Jacobian
  prev_var = -1;
  prev_lag = -999999999;
  count_col_endo = 0;
992
  for (map<pair< int, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
993
994
995
996
997
998
999
1000
       it != first_derivatives_reordered_endo.end(); it++)
    {
      unsigned int eq = it->first.second.second;
      int var = it->first.second.first;
      int lag = it->first.first;
      expr_t d1 = it->second;
      FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
      fnumexpr.write(code_file, instruction_number);
1001
      if (prev_var != var || prev_lag != lag)
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
      d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
      FSTPG3_ fstpg3(eq, var, lag, count_col_endo-1);
      fstpg3.write(code_file, instruction_number);
    }
  prev_var = -1;
  prev_lag = -999999999;
1013
  count_col_exo = 0;
1014
  for (map<pair< pair<int, int>, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_exo.begin();
1015
1016
1017
1018
       it != first_derivatives_reordered_exo.end(); it++)
    {
      unsigned int eq = it->first.second.second;
      int var = it->first.second.first;
1019
      int lag = it->first.first.first;
1020
1021
1022
      expr_t d1 = it->second;
      FNUMEXPR_ fnumexpr(FirstExoDerivative, eq, var, lag);
      fnumexpr.write(code_file, instruction_number);
1023
      if (prev_var != var || prev_lag != lag)
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
        {
          prev_var = var;
          prev_lag = lag;
          count_col_exo++;
        }
      d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
      FSTPG3_ fstpg3(eq, var, lag, count_col_exo-1);
      fstpg3.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);

1040
  FENDBLOCK_ fendblock;
1041
  fendblock.write(code_file, instruction_number);
1042
  FEND_ fend;
1043
  fend.write(code_file, instruction_number);
1044
1045
1046
1047
  code_file.close();
}

void
1048
DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin_basename, const map_idx_t &map_idx) const
1049
1050
{
  struct Uff_l
sebastien's avatar
sebastien committed
1051
  {
1052
1053
1054