DynamicModel.cc

/*
 * Copyright (C) 2003-2009 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 <cstdio>
#include <cerrno>
#include <algorithm>
#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,
                           NumericalConstants &num_constants_arg) :
    ModelTree(symbol_table_arg, num_constants_arg),
    max_lag(0), max_lead(0),
    max_endo_lag(0), max_endo_lead(0),
    max_exo_lag(0), max_exo_lead(0),
    max_exo_det_lag(0), max_exo_det_lead(0),
    dynJacobianColsNbr(0),
    global_temporary_terms(true),
    cutoff(1e-15),
    mfs(0)
{
}

VariableNode *
DynamicModel::AddVariable(int symb_id, int lag)
{
  return AddVariableInternal(symb_id, lag);
}

void
DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const
  {
    first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag)));
    if (it != first_derivatives.end())
      (it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
    else
      {
        FLDZ_ fldz;
        fldz.write(code_file);
      }
  }


void
DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, map_idx_type &map_idx) const
{
  map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
  if (it != first_chain_rule_derivatives.end())
    (it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
  else
    {
      FLDZ_ fldz;
      fldz.write(code_file);
    }
}


void
DynamicModel::computeTemporaryTermsOrdered()
{
  map<NodeID, pair<int, int> > first_occurence;
  map<NodeID, int> reference_count;
  BinaryOpNode *eq_node;
  first_derivatives_type::const_iterator it;
  first_chain_rule_derivatives_type::const_iterator it_chr;
  ostringstream tmp_s;
  v_temporary_terms.clear();
  map_idx.clear();

  unsigned int nb_blocks = getNbBlocks();
  v_temporary_terms = vector<vector<temporary_terms_type> >(nb_blocks);
  v_temporary_terms_inuse = vector<temporary_terms_inuse_type> (nb_blocks);
  temporary_terms.clear();

  if(!global_temporary_terms)
    {
      for (unsigned int block = 0; block < nb_blocks; block++)
        {
          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;
          v_temporary_terms[block] = vector<temporary_terms_type>(block_size);
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i<block_nb_recursives && isBlockEquationRenormalized( block, i))
                getBlockEquationRenormalizedNodeID( block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
              else
                {
                  eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
                }
            }
          for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
            {
              NodeID id=it->second.second;
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
            }
          for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
          for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
            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;
        }
    }
  else
    {
      for (unsigned int block = 0; block < nb_blocks; block++)
        {
          // Compute the temporary terms reordered
          unsigned int block_size = getBlockSize(block);
          unsigned int block_nb_mfs = getBlockMfs(block);
          unsigned int block_nb_recursives = block_size - block_nb_mfs;
          v_temporary_terms[block] = vector<temporary_terms_type>(block_size);
          for (unsigned int i = 0; i < block_size; i++)
            {
              if (i<block_nb_recursives && isBlockEquationRenormalized( block, i))
                getBlockEquationRenormalizedNodeID( block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
              else
                {
                  eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
                }
            }
          for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
            {
              NodeID id=it->second.second;
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
            }
          for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
          for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
        }
      for (unsigned int block = 0; block < nb_blocks; block++)
        {
          // 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++)
            {
              if (i<block_nb_recursives && isBlockEquationRenormalized( block, i))
                  getBlockEquationRenormalizedNodeID( block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
              else
                {
                  eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
                  eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
                }
            }
          for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
            {
              NodeID id=it->second.second;
              id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
            }
          for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
            it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
            it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          v_temporary_terms_inuse[block] = temporary_terms_in_use;
        }
      // Add a mapping form node ID to temporary terms order
      int j=0;
      for (temporary_terms_type::const_iterator it = temporary_terms.begin();
           it != temporary_terms.end(); it++)
        map_idx[(*it)->idx]=j++;
    }
}

void
DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
  {
    string tmp_s, sps;
    ostringstream tmp_output, tmp1_output, global_output;
    NodeID lhs=NULL, rhs=NULL;
    BinaryOpNode *eq_node;
    ostringstream Uf[symbol_table.endo_nbr()];
    map<NodeID, int> reference_count;
    temporary_terms_type local_temporary_terms;
    ofstream  output;
    int nze, nze_exo, nze_other_endo;
    vector<int> feedback_variables;
    ExprNodeOutputType local_output_type;

    if(global_temporary_terms)
      {
        local_output_type = oMatlabDynamicModelSparse;
        local_temporary_terms = temporary_terms;
      }
    else
      local_output_type = oMatlabDynamicModelSparseLocalTemporaryTerms;

    //----------------------------------------------------------------------
    //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();
        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;
        unsigned int block_exo_size = exo_block[block].size();
        unsigned int block_exo_det_size = exo_det_block[block].size();
        unsigned int block_other_endo_size = other_endo_block[block].size();
        int block_max_lag=max_leadlag_block[block].first;
        if(global_temporary_terms)
          {
            local_output_type = oMatlabDynamicModelSparse;
            local_temporary_terms = temporary_terms;
          }
        else
          local_output_type = oMatlabDynamicModelSparseLocalTemporaryTerms;

        tmp1_output.str("");
        tmp1_output << dynamic_basename << "_" << block+1 << ".m";
        output.open(tmp1_output.str().c_str(), ios::out | ios::binary);
        output << "%\n";
        output << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare\n";
        output << "%\n";
        output << "% Warning : this file is generated automatically by Dynare\n";
        output << "%           from model file (.mod)\n\n";
        output << "%/\n";
        if (simulation_type ==EVALUATE_BACKWARD || simulation_type ==EVALUATE_FORWARD)
          {
            output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, jacobian_eval, y_kmin, periods)\n";
          }
        else if (simulation_type ==SOLVE_FORWARD_COMPLETE || simulation_type ==SOLVE_BACKWARD_COMPLETE)
          output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, it_, jacobian_eval)\n";
        else if (simulation_type ==SOLVE_BACKWARD_SIMPLE || simulation_type ==SOLVE_FORWARD_SIMPLE)
          output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, it_, jacobian_eval)\n";
        else
          output << "function [residual, y, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, periods, jacobian_eval, y_kmin, y_size)\n";
        BlockType block_type;
        if(simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE ||simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
          block_type = SIMULTAN;
        else if(simulation_type == SOLVE_FORWARD_COMPLETE ||simulation_type == SOLVE_BACKWARD_COMPLETE)
          block_type = SIMULTANS;
        else if((simulation_type == SOLVE_FORWARD_SIMPLE ||simulation_type == SOLVE_BACKWARD_SIMPLE ||
                 simulation_type == EVALUATE_BACKWARD    || simulation_type == EVALUATE_FORWARD)
                && getBlockFirstEquation(block) < prologue)
          block_type = PROLOGUE;
        else if((simulation_type == SOLVE_FORWARD_SIMPLE ||simulation_type == SOLVE_BACKWARD_SIMPLE ||
                 simulation_type == EVALUATE_BACKWARD    || simulation_type == EVALUATE_FORWARD)
                && getBlockFirstEquation(block) >= equations.size() - epilogue)
          block_type = EPILOGUE;
        else
          block_type = SIMULTANS;
        output << "  % ////////////////////////////////////////////////////////////////////////" << endl
        << "  % //" << string("                     Block ").substr(int(log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
        << "          //" << endl
        << "  % //                     Simulation type "
        << BlockSim(simulation_type) << "  //" << endl
        << "  % ////////////////////////////////////////////////////////////////////////" << endl;
        output << "  global options_;" << endl;
        //The Temporary terms
        if (simulation_type ==EVALUATE_BACKWARD || simulation_type ==EVALUATE_FORWARD)
          {
            output << "  if(jacobian_eval)\n";
            output << "    g1 = spalloc(" << block_mfs  << ", " << block_mfs*(1+getBlockMaxLag(block)+getBlockMaxLead(block)) << ", " << nze << ");\n";
            output << "    g1_x=spalloc(" << block_size << ", " << (block_exo_size + block_exo_det_size)*
                      (1+max(exo_det_max_leadlag_block[block].first, exo_max_leadlag_block[block].first)+max(exo_det_max_leadlag_block[block].second, exo_max_leadlag_block[block].second))
                   << ", " << nze_exo << ");\n";
            output << "    g1_o=spalloc(" << block_size << ", " << block_other_endo_size*
                      (1+other_endo_max_leadlag_block[block].first+other_endo_max_leadlag_block[block].second)
                   << ", " << nze_other_endo << ");\n";
            output << "  end;\n";
          }
        else
          {
            output << "  if(jacobian_eval)\n";
            output << "    g1 = spalloc(" << block_size << ", " << block_size*(1+getBlockMaxLag(block)+getBlockMaxLead(block)) << ", " << nze << ");\n";
            output << "    g1_x=spalloc(" << block_size << ", " << (block_exo_size + block_exo_det_size)*
                      (1+max(exo_det_max_leadlag_block[block].first, exo_max_leadlag_block[block].first)+max(exo_det_max_leadlag_block[block].second, exo_max_leadlag_block[block].second))
                   << ", " << nze_exo << ");\n";
            output << "    g1_o=spalloc(" << block_size << ", " << block_other_endo_size*
                      (1+other_endo_max_leadlag_block[block].first+other_endo_max_leadlag_block[block].second)
                   << ", " << nze_other_endo << ");\n";
            output << "  else\n";
            if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
              {
                output << "    g1 = spalloc(" << block_mfs << "*options_.periods, "
                << block_mfs << "*(options_.periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
                << ", " << nze << "*options_.periods);\n";
              }
            else
              {
                output << "    g1 = spalloc(" << block_mfs
                << ", " << block_mfs << ", " << nze << ");\n";
                output << "    g1_tmp_r = spalloc(" << block_recursive
                << ", " << block_size << ", " << nze << ");\n";
                output << "    g1_tmp_b = spalloc(" << block_mfs
                << ", " << block_size << ", " << nze << ");\n";
              }
            output << "  end;\n";
          }

        output << "  g2=0;g3=0;\n";
        if (v_temporary_terms_inuse[block].size())
          {
            tmp_output.str("");
            for (temporary_terms_inuse_type::const_iterator it = v_temporary_terms_inuse[block].begin();
                 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)
          {
            temporary_terms_type tt2;
            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;
                    for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
                         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++)
          {
            temporary_terms_type tt2;
            tt2.clear();
            if (v_temporary_terms[block].size())
              {
                output << "  " << "% //Temporary variables" << endl;
                for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
                     it != v_temporary_terms[block][i].end(); it++)
                  {
                    output << "  " <<  sps;
                    (*it)->writeOutput(output, local_output_type, local_temporary_terms);
                    output << " = ";
                    (*it)->writeOutput(output, local_output_type, tt2);
                    // 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)) ;
            eq_node = (BinaryOpNode*)getBlockEquationNodeID(block,i);
            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("");
                        eq_node = (BinaryOpNode *)getBlockEquationRenormalizedNodeID(block, i);
                        lhs = eq_node->get_arg1();
                        rhs = eq_node->get_arg2();
                        lhs->writeOutput(output, local_output_type, local_temporary_terms);
                        output << " = ";
                        rhs->writeOutput(output, local_output_type, local_temporary_terms);
                      }
                  }
                else
                  {
                    cerr << "Type missmatch for equation " << equation_ID+1  << "\n";
                    exit(EXIT_FAILURE);
                  }
                output << ";\n";
                break;
              case SOLVE_BACKWARD_SIMPLE:
              case SOLVE_FORWARD_SIMPLE:
              case SOLVE_BACKWARD_COMPLETE:
              case SOLVE_FORWARD_COMPLETE:
                if (i<block_recursive)
                  goto evaluation;
                feedback_variables.push_back(variable_ID);
                output << "  % equation " << equation_ID+1 << " variable : " << sModel
                << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
                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
                << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
                Uf[equation_ID] << "    b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
                output << "    residual(" << i+1-block_recursive << ", it_) = (";
                goto end;
              default:
end:
                output << tmp_output.str();
                output << ") - (";
                rhs->writeOutput(output, local_output_type, local_temporary_terms);
                output << ");\n";
#ifdef CONDITION
                if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
                  output << "  condition(" << i+1 << ")=0;\n";
#endif
              }
          }
        // The Jacobian if we have to solve the block
        if (simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE)
          output << "  " << sps << "% Jacobian  " << endl;
        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;
          else
            output << "    % Jacobian  " << endl << "    if jacobian_eval" << endl;
        switch (simulation_type)
          {
          case EVALUATE_BACKWARD:
          case EVALUATE_FORWARD:
            for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
              {
                int lag = it->first.first;
                int eq = it->first.second.first;
                int var = it->first.second.second;
                int eqr = getBlockInitialEquationID(block, eq);
                int varr = getBlockInitialVariableID(block, var);

                NodeID id = it->second;

                output << "      g1(" << eqr+1 << ", " << varr+1+(lag+block_max_lag)*block_size << ") = ";
                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;
              }
            for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
              {
                int lag = it->first.first;
                int eq = it->first.second.first;
                int var = it->first.second.second;
                int eqr = getBlockInitialEquationID(block, eq);
                NodeID id = it->second;

                output << "      g1_o(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
                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_o;\n";
            output << "    end;" << endl;
            output << "  end;" << endl;
            break;
          case SOLVE_BACKWARD_SIMPLE:
          case SOLVE_FORWARD_SIMPLE:
          case SOLVE_BACKWARD_COMPLETE:
          case SOLVE_FORWARD_COMPLETE:
            for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
              {
                int lag = it->first.first;
                unsigned int eq = it->first.second.first;
                unsigned int var = it->first.second.second;
                NodeID id = it->second;

                output << "    g1(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
                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;
              }

            for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
              {
                int lag = it->first.first;
                unsigned int eq = it->first.second.first;
                unsigned int var = it->first.second.second;
                NodeID id = it->second;

                output << "    g1_o(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
                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_o;\n";
            output << "  else" << endl;
            for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
              {
                unsigned int eq = it->first.first;
                unsigned int var = it->first.second;
                unsigned int eqr = getBlockEquationID(block, eq);
                unsigned int varr = getBlockVariableID(block, var);
                NodeID id = it->second.second;
                int lag = it->second.first;
                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;
              }
            output << "  end;\n";
            break;
          case SOLVE_TWO_BOUNDARIES_SIMPLE:
          case SOLVE_TWO_BOUNDARIES_COMPLETE:
            output << "    if ~jacobian_eval" << endl;
            for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
              {
                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;
                NodeID id = it->second.second;
                int lag = it->second.first;
                if(eq>=block_recursive and var>=block_recursive)
                  {
                    if (lag==0)
                      Uf[eqr] << "+g1(" << eq+1-block_recursive
                        << "+Per_J_, " << var+1-block_recursive
                        << "+Per_K_)*y(it_, " << varr+1 << ")";
                    else if (lag==1)
                      Uf[eqr] << "+g1(" << eq+1-block_recursive
                        << "+Per_J_, " << var+1-block_recursive
                        << "+Per_y_)*y(it_+1, " << varr+1 << ")";
                    else if (lag>0)
                      Uf[eqr] << "+g1(" << eq+1-block_recursive
                        << "+Per_J_, " << var+1-block_recursive
                        << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")";
                    else if (lag<0)
                      Uf[eqr] << "+g1(" << eq+1-block_recursive
                        << "+Per_J_, " << var+1-block_recursive
                        << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")";
                    if (lag==0)
                      tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                        << var+1-block_recursive << "+Per_K_) = ";
                    else if (lag==1)
                      tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                        << var+1-block_recursive << "+Per_y_) = ";
                    else if (lag>0)
                      tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                        << var+1-block_recursive << "+y_size*(it_+" << lag-1 << ")) = ";
                    else if (lag<0)
                      tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                        << var+1-block_recursive << "+y_size*(it_" << lag-1 << ")) = ";
                    output << " " << tmp_output.str();
                    id->writeOutput(output, local_output_type, local_temporary_terms);
                    output << ";";
                    output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
                      << "(" << lag << ") " << varr+1
                      << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
                  }

#ifdef CONDITION
                output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
                output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
#endif
              }
            for (unsigned int i = 0; i < block_size; i++)
              {
                if (i>=block_recursive)
                  output << "  " << Uf[getBlockEquationID(block, i)].str() << ";\n";
#ifdef CONDITION
                output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
                output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
#endif
              }
#ifdef CONDITION
            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";
#endif

            output << "    else" << endl;

            for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
              {
                int lag = it->first.first;
                unsigned int eq = it->first.second.first;
                unsigned int var = it->first.second.second;
                NodeID id = it->second;
                output << "      g1(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
                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;
              }
            for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
              {
                int lag = it->first.first;
                unsigned int eq = it->first.second.first;
                unsigned int var = it->first.second.second;
                NodeID id = it->second;

                output << "      g1_o(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
                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_o;\n";
            output << "    end;\n";
            output << "  end;\n";
            break;
          default:
            break;
          }
        output.close();
      }
  }

void
DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) 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;
    NodeID lhs=NULL, rhs=NULL;
    BinaryOpNode *eq_node;
    Uff Uf[symbol_table.endo_nbr()];
    map<NodeID, int> reference_count;
    vector<int> feedback_variables;
    bool file_open=false;

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

    FDIMT_ fdimt(temporary_terms.size());
    fdimt.write(code_file);

    for (unsigned int block = 0; block < getNbBlocks(); block++)
      {
        feedback_variables.clear();
        if (block>0)
          {
            FENDBLOCK_ fendblock;
            fendblock.write(code_file);
          }
        int count_u;
        int u_count_int=0;
        BlockSimulationType simulation_type = getBlockSimulationType(block);
        unsigned int block_size = getBlockSize(block);
        unsigned int block_mfs = getBlockMfs(block);
        unsigned int block_recursive = block_size - block_mfs;
        int block_max_lag=max_leadlag_block[block].first;

        if (simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
            simulation_type==SOLVE_BACKWARD_COMPLETE || simulation_type==SOLVE_FORWARD_COMPLETE)
          {
            Write_Inf_To_Bin_File(file_name, bin_basename, block, u_count_int,file_open,
                                  simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE);
            file_open=true;
          }
        FBEGINBLOCK_ fbeginblock(block_mfs,
                                 simulation_type,
                                 getBlockFirstEquation(block),
                                 block_size,
                                 variable_reordered,
                                 equation_reordered,
                                 blocks_linear[block],
                                 symbol_table.endo_nbr(),
                                 block_max_lag,
                                 block_max_lag,
                                 u_count_int
                                 );
        fbeginblock.write(code_file);

        // The equations
        for (i = 0;i < (int) block_size;i++)
          {
            //The Temporary terms
            temporary_terms_type tt2;
            tt2.clear();
            if (v_temporary_terms[block][i].size())
              {
                for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
                     it != v_temporary_terms[block][i].end(); it++)
                  {
                    (*it)->compile(code_file, false, tt2, map_idx, true, false);
                    FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
                    fstpt.write(code_file);
                    // Insert current node into tt2
                    tt2.insert(*it);
#ifdef DEBUGC
                    cout << "FSTPT " << v << "\n";
                    code_file.write(&FOK, sizeof(FOK));
                    code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
                    ki++;
#endif

                  }
              }
#ifdef DEBUGC
            for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
                     it != v_temporary_terms[block][i].end(); it++)
              {
                map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
                cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
              }
#endif

            int variable_ID, equation_ID;
            EquationType equ_type;


            switch (simulation_type)
              {
evaluation:
              case EVALUATE_BACKWARD:
              case EVALUATE_FORWARD:
                equ_type = getBlockEquationType(block, i);
                if (equ_type == E_EVALUATE)
                  {
                    eq_node = (BinaryOpNode*)getBlockEquationNodeID(block,i);
                    lhs = eq_node->get_arg1();
                    rhs = eq_node->get_arg2();
                    rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
                    lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
                  }
                else if (equ_type == E_EVALUATE_S)
                  {
                    eq_node = (BinaryOpNode*)getBlockEquationRenormalizedNodeID(block,i);
                    lhs = eq_node->get_arg1();
                    rhs = eq_node->get_arg2();
                    rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
                    lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
                  }
                break;
              case SOLVE_BACKWARD_COMPLETE:
              case SOLVE_FORWARD_COMPLETE:
              case SOLVE_TWO_BOUNDARIES_COMPLETE:
              case SOLVE_TWO_BOUNDARIES_SIMPLE:
                if (i< (int) block_recursive)
                  goto evaluation;
                variable_ID = getBlockVariableID(block, i);
                equation_ID = getBlockEquationID(block, i);
                feedback_variables.push_back(variable_ID);
                Uf[equation_ID].Ufl=NULL;
                goto end;
              default:
end:
                eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
                lhs = eq_node->get_arg1();
                rhs = eq_node->get_arg2();
                lhs->compile(code_file, false, temporary_terms, map_idx, true, false);
                rhs->compile(code_file, false, temporary_terms, map_idx, true, false);

                FBINARY_ fbinary(oMinus);
                fbinary.write(code_file);
                FSTPR_ fstpr(i - block_recursive);
                fstpr.write(code_file);
              }
          }
        FENDEQU_ fendequ;
        fendequ.write(code_file);
        // The Jacobian if we have to solve the block
        if    (simulation_type != EVALUATE_BACKWARD
            && simulation_type != EVALUATE_FORWARD)
          {
            switch (simulation_type)
              {
              case SOLVE_BACKWARD_SIMPLE:
              case SOLVE_FORWARD_SIMPLE:
                compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
                  {
                    FSTPG_ fstpg(0);
                    fstpg.write(code_file);
                  }
                break;

              case SOLVE_BACKWARD_COMPLETE:
              case SOLVE_FORWARD_COMPLETE:
              case SOLVE_TWO_BOUNDARIES_COMPLETE:
              case SOLVE_TWO_BOUNDARIES_SIMPLE:
                count_u = feedback_variables.size();
                for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
                  {
                    unsigned int eq = it->first.first;
                    unsigned int var = it->first.second;
                    unsigned int eqr = getBlockEquationID(block, eq);
                    unsigned int varr = getBlockVariableID(block, var);
                    int lag = it->second.first;
                    if(eq>=block_recursive and var>=block_recursive)
                      {
                        if (!Uf[eqr].Ufl)
                          {
                            Uf[eqr].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
                            Uf[eqr].Ufl_First=Uf[eqr].Ufl;
                          }
                        else
                          {
                            Uf[eqr].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
                            Uf[eqr].Ufl=Uf[eqr].Ufl->pNext;
                          }
                        Uf[eqr].Ufl->pNext=NULL;
                        Uf[eqr].Ufl->u=count_u;
                        Uf[eqr].Ufl->var=varr;
                        Uf[eqr].Ufl->lag=lag;
                        compileChainRuleDerivative(code_file, eqr, varr, lag, map_idx);
                        FSTPU_ fstpu(count_u);
                        fstpu.write(code_file);
                        count_u++;
                      }
                  }
                for (i = 0;i < (int) block_size;i++)
                  {
                    if(i>= (int) block_recursive)
                      {
                        FLDR_ fldr(i-block_recursive);
                        fldr.write(code_file);

                        FLDZ_ fldz;
                        fldz.write(code_file);

                        v=getBlockEquationID(block, i);
                        for (Uf[v].Ufl=Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl=Uf[v].Ufl->pNext)
                          {
                            FLDU_ fldu(Uf[v].Ufl->u);
                            fldu.write(code_file);
                            FLDV_ fldv(eEndogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag);
                            fldv.write(code_file);

                            FBINARY_ fbinary(oTimes);
                            fbinary.write(code_file);

                            FCUML_ fcuml;
                            fcuml.write(code_file);
                          }
                        Uf[v].Ufl=Uf[v].Ufl_First;
                        while (Uf[v].Ufl)
                          {
                            Uf[v].Ufl_First=Uf[v].Ufl->pNext;
                            free(Uf[v].Ufl);
                            Uf[v].Ufl=Uf[v].Ufl_First;
                          }
                        FBINARY_ fbinary(oMinus);
                        fbinary.write(code_file);

                        FSTPU_ fstpu(i - block_recursive);
                        fstpu.write(code_file);
                      }
                  }
                break;
              default:
                break;
              }
          }
      }
    FENDBLOCK_ fendblock;
    fendblock.write(code_file);
    FEND_ fend;
    fend.write(code_file);
    code_file.close();
  }

void
DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
  {
    string filename = dynamic_basename + ".m";

    ofstream mDynamicModelFile;
    mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
    if (!mDynamicModelFile.is_open())
      {
        cerr << "Error: Can't open file " << filename << " for writing" << endl;
        exit(EXIT_FAILURE);
      }
    mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, it_)" << endl
    << "%" << endl
    << "% Status : Computes dynamic model for Dynare" << endl
    << "%" << endl
    << "% Warning : this file is generated automatically by Dynare" << endl
    << "%           from model file (.mod)" << endl << endl;

    if (containsSteadyStateOperator())
      mDynamicModelFile << "global oo_;" << endl << endl;

    writeDynamicModel(mDynamicModelFile, false);

    mDynamicModelFile.close();
  }

void
DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
  {
    string filename = dynamic_basename + ".c";
    ofstream mDynamicModelFile;

    mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
    if (!mDynamicModelFile.is_open())
      {
        cerr << "Error: Can't open file " << filename << " for writing" << endl;
        exit(EXIT_FAILURE);
      }
    mDynamicModelFile << "/*" << endl
                      << " * " << filename << " : Computes dynamic model for Dynare" << endl
                      << " *" << endl
                      << " * Warning : this file is generated automatically by Dynare" << endl
                      << " *           from model file (.mod)" << endl
                      << endl
                      << " */" << endl
                      << "#include <math.h>" << endl
                      << "#include \"mex.h\"" << endl
                      << endl
                      << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
                      << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;

    // Writing the function body
    writeDynamicModel(mDynamicModelFile, true);

    // Writing the gateway routine
    mDynamicModelFile << "/* The gateway routine */" << endl
    << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
    << "{" << endl
    << "  double *y, *x, *params;" << endl
    << "  double *residual, *g1, *v2, *v3;" << endl
    << "  int nb_row_x, it_;" << endl
    << endl
    << "  /* Create a pointer to the input matrix y. */" << endl
    << "  y = mxGetPr(prhs[0]);" << endl
    << endl
    << "  /* Create a pointer to the input matrix x. */" << endl
    << "  x = mxGetPr(prhs[1]);" << endl
    << endl
    << "  /* Create a pointer to the input matrix params. */" << endl
    << "  params = mxGetPr(prhs[2]);" << endl
    << endl
    << "  /* Fetch time index */" << endl
    << "  it_ = (int) mxGetScalar(prhs[3]) - 1;" << endl
    << endl
    << "  /* Gets number of rows of matrix x. */" << endl
    << "  nb_row_x = mxGetM(prhs[1]);" << endl
    << endl
    << "  residual = NULL;" << endl
    << "  if (nlhs >= 1)" << endl
    << "  {" << endl
    << "     /* Set the output pointer to the output matrix residual. */" << endl
    << "     plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
    << "     /* Create a C pointer to a copy of the output matrix residual. */" << endl
    << "     residual = mxGetPr(plhs[0]);" << endl
    << "  }" << endl
    << endl
    << "  g1 = NULL;" << endl
    << "  if (nlhs >= 2)" << endl
    << "  {" << endl
    << "     /* Set the output pointer to the output matrix g1. */" << endl

    << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl
    << "     /* Create a C pointer to a copy of the output matrix g1. */" << endl
    << "     g1 = mxGetPr(plhs[1]);" << endl
    << "  }" << endl
    << endl