Commit 6dc87072 authored by sebastien's avatar sebastien
Browse files

trunk preprocessor:

* created a distinct expression tree for the static model (thus giving better sharing of sub-expressions and better computation of temporary terms for the static model)
* for that purpose, created StaticModel and DynamicModel classes (ModelTree still persists, but only contains code shared between StaticModel and DynamicModel)
* removed sparse static file (to be later replaced by new algorithm for steady state computation on large models)


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2592 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 025fb5ae
......@@ -33,22 +33,8 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ M_
global options_
exist_block_structure=isfield(M_,'block_structure');
%disp(['exist_block_structure = ' int2str(exist_block_structure)]);
if exist_block_structure
%block decomposition is inplemented
[x, info]= feval([M_.fname '_static']);
%size(x)
%info
if info<=0
error('solve_one_boundary has failed in block %d\n',-info/10);
else
info = 0;
end;
return;
end;
options_ = set_default_option(options_,'solve_algo',2);
info = 0;
if options_.solve_algo == 0
......
......@@ -49,21 +49,6 @@ enum BlockType
SIMULTAN = 3 //<! Simultaneous time unseparable block
};
/*enum BlockSimulationType
{
UNKNOWN = -1, //!< Unknown simulation type
EVALUATE_FORWARD = 0, //!< Simple evaluation, normalized variable on left-hand side, forward
EVALUATE_BACKWARD = 1, //!< Simple evaluation, normalized variable on left-hand side, backward
SOLVE_FORWARD_SIMPLE = 2, //!< Block of one equation, newton solver needed, forward
SOLVE_BACKWARD_SIMPLE = 3, //!< Block of one equation, newton solver needed, backward
SOLVE_TWO_BOUNDARIES_SIMPLE = 4, //!< Block of one equation, newton solver needed, forward & ackward
SOLVE_FORWARD_COMPLETE = 5, //!< Block of several equations, newton solver needed, forward
SOLVE_BACKWARD_COMPLETE = 6, //!< Block of several equations, newton solver needed, backward
SOLVE_TWO_BOUNDARIES_COMPLETE = 7, //!< Block of several equations, newton solver needed, forward and backwar
EVALUATE_FORWARD_R = 8, //!< Simple evaluation, normalized variable on right-hand side, forward
EVALUATE_BACKWARD_R = 9 //!< Simple evaluation, normalized variable on right-hand side, backward
};
*/
enum BlockSimulationType
{
UNKNOWN, //!< Unknown simulation type
......@@ -131,4 +116,9 @@ enum BinaryOpcode
oDifferent
};
enum TrinaryOpcode
{
oNormcdf
};
#endif
......@@ -38,19 +38,6 @@ SteadyStatement::writeOutput(ostream &output, const string &basename) const
output << "steady;\n";
}
SteadySparseStatement::SteadySparseStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
}
void
SteadySparseStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
//output << basename << "_static;\n";
output << "steady;\n";
}
CheckStatement::CheckStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
......@@ -914,7 +901,7 @@ ModelComparisonStatement::writeOutput(ostream &output, const string &basename) c
output << "model_comparison(ModelNames_,ModelPriors_,oo_,options_,M_.fname);" << endl;
}
PlannerObjectiveStatement::PlannerObjectiveStatement(ModelTree *model_tree_arg) :
PlannerObjectiveStatement::PlannerObjectiveStatement(StaticModel *model_tree_arg) :
model_tree(model_tree_arg)
{
}
......@@ -938,7 +925,7 @@ void
PlannerObjectiveStatement::computingPass()
{
model_tree->computeStaticHessian = true;
model_tree->computingPass(eval_context_type(), false);
model_tree->computingPass(false);
}
void
......
......@@ -25,7 +25,7 @@
#include "SymbolList.hh"
#include "SymbolTable.hh"
#include "Statement.hh"
#include "ModelTree.hh"
#include "StaticModel.hh"
class SteadyStatement : public Statement
{
......@@ -36,15 +36,6 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const;
};
class SteadySparseStatement : public Statement
{
private:
const OptionsList options_list;
public:
SteadySparseStatement(const OptionsList &options_list_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class CheckStatement : public Statement
{
private:
......@@ -409,12 +400,12 @@ public:
class PlannerObjectiveStatement : public Statement
{
private:
ModelTree *model_tree;
StaticModel *model_tree;
public:
//! Constructor
/*! \param model_tree_arg the model tree used to store the objective function.
It is owned by the PlannerObjectiveStatement, and will be deleted by its destructor */
PlannerObjectiveStatement(ModelTree *model_tree_arg);
PlannerObjectiveStatement(StaticModel *model_tree_arg);
virtual ~PlannerObjectiveStatement();
/*! \todo check there are only endogenous variables at the current period in the objective
(no exogenous, no lead/lag) */
......
......@@ -75,8 +75,8 @@ DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
{
// Simplify x+(-y) in x-y
UnaryOpNode *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
if (uarg2 != NULL && uarg2->op_code == oUminus)
return AddMinus(iArg1, uarg2->arg);
if (uarg2 != NULL && uarg2->get_op_code() == oUminus)
return AddMinus(iArg1, uarg2->get_arg());
// To treat commutativity of "+"
// Nodes iArg1 and iArg2 are sorted by index
......@@ -116,8 +116,8 @@ DataTree::AddUMinus(NodeID iArg1)
{
// Simplify -(-x) in x
UnaryOpNode *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
if (uarg != NULL && uarg->op_code == oUminus)
return uarg->arg;
if (uarg != NULL && uarg->get_op_code() == oUminus)
return uarg->get_arg();
return AddUnaryOp(oUminus, iArg1);
}
......
/*
* 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 <cmath>
#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),
cutoff(1e-15),
markowitz(0.7),
computeJacobian(false),
computeJacobianExo(false),
computeHessian(false),
computeThirdDerivatives(false),
block_triangular(symbol_table_arg)
{
}
void
DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const
{
first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(symb_id, lag)));
if (it != first_derivatives.end())
(it->second)->compile(code_file,false, output_type, temporary_terms, map_idx);
else
code_file.write(&FLDZ, sizeof(FLDZ));
}
void
DynamicModel::BuildIncidenceMatrix()
{
set<pair<int, int> > endogenous, exogenous;
for (int eq = 0; eq < (int) equations.size(); eq++)
{
BinaryOpNode *eq_node = equations[eq];
endogenous.clear();
NodeID Id = eq_node->get_arg1();
Id->collectEndogenous(endogenous);
Id = eq_node->get_arg2();
Id->collectEndogenous(endogenous);
for (set<pair<int, int> >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++)
{
block_triangular.incidencematrix.fill_IM(eq, symbol_table.getTypeSpecificID(it_endogenous->first), it_endogenous->second, eEndogenous);
}
exogenous.clear();
Id = eq_node->get_arg1();
Id->collectExogenous(exogenous);
Id = eq_node->get_arg2();
Id->collectExogenous(exogenous);
for (set<pair<int, int> >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++)
{
block_triangular.incidencematrix.fill_IM(eq, symbol_table.getTypeSpecificID(it_exogenous->first), it_exogenous->second, eExogenous);
}
}
}
void
DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
{
map<NodeID, pair<int, int> > first_occurence;
map<NodeID, int> reference_count;
int i, j, m, eq, var, lag;
temporary_terms_type vect;
ostringstream tmp_output;
BinaryOpNode *eq_node;
first_derivatives_type::const_iterator it;
ostringstream tmp_s;
temporary_terms.clear();
map_idx.clear();
for (j = 0;j < ModelBlock->Size;j++)
{
// Compute the temporary terms reordered
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx);
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
//it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++)
{
eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eExogenous, var), lag)));
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
//jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
if (block_triangular.incidencematrix.Model_Max_Lag_Endo - ModelBlock->Block_List[j].Max_Lag +m >=0)
{
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_other_endo;i++)
{
eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
//it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
}
}
for (j = 0;j < ModelBlock->Size;j++)
{
// Compute the temporary terms reordered
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
eq_node->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
//it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++)
{
eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eExogenous, var), lag)));
//it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
//jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
if (block_triangular.incidencematrix.Model_Max_Lag_Endo - ModelBlock->Block_List[j].Max_Lag +m >=0)
{
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_other_endo;i++)
{
eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
//it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
}
}
// Add a mapping form node ID to temporary terms order
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( Model_Block *ModelBlock, const string &dynamic_basename) const
{
int i,j,k,m;
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;
int prev_Simulation_Type=-1, count_derivates=0;
int jacobian_max_endo_col;
ofstream output;
temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
int nze, nze_exo, nze_other_endo;
//----------------------------------------------------------------------
//For each block
for (j = 0;j < ModelBlock->Size;j++)
{
//For a block composed of a single equation determines wether we have to evaluate or to solve the equation
nze = nze_exo = nze_other_endo =0;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
nze+=ModelBlock->Block_List[j].IM_lead_lag[m].size;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead_Exo+ModelBlock->Block_List[j].Max_Lag_Exo;m++)
nze_exo+=ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead_Other_Endo+ModelBlock->Block_List[j].Max_Lag_Other_Endo;m++)
nze_other_endo+=ModelBlock->Block_List[j].IM_lead_lag[m].size_other_endo;
tmp1_output.str("");
tmp1_output << dynamic_basename << "_" << j+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 (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
{
output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, jacobian_eval, y_kmin, periods)\n";
}
else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE)
output << "function [residual, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE)
output << "function [residual, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
else
output << "function [residual, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, periods, jacobian_eval, y_kmin, y_size)\n";
output << " % ////////////////////////////////////////////////////////////////////////" << endl
<< " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type)
<< " //" << endl
<< " % // Simulation type "
<< BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
//The Temporary terms
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
{
output << " if(jacobian_eval)\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size*(1+ModelBlock->Block_List[j].Max_Lag_Endo+ModelBlock->Block_List[j].Max_Lead_Endo) << ", " << nze << ");\n";
output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
output << " end;\n";
}
else
{
output << " if(jacobian_eval)\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size*(1+ModelBlock->Block_List[j].Max_Lag_Endo+ModelBlock->Block_List[j].Max_Lead_Endo) << ", " << nze << ");\n";
output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
output << " else\n";
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size*ModelBlock->Periods << ", " << ModelBlock->Block_List[j].Size*(ModelBlock->Periods+ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead) << ", " << nze*ModelBlock->Periods << ");\n";
else
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
output << " end;\n";
}
output << " g2=0;g3=0;\n";
if(ModelBlock->Block_List[j].Temporary_InUse->size())
{
tmp_output.str("");
for (temporary_terms_inuse_type::const_iterator it = ModelBlock->Block_List[j].Temporary_InUse->begin();
it != ModelBlock->Block_List[j].Temporary_InUse->end(); it++)
tmp_output << " T" << *it;
output << " global" << tmp_output.str() << ";\n";
}
output << " residual=zeros(" << ModelBlock->Block_List[j].Size << ",1);\n";
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
output << " for it_ = y_kmin+1:(y_kmin+periods)\n";
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
output << " b = [];\n";
output << " for it_ = y_kmin+1:(periods+y_kmin)\n";
output << " Per_y_=it_*y_size;\n";
output << " Per_J_=(it_-y_kmin-1)*y_size;\n";
output << " Per_K_=(it_-1)*y_size;\n";
sps=" ";
}
else
if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD ||
ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
sps = " ";
else
sps="";
// The equations
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
temporary_terms_type tt2;
tt2.clear();
if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size())
output << " " << sps << "% //Temporary variables" << endl;
for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->begin();
it != ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->end(); it++)
{
output << " " << sps;
(*it)->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " = ";
(*it)->writeOutput(output, oMatlabDynamicModelSparse, tt2);
// Insert current node into tt2
tt2.insert(*it);
output << ";" << endl;
}
string sModel = symbol_table.getName(ModelBlock->Block_List[j].Variable[i]) ;
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
tmp_output.str("");
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
switch (ModelBlock->Block_List[j].Simulation_Type)
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
output << " ";
output << tmp_output.str();
output << " = ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << ";\n";
break;
case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
output << " ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " = ";
lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << ";\n";
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
output << " " << "residual(" << i+1 << ") = (";
goto end;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << "+Per_J_) = -residual(" << i+1 << ", it_)";
output << " residual(" << i+1 << ", it_) = (";
goto end;
default:
end:
output << tmp_output.str();
output << ") - (";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << ");\n";
#ifdef CONDITION
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
output << " condition(" << i+1 << ")=0;\n";
#endif
}
}
// The Jacobian if we have to solve the block
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
output << " " << sps << "% Jacobian " << endl;
else
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE ||
ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
output << " % Jacobian " << endl << " if jacobian_eval" << endl;
else
output << " % Jacobian " << endl << " if jacobian_eval" << endl;
switch (ModelBlock->Block_List[j].Simulation_Type)
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
count_derivates++;
for (m=0;m<ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
output << " g1(" << eqr+1 << ", " << /*varr+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr*/
varr+1+m*ModelBlock->Block_List[j].Size << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
<< "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
<< ") " << var+