DynamicModel.cc 239 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
2
 * Copyright (C) 2003-2017 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 47
  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),
48
  global_temporary_terms(true)
sebastien's avatar
sebastien committed
49 50 51
{
}

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

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

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

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

96
  unsigned int nb_blocks = getNbBlocks();
97 98
  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
99
  temporary_terms.clear();
100

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

200 201 202 203 204
void
DynamicModel::computeTemporaryTermsMapping()
{
  // Add a mapping form node ID to temporary terms order
  int j = 0;
205
  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
206
       it != temporary_terms.end(); it++)
207 208 209
    map_idx[(*it)->idx] = j++;
}

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

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

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

252 253 254 255
      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++)
256
        tmp_block_endo_derivative[make_pair(it->second.first, make_pair(it->first.second, it->first.first))] = it->second.second;
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
      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;
          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++)
273
        tmp_block_exo_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first))] = it->second;
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
      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;
          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++)
290
        tmp_block_exo_det_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first))] = it->second;
291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
      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;
          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++)
307
        tmp_block_other_endo_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first))] = it->second;
308 309 310
      prev_var = 999999999;
      prev_lag = -9999999;
      count_col_other_endo = 0;
311
      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++)
312 313 314 315 316 317 318 319 320 321 322
        {
          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++;
            }
        }

323 324 325 326 327 328 329 330 331 332 333
      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)
        {
334
          output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)\n";
335 336
        }
      else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
337
        output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
338
      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
339
        output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
340
      else
341
        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";
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
      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";
367 368
          output << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
369
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
370
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
371 372 373 374 375
          output << "  end;\n";
        }
      else
        {
          output << "  if(jacobian_eval)\n";
376 377
          output << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");\n";
          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
378
          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
379
          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
380 381 382
          output << "  else\n";
          if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
            {
383 384 385
              output << "    g1 = spalloc(" << block_mfs << "*Periods, "
                     << block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
                     << ", " << nze << "*Periods);\n";
386
            }
ferhat's avatar
ferhat committed
387
          else
388 389 390 391 392 393
            {
              output << "    g1 = spalloc(" << block_mfs
                     << ", " << block_mfs << ", " << nze << ");\n";
            }
          output << "  end;\n";
        }
394

395 396 397 398
      output << "  g2=0;g3=0;\n";
      if (v_temporary_terms_inuse[block].size())
        {
          tmp_output.str("");
399
          for (temporary_terms_inuse_t::const_iterator it = v_temporary_terms_inuse[block].begin();
400 401 402 403 404 405
               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)
        {
406
          temporary_terms_t tt2;
407 408 409 410 411 412 413
          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;
414
                  for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
415 416 417 418 419 420 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
                       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++)
        {
450
          temporary_terms_t tt2;
451 452 453 454
          tt2.clear();
          if (v_temporary_terms[block].size())
            {
              output << "  " << "% //Temporary variables" << endl;
455
              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
456 457
                   it != v_temporary_terms[block][i].end(); it++)
                {
458
                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
459 460
                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);

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

579
          expr_t id = it->second;
580

581 582
          output << "      g1(" << eq+1 << ", " << count_col << ") = ";
          id->writeOutput(output, local_output_type, local_temporary_terms);
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
583
          output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
584
                 << "(" << lag
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
585 586
                 << ") " << varr+1 << ", " << var+1
                 << ", equation=" << eqr+1 << ", " << eq+1 << endl;
587 588 589 590 591 592 593 594 595 596 597
        }
      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)
598
            {
599 600 601
              prev_var = var;
              prev_lag = lag;
              count_col++;
602
            }
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620
          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)
621
            {
622 623 624
              prev_var = var;
              prev_lag = lag;
              count_col++;
625
            }
626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
          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)
644
            {
645 646 647
              prev_var = var;
              prev_lag = lag;
              count_col++;
648
            }
649 650
          expr_t id = it->second;

651
          output << "      g1_o(" << eqr+1 << ", " << /*var+1+(lag+block_max_lag)*block_size*/ count_col << ") = ";
652 653 654 655 656 657 658 659 660 661 662 663 664 665
          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:
666 667 668 669 670 671 672 673
          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;
674
          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
675 676 677 678 679
            {
              unsigned int eq = it->first.first;
              unsigned int var = it->first.second;
              unsigned int eqr = getBlockEquationID(block, eq);
              unsigned int varr = getBlockVariableID(block, var);
680
              expr_t id = it->second.second;
681
              int lag = it->second.first;
682 683 684 685 686 687 688 689 690 691
              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;
                }

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

728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746
                  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;
                }
747

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

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

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

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

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

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

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

831 832
  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;
833 834 835 836 837 838 839 840 841 842
  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;
843
      else if (getTypeByDerivID(deriv_id) == eExogenous || getTypeByDerivID(deriv_id) == eExogenousDet)
844
        first_derivatives_reordered_exo[make_pair(make_pair(lag, getTypeByDerivID(deriv_id)), make_pair(var, eq))] = it->second;
845 846 847 848 849 850 851 852 853
    }
  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;
854
      if (prev_var != var || prev_lag != lag)
855 856 857 858 859 860
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
    }
861 862
  prev_var = -1;
  prev_lag = -999999999;
863
  int prev_type = -1;
864
  int count_col_exo = 0;
865
  int count_col_det_exo = 0;
866

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

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

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

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

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

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

922 923 924
  vector<vector<pair<pair<int, int>, int > > > derivatives;
  derivatives.resize(symbol_table.endo_nbr());
  count_u = symbol_table.endo_nbr();
925
  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
926 927 928 929 930
       it != first_derivatives.end(); it++)
    {
      int deriv_id = it->first.second;
      if (getTypeByDerivID(deriv_id) == eEndogenous)
        {
931
          expr_t d1 = it->second;
932 933 934 935
          unsigned int eq = it->first.first;
          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 940
          if (!derivatives[eq].size())
            derivatives[eq].clear();
          derivatives[eq].push_back(make_pair(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 955
          for (vector<pair<pair<int, int>, int> >::const_iterator it = derivatives[i].begin();
               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 992 993 994 995 996 997 998 999 1000

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