DynamicModel.cc

/*
 * Copyright © 2003-2021 Dynare Team
 *
 * 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/>.
 */

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <numeric>
#include <regex>

#include "DynamicModel.hh"

void
DynamicModel::copyHelper(const DynamicModel &m)
{
  auto f = [this](const ExprNode *e) { return e->clone(*this); };

  for (const auto &it : m.static_only_equations)
    static_only_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));

    auto convert_block_derivative = [f](const map<tuple<int, int, int>, expr_t> &dt)
                                    {
                                      map<tuple<int, int, int>, expr_t> dt2;
                                      for (const auto &it : dt)
                                        dt2[it.first] = f(it.second);
                                      return dt2;
                                    };
  for (const auto &it : m.blocks_derivatives_other_endo)
    blocks_derivatives_other_endo.emplace_back(convert_block_derivative(it));
  for (const auto &it : m.blocks_derivatives_exo)
    blocks_derivatives_exo.emplace_back(convert_block_derivative(it));
  for (const auto &it : m.blocks_derivatives_exo_det)
    blocks_derivatives_exo_det.emplace_back(convert_block_derivative(it));

  for (const auto &[key, expr] : m.pac_expectation_substitution)
    pac_expectation_substitution.emplace(key, f(expr));
}

DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
                           NumericalConstants &num_constants_arg,
                           ExternalFunctionsTable &external_functions_table_arg,
                           TrendComponentModelTable &trend_component_model_table_arg,
                           VarModelTable &var_model_table_arg) :
  ModelTree{symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
  trend_component_model_table{trend_component_model_table_arg},
  var_model_table{var_model_table_arg}
{
}

DynamicModel::DynamicModel(const DynamicModel &m) :
  ModelTree{m},
  trend_component_model_table{m.trend_component_model_table},
  var_model_table{m.var_model_table},
  balanced_growth_test_tol{m.balanced_growth_test_tol},
  static_only_equations_lineno{m.static_only_equations_lineno},
  static_only_equations_equation_tags{m.static_only_equations_equation_tags},
  deriv_id_table{m.deriv_id_table},
  inv_deriv_id_table{m.inv_deriv_id_table},
  dyn_jacobian_cols_table{m.dyn_jacobian_cols_table},
  max_lag{m.max_lag},
  max_lead{m.max_lead},
  max_endo_lag{m.max_endo_lag},
  max_endo_lead{m.max_endo_lead},
  max_exo_lag{m.max_exo_lag},
  max_exo_lead{m.max_exo_lead},
  max_exo_det_lag{m.max_exo_det_lag},
  max_exo_det_lead{m.max_exo_det_lead},
  max_lag_orig{m.max_lag_orig},
  max_lead_orig{m.max_lead_orig},
  max_lag_with_diffs_expanded_orig{m.max_lag_with_diffs_expanded_orig},
  max_endo_lag_orig{m.max_endo_lag_orig},
  max_endo_lead_orig{m.max_endo_lead_orig},
  max_exo_lag_orig{m.max_exo_lag_orig},
  max_exo_lead_orig{m.max_exo_lead_orig},
  max_exo_det_lag_orig{m.max_exo_det_lag_orig},
  max_exo_det_lead_orig{m.max_exo_det_lead_orig},
  xrefs{m.xrefs},
  xref_param{m.xref_param},
  xref_endo{m.xref_endo},
  xref_exo{m.xref_exo},
  xref_exo_det{m.xref_exo_det},
  nonzero_hessian_eqs{m.nonzero_hessian_eqs},
  dynJacobianColsNbr{m.dynJacobianColsNbr},
  variableMapping{m.variableMapping},
  blocks_other_endo{m.blocks_other_endo},
  blocks_exo{m.blocks_exo},
  blocks_exo_det{m.blocks_exo_det},
  blocks_jacob_cols_endo{m.blocks_jacob_cols_endo},
  blocks_jacob_cols_other_endo{m.blocks_jacob_cols_other_endo},
  blocks_jacob_cols_exo{m.blocks_jacob_cols_exo},
  blocks_jacob_cols_exo_det{m.blocks_jacob_cols_exo_det},
  var_expectation_functions_to_write{m.var_expectation_functions_to_write},
  pac_mce_alpha_symb_ids{m.pac_mce_alpha_symb_ids},
  pac_h0_indices{m.pac_h0_indices},
  pac_h1_indices{m.pac_h1_indices},
  pac_mce_z1_symb_ids{m.pac_mce_z1_symb_ids},
  pac_eqtag_and_lag{m.pac_eqtag_and_lag},
  pac_model_info{m.pac_model_info},
  pac_equation_info{m.pac_equation_info}
{
  copyHelper(m);
}

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

  assert(&trend_component_model_table == &m.trend_component_model_table);
  assert(&var_model_table == &m.var_model_table);
  balanced_growth_test_tol = m.balanced_growth_test_tol;

  static_only_equations_lineno = m.static_only_equations_lineno;
  static_only_equations_equation_tags = m.static_only_equations_equation_tags;
  deriv_id_table = m.deriv_id_table;
  inv_deriv_id_table = m.inv_deriv_id_table;
  dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
  max_lag = m.max_lag;
  max_lead = m.max_lead;
  max_endo_lag = m.max_endo_lag;
  max_endo_lead = m.max_endo_lead;
  max_exo_lag = m.max_exo_lag;
  max_exo_lead = m.max_exo_lead;
  max_exo_det_lag = m.max_exo_det_lag;
  max_exo_det_lead = m.max_exo_det_lead;
  max_lag_orig = m.max_lag_orig;
  max_lead_orig = m.max_lead_orig;
  max_lag_with_diffs_expanded_orig = m.max_lag_with_diffs_expanded_orig;
  max_endo_lag_orig = m.max_endo_lag_orig;
  max_endo_lead_orig = m.max_endo_lead_orig;
  max_exo_lag_orig = m.max_exo_lag_orig;
  max_exo_lead_orig = m.max_exo_lead_orig;
  max_exo_det_lag_orig = m.max_exo_det_lag_orig;
  max_exo_det_lead_orig = m.max_exo_det_lead_orig;
  xrefs = m.xrefs;
  xref_param = m.xref_param;
  xref_endo = m.xref_endo;
  xref_exo = m.xref_exo;
  xref_exo_det = m.xref_exo_det;
  nonzero_hessian_eqs = m.nonzero_hessian_eqs;
  dynJacobianColsNbr = m.dynJacobianColsNbr;
  variableMapping = m.variableMapping;
  blocks_derivatives_other_endo.clear();
  blocks_derivatives_exo.clear();
  blocks_derivatives_exo_det.clear();
  blocks_other_endo = m.blocks_other_endo;
  blocks_exo = m.blocks_exo;
  blocks_exo_det = m.blocks_exo_det;
  blocks_jacob_cols_endo = m.blocks_jacob_cols_endo;
  blocks_jacob_cols_other_endo = m.blocks_jacob_cols_other_endo;
  blocks_jacob_cols_exo = m.blocks_jacob_cols_exo;
  blocks_jacob_cols_exo_det = m.blocks_jacob_cols_exo_det;

  var_expectation_functions_to_write = m.var_expectation_functions_to_write;

  pac_mce_alpha_symb_ids = m.pac_mce_alpha_symb_ids;
  pac_h0_indices = m.pac_h0_indices;
  pac_h1_indices = m.pac_h1_indices;
  pac_mce_z1_symb_ids = m.pac_mce_z1_symb_ids;
  pac_eqtag_and_lag = m.pac_eqtag_and_lag;
  pac_expectation_substitution.clear();
  pac_model_info = m.pac_model_info;
  pac_equation_info = m.pac_equation_info;

  copyHelper(m);

  return *this;
}

void
DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs) const
{
  if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
      it != derivatives[1].end())
    it->second->compile(code_file, instruction_number, false, temporary_terms, temporary_terms_idxs, true, false);
  else
    {
      FLDZ_ fldz;
      fldz.write(code_file, instruction_number);
    }
}

void
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs) const
{
  if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
      it != blocks_derivatives[blk].end())
    it->second->compile(code_file, instruction_number, false, temporary_terms, temporary_terms_idxs, true, false);
  else
    {
      FLDZ_ fldz;
      fldz.write(code_file, instruction_number);
    }
}

void
DynamicModel::additionalBlockTemporaryTerms(int blk,
                                            vector<vector<temporary_terms_t>> &blocks_temporary_terms,
                                            map<expr_t, tuple<int, int, int>> &reference_count) const
{
  for (const auto &[ignore, d] : blocks_derivatives_exo[blk])
    d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
  for (const auto &[ignore, d] : blocks_derivatives_exo_det[blk])
    d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
  for (const auto &[ignore, d] : blocks_derivatives_other_endo[blk])
    d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
}

void
DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const
{
  BlockSimulationType simulation_type = blocks[blk].simulation_type;
  int block_size = blocks[blk].size;
  int block_mfs_size = blocks[blk].mfs_size;
  int block_recursive_size = blocks[blk].getRecursiveSize();

  deriv_node_temp_terms_t tef_terms;

  auto write_eq_tt = [&](int eq)
                     {
                       for (auto it : blocks_temporary_terms[blk][eq])
                         {
                           if (dynamic_cast<AbstractExternalFunctionNode *>(it))
                             it->writeExternalFunctionOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);

                           output << "  ";
                           it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms);
                           output << '=';
                           it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
                           temporary_terms.insert(it);
                           output << ';' << endl;
                         }
                     };

  // The equations
  for (int eq = 0; eq < block_size; eq++)
    {
      write_eq_tt(eq);

      EquationType equ_type = getBlockEquationType(blk, eq);
      BinaryOpNode *e = getBlockEquationExpr(blk, eq);
      expr_t lhs = e->arg1, rhs = e->arg2;
      switch (simulation_type)
        {
        case BlockSimulationType::evaluateBackward:
        case BlockSimulationType::evaluateForward:
          evaluation:
          if (equ_type == EquationType::evaluateRenormalized)
            {
              e = getBlockEquationRenormalizedExpr(blk, eq);
              lhs = e->arg1;
              rhs = e->arg2;
            }
          else if (equ_type != EquationType::evaluate)
            {
              cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1  << endl;
              exit(EXIT_FAILURE);
            }
          output << "  ";
          lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
          output << '=';
          rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
          output << ';' << endl;
          break;
        case BlockSimulationType::solveBackwardSimple:
        case BlockSimulationType::solveForwardSimple:
        case BlockSimulationType::solveBackwardComplete:
        case BlockSimulationType::solveForwardComplete:
        case BlockSimulationType::solveTwoBoundariesComplete:
        case BlockSimulationType::solveTwoBoundariesSimple:
          if (eq < block_recursive_size)
            goto evaluation;
          output << "  residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
                 << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
          goto end;
        default:
        end:
          lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
          output << ")-(";
          rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
          output << ");" << endl;
        }
    }

  // The Jacobian if we have to solve the block

  // Write temporary terms for derivatives
  write_eq_tt(blocks[blk].size);

  if (isCOutput(output_type))
    output << "  if (stochastic_mode) {" << endl;
  else
    output << "  if stochastic_mode" << endl;

  ostringstream i_output, j_output, v_output;
  int line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
  for (const auto &[indices, d] : blocks_derivatives[blk])
    {
      auto [eq, var, lag] = indices;
      int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag });
      i_output << "    g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
      j_output << "    g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
      v_output << "    g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
      v_output << ';' << endl;
      line_counter++;
    }
  assert(line_counter == nze_stochastic+ARRAY_SUBSCRIPT_OFFSET(output_type));
  output << i_output.str() << j_output.str() << v_output.str();

  i_output.str("");
  j_output.str("");
  v_output.str("");
  line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
  for (const auto &[indices, d] : blocks_derivatives_exo[blk])
    {
      auto [eq, var, lag] = indices;
      int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag });
      i_output << "    g1_x_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
      j_output << "    g1_x_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
      v_output << "    g1_x_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
      v_output << ';' << endl;
      line_counter++;
    }
  assert(line_counter == nze_exo+ARRAY_SUBSCRIPT_OFFSET(output_type));
  output << i_output.str() << j_output.str() << v_output.str();

  i_output.str("");
  j_output.str("");
  v_output.str("");
  line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
  for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
    {
      auto [eq, var, lag] = indices;
      int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag });
      i_output << "    g1_xd_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
      j_output << "    g1_xd_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
      v_output << "    g1_xd_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
      v_output << ';' << endl;
      line_counter++;
    }
  assert(line_counter == nze_exo_det+ARRAY_SUBSCRIPT_OFFSET(output_type));
  output << i_output.str() << j_output.str() << v_output.str();

  i_output.str("");
  j_output.str("");
  v_output.str("");
  line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
  for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
    {
      auto [eq, var, lag] = indices;
      int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag });
      i_output << "    g1_o_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
      j_output << "    g1_o_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
      v_output << "    g1_o_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
      v_output << ';' << endl;
      line_counter++;
    }
  assert(line_counter == nze_other_endo+ARRAY_SUBSCRIPT_OFFSET(output_type));
  output << i_output.str() << j_output.str() << v_output.str();

  // Deterministic mode
  if (simulation_type != BlockSimulationType::evaluateForward
      && simulation_type != BlockSimulationType::evaluateBackward)
    {
      if (isCOutput(output_type))
        output << "  } else {" << endl;
      else
        output << "  else" << endl;
      i_output.str("");
      j_output.str("");
      v_output.str("");
      line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
      if (simulation_type == BlockSimulationType::solveBackwardSimple
          || simulation_type == BlockSimulationType::solveForwardSimple
          || simulation_type == BlockSimulationType::solveBackwardComplete
          || simulation_type == BlockSimulationType::solveForwardComplete)
        for (const auto &[indices, d] : blocks_derivatives[blk])
          {
            auto [eq, var, lag] = indices;
            if (lag == 0 && eq >= block_recursive_size && var >= block_recursive_size)
              {
                i_output << "    g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
                         << eq+1-block_recursive_size << ';' << endl;
                j_output << "    g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
                         << var+1-block_recursive_size << ';' << endl;
                v_output << "    g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
                d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
                v_output << ';' << endl;
                line_counter++;
              }
          }
      else // solveTwoBoundariesSimple || solveTwoBoundariesComplete
        for (const auto &[indices, d] : blocks_derivatives[blk])
        {
          auto [eq, var, lag] = indices;
          assert(lag >= -1 && lag <= 1);
          if (eq >= block_recursive_size && var >= block_recursive_size)
            {
              i_output << "    g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
                       << eq+1-block_recursive_size << ';' << endl;
              j_output << "    g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
                       << var+1-block_recursive_size+block_mfs_size*(lag+1) << ';' << endl;
              v_output << "    g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
              d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
              v_output << ';' << endl;
              line_counter++;
            }
        }
      assert(line_counter == nze_deterministic+ARRAY_SUBSCRIPT_OFFSET(output_type));
      output << i_output.str() << j_output.str() << v_output.str();
    }
  if (isCOutput(output_type))
    output << "  }" << endl;
  else
    output << "  end" << endl;
}

int
DynamicModel::nzeDeterministicJacobianForBlock(int blk) const
{
  BlockSimulationType simulation_type = blocks[blk].simulation_type;
  int block_recursive_size = blocks[blk].getRecursiveSize();

  int nze_deterministic = 0;
  if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
      || simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
    nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(),
                                 [=](const auto &kv) {
                                   auto [eq, var, lag] = kv.first;
                                   return eq >= block_recursive_size && var >= block_recursive_size;
                                 });
  else if (simulation_type == BlockSimulationType::solveBackwardSimple
           || simulation_type == BlockSimulationType::solveForwardSimple
           || simulation_type == BlockSimulationType::solveBackwardComplete
           || simulation_type == BlockSimulationType::solveForwardComplete)
    nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(),
                                 [=](const auto &kv) {
                                   auto [eq, var, lag] = kv.first;
                                   return lag == 0 && eq >= block_recursive_size && var >= block_recursive_size;
                                 });
  return nze_deterministic;
}

void
DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
{
  temporary_terms_t temporary_terms; // Temp terms written so far

  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
    {
      BlockSimulationType simulation_type = blocks[blk].simulation_type;
      int block_size = blocks[blk].size;
      int block_mfs_size = blocks[blk].mfs_size;

      // Number of nonzero derivatives for the various Jacobians
      int nze_stochastic = blocks_derivatives[blk].size();
      int nze_deterministic = nzeDeterministicJacobianForBlock(blk);
      int nze_other_endo = blocks_derivatives_other_endo[blk].size();
      int nze_exo = blocks_derivatives_exo[blk].size();
      int nze_exo_det = blocks_derivatives_exo_det[blk].size();

      string filename = packageDir(basename + ".block") + "/dynamic_" + to_string(blk+1) + ".m";
      ofstream output;
      output.open(filename, ios::out | ios::binary);
      if (!output.is_open())
        {
          cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
          exit(EXIT_FAILURE);
        }

      output << "%" << endl
             << "% " << filename << " : Computes dynamic version of one block" << endl
             << "%" << endl
             << "% Warning : this file is generated automatically by Dynare" << endl
             << "%           from model file (.mod)" << endl << endl
             << "%" << endl;
      if (simulation_type == BlockSimulationType::evaluateBackward
          || simulation_type == BlockSimulationType::evaluateForward)
        output << "function [y, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode)" << endl;
      else
        output << "function [residual, y, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode)" << endl;

      output << "  % ////////////////////////////////////////////////////////////////////////" << endl
             << "  % //" << string("                     Block ").substr(static_cast<int>(log10(blk + 1))) << blk+1
             << "                                        //" << endl
             << "  % //                     Simulation type "
             << BlockSim(simulation_type) << "  //" << endl
             << "  % ////////////////////////////////////////////////////////////////////////" << endl;

      if (simulation_type != BlockSimulationType::evaluateForward
          && simulation_type != BlockSimulationType::evaluateBackward)
        output << "  residual=zeros(" << block_mfs_size << ",1);" << endl;

      output << "  if stochastic_mode" << endl
             << "    g1_i=zeros(" << nze_stochastic << ",1);" << endl
             << "    g1_j=zeros(" << nze_stochastic << ",1);" << endl
             << "    g1_v=zeros(" << nze_stochastic << ",1);" << endl
             << "    g1_x_i=zeros(" << nze_exo << ",1);" << endl
             << "    g1_x_j=zeros(" << nze_exo << ",1);" << endl
             << "    g1_x_v=zeros(" << nze_exo << ",1);" << endl
             << "    g1_xd_i=zeros(" << nze_exo_det << ",1);" << endl
             << "    g1_xd_j=zeros(" << nze_exo_det << ",1);" << endl
             << "    g1_xd_v=zeros(" << nze_exo_det << ",1);" << endl
             << "    g1_o_i=zeros(" << nze_other_endo << ",1);" << endl
             << "    g1_o_j=zeros(" << nze_other_endo << ",1);" << endl
             << "    g1_o_v=zeros(" << nze_other_endo << ",1);" << endl;
      if (simulation_type != BlockSimulationType::evaluateForward
          && simulation_type != BlockSimulationType::evaluateBackward)
        output << "  else" << endl
               << "    g1_i=zeros(" << nze_deterministic << ",1);" << endl
               << "    g1_j=zeros(" << nze_deterministic << ",1);" << endl
               << "    g1_v=zeros(" << nze_deterministic << ",1);" << endl;
      output << "  end" << endl
             << endl;

      writeDynamicPerBlockHelper(blk, output, ExprNodeOutputType::matlabDynamicModel, temporary_terms,
                                 nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);

      output << endl
             << "  if stochastic_mode" << endl
             << "    g1=sparse(g1_i, g1_j, g1_v, " << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ");" << endl
             << "    varargout{1}=sparse(g1_x_i, g1_x_j, g1_x_v, " << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ");" << endl
             << "    varargout{2}=sparse(g1_xd_i, g1_xd_j, g1_xd_v, " << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
             << "    varargout{3}=sparse(g1_o_i, g1_o_j, g1_o_v, " << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ");" << endl
             << "  else" << endl;
      switch (simulation_type)
        {
        case BlockSimulationType::evaluateForward:
        case BlockSimulationType::evaluateBackward:
          output << "    g1=[];" << endl;
          break;
        case BlockSimulationType::solveBackwardSimple:
        case BlockSimulationType::solveForwardSimple:
        case BlockSimulationType::solveBackwardComplete:
        case BlockSimulationType::solveForwardComplete:
          output << "    g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
                 << ", " << block_mfs_size << ");" << endl;
          break;
        case BlockSimulationType::solveTwoBoundariesSimple:
        case BlockSimulationType::solveTwoBoundariesComplete:
          output << "    g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
                 << ", " << 3*block_mfs_size << ");" << endl;
          break;
        default:
          break;
        }
      output << "  end" << endl
             << "end" << endl;
      output.close();
    }
}

void
DynamicModel::writeDynamicPerBlockCFiles(const string &basename) const
{
  temporary_terms_t temporary_terms; // Temp terms written so far

  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
    {
      BlockSimulationType simulation_type = blocks[blk].simulation_type;
      int block_size = blocks[blk].size;
      int block_mfs_size = blocks[blk].mfs_size;

      // Number of nonzero derivatives for the various Jacobians
      int nze_stochastic = blocks_derivatives[blk].size();
      int nze_deterministic = nzeDeterministicJacobianForBlock(blk);
      int nze_other_endo = blocks_derivatives_other_endo[blk].size();
      int nze_exo = blocks_derivatives_exo[blk].size();
      int nze_exo_det = blocks_derivatives_exo_det[blk].size();

      string filename = basename + "/model/src/dynamic_" + to_string(blk+1) + ".c";
      ofstream output;
      output.open(filename, ios::out | ios::binary);
      if (!output.is_open())
        {
          cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
          exit(EXIT_FAILURE);
        }

      output << "/* Block " << blk+1 << endl
             << "   " << BlockSim(simulation_type) << " */" << endl
             << endl
             << "#include <math.h>" << endl
             << "#include <stdlib.h>" << endl
             << "#include <stdbool.h>" << endl
             << R"(#include "mex.h")" << endl
             << endl;

      // Write function definition if BinaryOpcode::powerDeriv is used
      writePowerDerivHeader(output);

      output << endl;

      if (simulation_type == BlockSimulationType::evaluateBackward
          || simulation_type == BlockSimulationType::evaluateForward)
        output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
      else
        output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
      output << '{' << endl;

      writeDynamicPerBlockHelper(blk, output, ExprNodeOutputType::CDynamicModel, temporary_terms,
                                 nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);

      output << '}' << endl
             << endl;

      ostringstream header;
      if (simulation_type == BlockSimulationType::evaluateBackward
          || simulation_type == BlockSimulationType::evaluateForward)
        header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)";
      else
        header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **residual, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)";
      output << header.str() << endl
             << '{' << endl
             << "  int nb_row_x = mxGetM(x);" << endl;

      if (simulation_type != BlockSimulationType::evaluateForward
          && simulation_type != BlockSimulationType::evaluateBackward)
        output << "  *residual = mxCreateDoubleMatrix(" << block_mfs_size << ",1,mxREAL);" << endl;

      output << "  mxArray *g1_i = NULL, *g1_j = NULL, *g1_v = NULL;" << endl
             << "  mxArray *g1_x_i = NULL, *g1_x_j = NULL, *g1_x_v = NULL;" << endl
             << "  mxArray *g1_xd_i = NULL, *g1_xd_j = NULL, *g1_xd_v = NULL;" << endl
             << "  mxArray *g1_o_i = NULL, *g1_o_j = NULL, *g1_o_v = NULL;" << endl
             << "  if (mxGetScalar(stochastic_mode)) {" << endl
             << "    g1_i=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
             << "    g1_j=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
             << "    g1_v=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
             << "    g1_x_i=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
             << "    g1_x_j=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
             << "    g1_x_v=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
             << "    g1_xd_i=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
             << "    g1_xd_j=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
             << "    g1_xd_v=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
             << "    g1_o_i=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl
             << "    g1_o_j=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl
             << "    g1_o_v=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl;
      if (simulation_type != BlockSimulationType::evaluateForward
          && simulation_type != BlockSimulationType::evaluateBackward)
        output << "  } else {" << endl
               << "    g1_i=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl
               << "    g1_j=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl
               << "    g1_v=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl;
      output << "  }" << endl
             << endl;

      // N.B.: In the following, it_ is decreased by 1, to follow C convention
      if (simulation_type == BlockSimulationType::evaluateBackward
          || simulation_type == BlockSimulationType::evaluateForward)
        output << "  dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_)-1, mxGetScalar(stochastic_mode), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl;
      else
        output << "  dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_)-1, mxGetScalar(stochastic_mode), mxGetPr(*residual), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl;

      output << endl
             << "  if (mxGetScalar(stochastic_mode)) {" << endl
             << "    mxArray *m = mxCreateDoubleScalar(" << block_size << ");" << endl
             << "    mxArray *n = mxCreateDoubleScalar(" << blocks_jacob_cols_endo[blk].size() << ");" << endl
             << "    mxArray *plhs[1];" << endl
             << "    mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
             << R"(    mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
             << "    *g1=plhs[0];" << endl
             << "    mxDestroyArray(g1_i);" << endl
             << "    mxDestroyArray(g1_j);" << endl
             << "    mxDestroyArray(g1_v);" << endl
             << "    mxDestroyArray(n);" << endl
             << "    n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo[blk].size() << ");" << endl
             << "    mxArray *prhs_x[5] = { g1_x_i, g1_x_j, g1_x_v, m, n };" << endl
             << R"(    mexCallMATLAB(1, plhs, 5, prhs_x, "sparse");)" << endl
             << "    *g1_x=plhs[0];" << endl
             << "    mxDestroyArray(g1_x_i);" << endl
             << "    mxDestroyArray(g1_x_j);" << endl
             << "    mxDestroyArray(g1_x_v);" << endl
             << "    mxDestroyArray(n);" << endl
             << "    n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
             << "    mxArray *prhs_xd[5] = { g1_xd_i, g1_xd_j, g1_xd_v, m, n };" << endl
             << R"(    mexCallMATLAB(1, plhs, 5, prhs_xd, "sparse");)" << endl
             << "    *g1_xd=plhs[0];" << endl
             << "    mxDestroyArray(g1_xd_i);" << endl
             << "    mxDestroyArray(g1_xd_j);" << endl
             << "    mxDestroyArray(g1_xd_v);" << endl
             << "    mxDestroyArray(n);" << endl
             << "    n = mxCreateDoubleScalar(" << blocks_jacob_cols_other_endo[blk].size() << ");" << endl
             << "    mxArray *prhs_o[5] = { g1_o_i, g1_o_j, g1_o_v, m, n };" << endl
             << R"(    mexCallMATLAB(1, plhs, 5, prhs_o, "sparse");)" << endl
             << "    *g1_o=plhs[0];" << endl
             << "    mxDestroyArray(g1_o_i);" << endl
             << "    mxDestroyArray(g1_o_j);" << endl
             << "    mxDestroyArray(g1_o_v);" << endl
             << "    mxDestroyArray(n);" << endl
             << "    mxDestroyArray(m);" << endl
             << "  } else {" << endl;
      switch (simulation_type)
        {
        case BlockSimulationType::evaluateForward:
        case BlockSimulationType::evaluateBackward:
          output << "    *g1=mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
          break;
        case BlockSimulationType::solveBackwardSimple:
        case BlockSimulationType::solveForwardSimple:
        case BlockSimulationType::solveBackwardComplete:
        case BlockSimulationType::solveForwardComplete:
          output << "    mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
                 << "    mxArray *n = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
                 << "    mxArray *plhs[1];" << endl
                 << "    mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
                 << R"(    mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
                 << "    *g1=plhs[0];" << endl
                 << "    mxDestroyArray(g1_i);" << endl
                 << "    mxDestroyArray(g1_j);" << endl
                 << "    mxDestroyArray(g1_v);" << endl
                 << "    mxDestroyArray(n);" << endl
                 << "    mxDestroyArray(m);" << endl;
          break;
        case BlockSimulationType::solveTwoBoundariesSimple:
        case BlockSimulationType::solveTwoBoundariesComplete:
          output << "    mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
                 << "    mxArray *n = mxCreateDoubleScalar(" << 3*block_mfs_size << ");" << endl
                 << "    mxArray *plhs[1];" << endl
                 << "    mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
                 << R"(    mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
                 << "    *g1=plhs[0];" << endl
                 << "    mxDestroyArray(g1_i);" << endl
                 << "    mxDestroyArray(g1_j);" << endl
                 << "    mxDestroyArray(g1_v);" << endl
                 << "    mxDestroyArray(n);" << endl
                 << "    mxDestroyArray(m);" << endl;
          break;
        default:
          break;
        }
      output << "    *g1_x=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
             << "    *g1_xd=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
             << "    *g1_o=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
             << "  }" << endl
             << "}" << endl;
      output.close();

      filename = basename + "/model/src/dynamic_" + to_string(blk+1) + ".h";
      ofstream header_output;
      header_output.open(filename, ios::out | ios::binary);
      if (!header_output.is_open())
        {
          cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
          exit(EXIT_FAILURE);
        }
      header_output << header.str() << ';' << endl;
      header_output.close();
    }
}

void
DynamicModel::writeDynamicBytecode(const string &basename) const
{
  ostringstream tmp_output;
  ofstream code_file;
  unsigned int instruction_number = 0;
  bool file_open = false;

  string main_name = basename + "/model/bytecode/dynamic.cod";
  code_file.open(main_name, ios::out | ios::binary | ios::ate);
  if (!code_file.is_open())
    {
      cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
      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 = BlockSimulationType::solveTwoBoundariesComplete;
  else if ((max_endo_lag >= 0) && (max_endo_lead == 0))
    simulation_type = BlockSimulationType::solveForwardComplete;
  else
    simulation_type = BlockSimulationType::solveBackwardComplete;

  writeBytecodeBinFile(basename + "/model/bytecode/dynamic.bin", u_count_int, file_open, simulation_type == BlockSimulationType::solveTwoBoundariesComplete);
  file_open = true;

  //Temporary variables declaration
  FDIMT_ fdimt(temporary_terms_idxs.size());
  fdimt.write(code_file, instruction_number);

  vector<int> exo, exo_det, other_endo;

  for (int i = 0; i < symbol_table.exo_det_nbr(); i++)
    exo_det.push_back(i);
  for (int i = 0; i < symbol_table.exo_nbr(); i++)
    exo.push_back(i);

  map<tuple<int, int, int>, expr_t> first_derivatives_reordered_endo;
  map<tuple<int, SymbolType, int, int>, expr_t> first_derivatives_reordered_exo;
  for (const auto & [indices, d1] : derivatives[1])
    {
      int deriv_id = indices[1];
      int eq = indices[0];
      int symb = getSymbIDByDerivID(deriv_id);
      int var = symbol_table.getTypeSpecificID(symb);
      int lag = getLagByDerivID(deriv_id);
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
        first_derivatives_reordered_endo[{ lag, var, eq }] = d1;
      else if (getTypeByDerivID(deriv_id) == SymbolType::exogenous || getTypeByDerivID(deriv_id) == SymbolType::exogenousDet)
        first_derivatives_reordered_exo[{ lag, getTypeByDerivID(deriv_id), var, eq }] = d1;
    }
  int prev_var = -1;
  int prev_lag = -999999999;
  int count_col_endo = 0;
  for (const auto &it : first_derivatives_reordered_endo)
    {
      int var, lag;
      tie(lag, var, ignore) = it.first;
      if (prev_var != var || prev_lag != lag)
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
    }
  prev_var = -1;
  prev_lag = -999999999;
  SymbolType prev_type{SymbolType::unusedEndogenous}; // Any non-exogenous type would do here
  int count_col_exo = 0;
  int count_col_det_exo = 0;

  for (const auto &it : first_derivatives_reordered_exo)
    {
      int var, lag;
      SymbolType type;
      tie(lag, type, var, ignore) = it.first;
      if (prev_var != var || prev_lag != lag || prev_type != type)
        {
          prev_var = var;
          prev_lag = lag;
          prev_type = type;
          if (type == SymbolType::exogenous)
            count_col_exo++;
          else if (type == SymbolType::exogenousDet)
            count_col_det_exo++;
        }
    }

  FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
                           simulation_type,
                           0,
                           symbol_table.endo_nbr(),
                           endo_idx_block2orig,
                           eq_idx_block2orig,
                           false,
                           symbol_table.endo_nbr(),
                           max_endo_lag,
                           max_endo_lead,
                           u_count_int,
                           count_col_endo,
                           symbol_table.exo_det_nbr(),
                           count_col_det_exo,
                           symbol_table.exo_nbr(),
                           count_col_exo,
                           0,
                           0,
                           exo_det,
                           exo,
                           other_endo);
  fbeginblock.write(code_file, instruction_number);

  temporary_terms_t temporary_terms_union;
  compileTemporaryTerms(code_file, instruction_number, true, false, temporary_terms_union, temporary_terms_idxs);

  compileModelEquations(code_file, instruction_number, true, false, temporary_terms_union, temporary_terms_idxs);

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

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

  vector<vector<tuple<int, int, int>>> my_derivatives(symbol_table.endo_nbr());;
  count_u = symbol_table.endo_nbr();
  for (const auto & [indices, d1] : derivatives[1])
    {
      int deriv_id = indices[1];
      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
        {
          int eq = indices[0];
          int symb = getSymbIDByDerivID(deriv_id);
          int var = symbol_table.getTypeSpecificID(symb);
          int lag = getLagByDerivID(deriv_id);
          FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eq, var, lag);
          fnumexpr.write(code_file, instruction_number);
          if (!my_derivatives[eq].size())
            my_derivatives[eq].clear();
          my_derivatives[eq].emplace_back(var, lag, count_u);
          d1->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, true, false);

          FSTPU_ fstpu(count_u);
          fstpu.write(code_file, instruction_number);
          count_u++;
        }
    }
  for (int i = 0; i < symbol_table.endo_nbr(); i++)
    {
      FLDR_ fldr(i);
      fldr.write(code_file, instruction_number);
      if (my_derivatives[i].size())
        {
          for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); ++it)
            {
              FLDU_ fldu(get<2>(*it));
              fldu.write(code_file, instruction_number);
              FLDV_ fldv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(get<0>(*it)), get<1>(*it)};
              fldv.write(code_file, instruction_number);
              FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
              fbinary.write(code_file, instruction_number);
              if (it != my_derivatives[i].begin())
                {
                  FBINARY_ fbinary{static_cast<int>(BinaryOpcode::plus)};
                  fbinary.write(code_file, instruction_number);
                }
            }
          FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
          fbinary.write(code_file, instruction_number);
        }
      FSTPU_ fstpu(i);
      fstpu.write(code_file, instruction_number);
    }

  // 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);
  prev_instruction_number = instruction_number;

  // The Jacobian
  prev_var = -1;
  prev_lag = -999999999;
  count_col_endo = 0;
  for (const auto &it : first_derivatives_reordered_endo)
    {
      auto [lag, var, eq] = it.first;
      expr_t d1 = it.second;
      FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eq, var, lag);
      fnumexpr.write(code_file, instruction_number);
      if (prev_var != var || prev_lag != lag)
        {
          prev_var = var;
          prev_lag = lag;
          count_col_endo++;
        }
      d1->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, true, false);
      FSTPG3_ fstpg3(eq, var, lag, count_col_endo-1);
      fstpg3.write(code_file, instruction_number);
    }
  prev_var = -1;
  prev_lag = -999999999;
  count_col_exo = 0;
  for (const auto &it : first_derivatives_reordered_exo)
    {
      auto [lag, ignore, var, eq] = it.first;
      expr_t d1 = it.second;
      FNUMEXPR_ fnumexpr(ExpressionType::FirstExoDerivative, eq, var, lag);
      fnumexpr.write(code_file, instruction_number);
      if (prev_var != var || prev_lag != lag)
        {
          prev_var = var;
          prev_lag = lag;
          count_col_exo++;
        }
      d1->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, 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);

  FENDBLOCK_ fendblock;
  fendblock.write(code_file, instruction_number);
  FEND_ fend;
  fend.write(code_file, instruction_number);
  code_file.close();
}

void
DynamicModel::writeDynamicBlockBytecode(const string &basename) const
{
  struct Uff_l
  {
    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;
  unsigned int instruction_number = 0;
  expr_t lhs = nullptr, rhs = nullptr;
  BinaryOpNode *eq_node;
  Uff Uf[symbol_table.endo_nbr()];
  map<expr_t, int> reference_count;
  vector<int> feedback_variables;
  bool file_open = false;
  string main_name = basename + "/model/bytecode/dynamic.cod";
  code_file.open(main_name, ios::out | ios::binary | ios::ate);
  if (!code_file.is_open())
    {
      cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
      exit(EXIT_FAILURE);
    }
  //Temporary variables declaration

  FDIMT_ fdimt(blocks_temporary_terms_idxs.size());
  fdimt.write(code_file, instruction_number);

  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
    {
      feedback_variables.clear();
      if (block > 0)
        {
          FENDBLOCK_ fendblock;
          fendblock.write(code_file, instruction_number);
        }
      int count_u;
      int u_count_int = 0;
      BlockSimulationType simulation_type = blocks[block].simulation_type;
      int block_size = blocks[block].size;
      int block_mfs = blocks[block].mfs_size;
      int block_recursive = blocks[block].getRecursiveSize();
      int block_max_lag = blocks[block].max_lag;
      int block_max_lead = blocks[block].max_lead;

      if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
          || simulation_type == BlockSimulationType::solveTwoBoundariesComplete
          || simulation_type == BlockSimulationType::solveBackwardComplete
          || simulation_type == BlockSimulationType::solveForwardComplete)
        {
          writeBlockBytecodeBinFile(basename, block, u_count_int, file_open,
                                    simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple);
          file_open = true;
        }

      FBEGINBLOCK_ fbeginblock(block_mfs,
                               simulation_type,
                               blocks[block].first_equation,
                               block_size,
                               endo_idx_block2orig,
                               eq_idx_block2orig,
                               blocks[block].linear,
                               symbol_table.endo_nbr(),
                               block_max_lag,
                               block_max_lead,
                               u_count_int,
                               blocks_jacob_cols_endo[block].size(),
                               blocks_exo_det[block].size(),
                               blocks_jacob_cols_exo_det[block].size(),
                               blocks_exo[block].size(),
                               blocks_jacob_cols_exo[block].size(),
                               blocks_other_endo[block].size(),
                               blocks_jacob_cols_other_endo[block].size(),
                               vector<int>(blocks_exo_det[block].begin(), blocks_exo_det[block].end()),
                               vector<int>(blocks_exo[block].begin(), blocks_exo[block].end()),
                               vector<int>(blocks_other_endo[block].begin(), blocks_other_endo[block].end()));
      fbeginblock.write(code_file, instruction_number);

      temporary_terms_t temporary_terms_union;

      //The Temporary terms
      deriv_node_temp_terms_t tef_terms;

      auto write_eq_tt = [&](int eq)
                         {
                           for (auto it : blocks_temporary_terms[block][eq])
                             {
                               if (dynamic_cast<AbstractExternalFunctionNode *>(it))
                                 it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);

                               FNUMEXPR_ fnumexpr(ExpressionType::TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
                               fnumexpr.write(code_file, instruction_number);
                               it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
                               FSTPT_ fstpt(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
                               fstpt.write(code_file, instruction_number);
                               temporary_terms_union.insert(it);
#ifdef DEBUGC
                               cout << "FSTPT " << v << endl;
                               instruction_number++;
                               code_file.write(&FOK, sizeof(FOK));
                               code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
                               ki++;
#endif
                             }
                         };

      // The equations
      for (i = 0; i < block_size; i++)
        {
          write_eq_tt(i);

          int variable_ID, equation_ID;
          EquationType equ_type;

          switch (simulation_type)
            {
            evaluation:
            case BlockSimulationType::evaluateBackward:
            case BlockSimulationType::evaluateForward:
              equ_type = getBlockEquationType(block, i);
              {
                FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
                fnumexpr.write(code_file, instruction_number);
              }
              if (equ_type == EquationType::evaluate)
                {
                  eq_node = getBlockEquationExpr(block, i);
                  lhs = eq_node->arg1;
                  rhs = eq_node->arg2;
                  rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false);
                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, true, false);
                }
              else if (equ_type == EquationType::evaluateRenormalized)
                {
                  eq_node = getBlockEquationRenormalizedExpr(block, i);
                  lhs = eq_node->arg1;
                  rhs = eq_node->arg2;
                  rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false);
                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, true, false);
                }
              break;
            case BlockSimulationType::solveBackwardComplete:
            case BlockSimulationType::solveForwardComplete:
            case BlockSimulationType::solveTwoBoundariesComplete:
            case BlockSimulationType::solveTwoBoundariesSimple:
              if (i < block_recursive)
                goto evaluation;
              variable_ID = getBlockVariableID(block, i);
              equation_ID = getBlockEquationID(block, i);
              feedback_variables.push_back(variable_ID);
              Uf[equation_ID].Ufl = nullptr;
              goto end;
            default:
            end:
              FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
              fnumexpr.write(code_file, instruction_number);
              eq_node = getBlockEquationExpr(block, i);
              lhs = eq_node->arg1;
              rhs = eq_node->arg2;
              lhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false);
              rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false);

              FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
              fbinary.write(code_file, instruction_number);
              FSTPR_ fstpr(i - block_recursive);
              fstpr.write(code_file, instruction_number);
            }
        }
      FENDEQU_ fendequ;
      fendequ.write(code_file, instruction_number);

      // 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;
      // The Jacobian if we have to solve the block determinsitic block
      if (simulation_type != BlockSimulationType::evaluateBackward
          && simulation_type != BlockSimulationType::evaluateForward)
        {
          // Write temporary terms for derivatives
          write_eq_tt(blocks[block].size);

          switch (simulation_type)
            {
            case BlockSimulationType::solveBackwardSimple:
            case BlockSimulationType::solveForwardSimple:
              {
                FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
                fnumexpr.write(code_file, instruction_number);