DynamicModel.cc 162 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2011 Dynare Team
sebastien's avatar
sebastien committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

20
#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>
sebastien's avatar
sebastien committed
27
28
29
30
31
32
33
34
35
36
37
38
#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,
39
40
41
                           NumericalConstants &num_constants_arg,
                           ExternalFunctionsTable &external_functions_table_arg) :
  ModelTree(symbol_table_arg, num_constants_arg, external_functions_table_arg),
42
43
44
45
46
47
48
49
  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),
  dynJacobianColsNbr(0),
  global_temporary_terms(true),
  cutoff(1e-15),
  mfs(0)
sebastien's avatar
sebastien committed
50
51
52
{
}

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

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

void
73
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
74
{
75
  map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
76
  if (it != first_chain_rule_derivatives.end())
77
    (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
78
  else
79
80
    {
      FLDZ_ fldz;
81
      fldz.write(code_file, instruction_number);
82
    }
83
84
}

85
86
87
void
DynamicModel::initializeVariablesAndEquations()
{
88
  for (int j = 0; j < equation_number(); j++)
89
90
91
92
93
94
    {
      equation_reordered.push_back(j);
      variable_reordered.push_back(j);
    }
}

sebastien's avatar
sebastien committed
95
void
96
DynamicModel::computeTemporaryTermsOrdered()
sebastien's avatar
sebastien committed
97
{
98
99
  map<expr_t, pair<int, int> > first_occurence;
  map<expr_t, int> reference_count;
sebastien's avatar
sebastien committed
100
  BinaryOpNode *eq_node;
101
102
  first_derivatives_t::const_iterator it;
  first_chain_rule_derivatives_t::const_iterator it_chr;
sebastien's avatar
sebastien committed
103
  ostringstream tmp_s;
104
105
  v_temporary_terms.clear();
  map_idx.clear();
sebastien's avatar
sebastien committed
106

107
  unsigned int nb_blocks = getNbBlocks();
108
109
  v_temporary_terms = vector<vector<temporary_terms_t> >(nb_blocks);
  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
sebastien's avatar
sebastien committed
110
  temporary_terms.clear();
111

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

207
208
209
210
211
void
DynamicModel::computeTemporaryTermsMapping()
{
  // Add a mapping form node ID to temporary terms order
  int j = 0;
212
  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
213
       it != temporary_terms.end(); it++)
214
215
216
    map_idx[(*it)->idx] = j++;
}

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

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

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

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

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
      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)
        {
          output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, jacobian_eval, y_kmin, periods)\n";
        }
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
        output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, it_, jacobian_eval)\n";
      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
        output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, it_, jacobian_eval)\n";
      else
        output << "function [residual, y, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, periods, jacobian_eval, y_kmin, y_size)\n";
      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;
376
      output << "  global options_ oo_;" << endl;
377
378
379
380
      //The Temporary terms
      if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
        {
          output << "  if(jacobian_eval)\n";
381
382
          output << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
383
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
384
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
385
386
387
388
389
          output << "  end;\n";
        }
      else
        {
          output << "  if(jacobian_eval)\n";
390
391
          output << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
392
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
393
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
394
395
396
397
398
399
400
          output << "  else\n";
          if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
            {
              output << "    g1 = spalloc(" << block_mfs << "*options_.periods, "
                     << block_mfs << "*(options_.periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
                     << ", " << nze << "*options_.periods);\n";
            }
ferhat's avatar
ferhat committed
401
          else
402
403
404
405
406
407
            {
              output << "    g1 = spalloc(" << block_mfs
                     << ", " << block_mfs << ", " << nze << ");\n";
            }
          output << "  end;\n";
        }
408

409
410
411
412
      output << "  g2=0;g3=0;\n";
      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
413
          for (temporary_terms_inuse_t::const_iterator it = v_temporary_terms_inuse[block].begin();
414
415
416
417
418
419
               it != v_temporary_terms_inuse[block].end(); it++)
            tmp_output << " T" << *it;
          output << "  global" << tmp_output.str() << ";\n";
        }
      if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        {
420
          temporary_terms_t tt2;
421
422
423
424
425
426
427
          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;
428
                  for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
                       it != v_temporary_terms[block][i].end(); it++)
                    {
                      output << "  ";
                      (*it)->writeOutput(output, oMatlabDynamicModel, local_temporary_terms);
                      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
      for (unsigned int i = 0; i < block_size; i++)
        {
464
          temporary_terms_t tt2;
465
466
467
468
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
469
              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
470
471
                   it != v_temporary_terms[block][i].end(); it++)
                {
472
473
474
                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);

475
                  output << "  " <<  sps;
476
                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
477
                  output << " = ";
478
                  (*it)->writeOutput(output, local_output_type, tt2, tef_terms);
479
480
481
482
483
484
485
486
487
488
                  // Insert current node into tt2
                  tt2.insert(*it);
                  output << ";" << endl;
                }
            }

          int variable_ID = getBlockVariableID(block, i);
          int equation_ID = getBlockEquationID(block, i);
          EquationType equ_type = getBlockEquationType(block, i);
          string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, variable_ID));
489
          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
          lhs = eq_node->get_arg1();
          rhs = eq_node->get_arg2();
          tmp_output.str("");
          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms);
          switch (simulation_type)
            {
            case EVALUATE_BACKWARD:
            case EVALUATE_FORWARD:
            evaluation:     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 << " = ";
                  rhs->writeOutput(output, local_output_type, local_temporary_terms);
                }
              else if (equ_type == E_EVALUATE_S)
                {
                  output << "%" << tmp_output.str();
                  output << " = ";
                  if (isBlockEquationRenormalized(block, i))
                    {
                      rhs->writeOutput(output, local_output_type, local_temporary_terms);
                      output << "\n    ";
                      tmp_output.str("");
517
                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
                      lhs = eq_node->get_arg1();
                      rhs = eq_node->get_arg2();
                      lhs->writeOutput(output, local_output_type, local_temporary_terms);
                      output << " = ";
                      rhs->writeOutput(output, local_output_type, local_temporary_terms);
                    }
                }
              else
                {
                  cerr << "Type missmatch for equation " << equation_ID+1  << "\n";
                  exit(EXIT_FAILURE);
                }
              output << ";\n";
              break;
            case SOLVE_BACKWARD_SIMPLE:
            case SOLVE_FORWARD_SIMPLE:
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_COMPLETE:
              if (i < block_recursive)
                goto evaluation;
              feedback_variables.push_back(variable_ID);
              output << "  % equation " << equation_ID+1 << " variable : " << sModel
540
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
541
542
543
544
545
546
547
548
              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
549
                     << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
550
551
552
553
554
555
556
557
558
              Uf[equation_ID] << "    b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
              output << "    residual(" << i+1-block_recursive << ", it_) = (";
              goto end;
            default:
            end:
              output << tmp_output.str();
              output << ") - (";
              rhs->writeOutput(output, local_output_type, local_temporary_terms);
              output << ");\n";
sebastien's avatar
sebastien committed
559
#ifdef CONDITION
560
561
              if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
                output << "  condition(" << i+1 << ")=0;\n";
sebastien's avatar
sebastien committed
562
#endif
563
564
565
566
            }
        }
      // The Jacobian if we have to solve the block
      if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
567
        output << "  " << sps << "% Jacobian  " << endl << "    if jacobian_eval" << endl;
568
569
570
571
      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
572
        else
573
          output << "    % Jacobian  " << endl << "    if jacobian_eval" << endl;
574
575
576
577
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
      for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
578
        {
579
580
581
582
          int lag = it->first.first;
          unsigned int var = it->first.second.first;
          unsigned int eq = it->first.second.second;
          if (var != prev_var || lag != prev_lag)
583
            {
584
585
586
587
              prev_var = var;
              prev_lag = lag;
              count_col++;
            }
588

589
          expr_t id = it->second;
590

591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
          output << "      g1(" << eq+1 << ", " << count_col << ") = ";
          id->writeOutput(output, local_output_type, local_temporary_terms);
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
                 << "(" << lag
                 << ") " << var+1
                 << ", equation=" << eq+1 << endl;
        }
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col = 0;
      for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
        {
          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)
608
            {
609
610
611
              prev_var = var;
              prev_lag = lag;
              count_col++;
612
            }
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
          expr_t id = it->second;
          output << "      g1_x(" << eqr+1 << ", " << count_col << ") = ";
          id->writeOutput(output, local_output_type, local_temporary_terms);
          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;
      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++)
        {
          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)
631
            {
632
633
634
              prev_var = var;
              prev_lag = lag;
              count_col++;
635
            }
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
          expr_t id = it->second;
          output << "      g1_xd(" << eqr+1 << ", " << count_col << ") = ";
          id->writeOutput(output, local_output_type, local_temporary_terms);
          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;
      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++)
        {
          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)
654
            {
655
656
657
              prev_var = var;
              prev_lag = lag;
              count_col++;
658
            }
659
660
          expr_t id = it->second;

661
          output << "      g1_o(" << eqr+1 << ", " << /*var+1+(lag+block_max_lag)*block_size*/ count_col << ") = ";
662
663
664
665
666
667
668
669
670
671
672
673
674
675
          id->writeOutput(output, local_output_type, local_temporary_terms);
          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:
676
677
678
679
680
681
682
683
          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;
684
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
685
686
687
688
689
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
690
              expr_t id = it->second.second;
691
              int lag = it->second.first;
692
693
694
695
696
697
698
699
700
701
              if (lag == 0)
                {
                  output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
                  id->writeOutput(output, local_output_type, local_temporary_terms);
                  output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
                         << "(" << lag
                         << ") " << varr+1
                         << ", equation=" << eqr+1 << endl;
                }

702
703
704
705
706
            }
          output << "  end;\n";
          break;
        case SOLVE_TWO_BOUNDARIES_SIMPLE:
        case SOLVE_TWO_BOUNDARIES_COMPLETE:
707
          output << "    else" << endl;
708
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
709
710
711
712
713
714
            {
              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;
715
              expr_t id = it->second.second;
716
              int lag = it->second.first;
717
              if (eq >= block_recursive && var >= block_recursive)
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
                {
                  if (lag == 0)
                    Uf[eqr] << "+g1(" << eq+1-block_recursive
                            << "+Per_J_, " << var+1-block_recursive
                            << "+Per_K_)*y(it_, " << varr+1 << ")";
                  else if (lag == 1)
                    Uf[eqr] << "+g1(" << eq+1-block_recursive
                            << "+Per_J_, " << var+1-block_recursive
                            << "+Per_y_)*y(it_+1, " << varr+1 << ")";
                  else if (lag > 0)
                    Uf[eqr] << "+g1(" << eq+1-block_recursive
                            << "+Per_J_, " << var+1-block_recursive
                            << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")";
                  else if (lag < 0)
                    Uf[eqr] << "+g1(" << eq+1-block_recursive
                            << "+Per_J_, " << var+1-block_recursive
                            << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")";
                  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();
                  id->writeOutput(output, local_output_type, local_temporary_terms);
                  output << ";";
                  output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
                         << "(" << lag << ") " << varr+1
                         << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
                }
754

sebastien's avatar
sebastien committed
755
#ifdef CONDITION
756
757
              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
              output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
sebastien's avatar
sebastien committed
758
#endif
759
760
761
762
763
            }
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i >= block_recursive)
                output << "  " << Uf[getBlockEquationID(block, i)].str() << ";\n";
sebastien's avatar
sebastien committed
764
#ifdef CONDITION
765
766
              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
767
#endif
768
            }
sebastien's avatar
sebastien committed
769
#ifdef CONDITION
770
771
772
773
774
775
776
777
778
779
780
781
782
783
          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
784
#endif
785
786
          output << "    end;" << endl;
          output << "  end;" << endl;
787
788
789
790
          break;
        default:
          break;
        }
791
      output << "end" << endl;
792
      writePowerDeriv(output, false);
793
794
795
      output.close();
    }
}
sebastien's avatar
sebastien committed
796
797

void
798
DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basename, const map_idx_t &map_idx) const
799
{
800

801
802
  ostringstream tmp_output;
  ofstream code_file;
803
  unsigned int instruction_number = 0;
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
  bool file_open = false;
  string main_name = file_name;

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

  int count_u;
  int u_count_int = 0;
  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;

825
  Write_Inf_To_Bin_File(file_name, u_count_int, file_open, simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE, symbol_table.endo_nbr());
826
827
828
829
  file_open = true;

  //Temporary variables declaration
  FDIMT_ fdimt(temporary_terms.size());
830
831
832
833
  fdimt.write(code_file, instruction_number);
  int other_endo_size = 0;

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

835
  for (int i = 0; i < symbol_table.exo_det_nbr(); i++)
836
    exo_det.push_back(i);
837
  for (int i = 0; i < symbol_table.exo_nbr(); i++)
838
    exo.push_back(i);
839
840
841
842
843
844
845
846
847
848
849
850

  map<pair< int, pair<int, int> >, expr_t> first_derivatives_reordered_endo, first_derivatives_reordered_exo;
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
    {
      int deriv_id = it->first.second;
      unsigned int eq = it->first.first;
      int symb = getSymbIDByDerivID(deriv_id);
      unsigned int var = symbol_table.getTypeSpecificID(symb);
      int lag = getLagByDerivID(deriv_id);
      if (getTypeByDerivID(deriv_id) == eEndogenous)
        first_derivatives_reordered_endo[make_pair(lag, make_pair(var, eq))] = it->second;
851
      else if (getTypeByDerivID(deriv_id) == eExogenous || getTypeByDerivID(deriv_id) == eExogenousDet)
852
853
854
855
856
857
858
859
860
861
        first_derivatives_reordered_exo[make_pair(lag, make_pair(var, eq))] = it->second;
    }
  int prev_var = -1;
  int prev_lag = -999999999;
  int count_col_endo = 0;
  for (map<pair< int, pair<int, int> >, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
       it != first_derivatives_reordered_endo.end(); it++)
    {
      int var = it->first.second.first;
      int lag = it->first.first;
862
      if (prev_var != var || prev_lag != lag)
863
864
865
866
867
868
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
    }
869
870
871
872
873
874
875
876
877
  prev_var = -1;
  prev_lag = -999999999;
  int count_col_exo = 0;

  for (map<pair< int, pair<int, int> >, expr_t>::const_iterator it = first_derivatives_reordered_exo.begin();
       it != first_derivatives_reordered_exo.end(); it++)
    {
      int var = it->first.second.first;
      int lag = it->first.first;
878
      if (prev_var != var || prev_lag != lag)
879
880
881
882
883
884
885
        {
          prev_var = var;
          prev_lag = lag;
          count_col_exo++;
        }
    }

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

908
  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
909

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

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

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

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

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

  // 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);
985
  prev_instruction_number = instruction_number;
986
987
988
989
990
991
992
993
994
995
996
997
998
999

  // The Jacobian
  prev_var = -1;
  prev_lag = -999999999;
  count_col_endo = 0;
  for (map<pair< int, pair<int, int> >, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
       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);
1000
      if (prev_var != var || prev_lag != lag)
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
        {
          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;
1012
  count_col_exo = 0;
1013
1014
1015
1016
1017
1018
1019
1020
1021
  for (map<pair< int, pair<int, int> >, expr_t>::const_iterator it = first_derivatives_reordered_exo.begin();
       it != first_derivatives_reordered_exo.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(FirstExoDerivative, eq, var, lag);
      fnumexpr.write(code_file, instruction_number);
1022
      if (prev_var != var || prev_lag != lag)
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
        {
          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);

1039
  FENDBLOCK_ fendblock;
1040
  fendblock.write(code_file, instruction_number);
1041
  FEND_ fend;
1042
  fend.write(code_file, instruction_number);
1043
  writePowerDeriv(code_file, false);
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
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
    int u, var, lag;
    Uff_l *pNext;
  };

  struct Uff
  {
    Uff_l *Ufl, *Ufl_First;
  };

  int i, v;
  string tmp_s;
  ostringstream tmp_output;
  ofstream code_file;
1065
  unsigned int instruction_number = 0;
1066
  expr_t lhs = NULL, rhs = NULL;
1067
1068
  BinaryOpNode *eq_node;
  Uff Uf[symbol_table.endo_nbr()];
1069
  map<expr_t, int> reference_count;
Ferhat Mihoubi's avatar