diff --git a/CodeInterpreter.hh b/CodeInterpreter.hh index 86d4218c1b15e9db798ada2365db9b046bfdd2be..aa043499a3fc9b5fae9351f652ad44a07a325e43 100644 --- a/CodeInterpreter.hh +++ b/CodeInterpreter.hh @@ -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 diff --git a/ComputingTasks.cc b/ComputingTasks.cc index b2b42662c6596c24cd1079c2ae398ad6fc15c6f3..0aa27b6dca28f4cf16690f0cbc153fecd2430adc 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -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 diff --git a/ComputingTasks.hh b/ComputingTasks.hh index 0f268d2b0f1b8ff9ce67a6c3a18e0fb2688f01ce..6cd9581f5fb17df9a9e9c0d61c823bfd3ff4b019 100644 --- a/ComputingTasks.hh +++ b/ComputingTasks.hh @@ -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) */ diff --git a/DataTree.cc b/DataTree.cc index 69adc2e68536fa9421da4c22d399c7c6ada8b25e..0c8d9877118488d8117ae0d1ceb40d7dd67f61df 100644 --- a/DataTree.cc +++ b/DataTree.cc @@ -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); } diff --git a/DynamicModel.cc b/DynamicModel.cc new file mode 100644 index 0000000000000000000000000000000000000000..f41066f1c650a80a3f7a8de3db9e558a2c791ddf --- /dev/null +++ b/DynamicModel.cc @@ -0,0 +1,2249 @@ +/* + * 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+1 + << ", equation=" << eq+1 << endl; + } + } + //jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr; + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + output << " g1_x(" << eqr+1 << ", " + << varr+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr() << ") = "; + writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=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++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i]; + output << " g1_o(" << eqr+1 << ", " + << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = "; + writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << 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: + 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 + << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + output << " g1_x(" << eqr+1 << ", " << varr+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = "; + writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=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++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i]; + output << " g1_o(" << eqr+1 << ", " + << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = "; + writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + } + output << " varargout{1}=g1_x;\n"; + output << " varargout{2}=g1_o;\n"; + output << " else" << endl; + + 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 << ") = "; + writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + output << " end;\n"; + break; + case SOLVE_TWO_BOUNDARIES_SIMPLE: + case SOLVE_TWO_BOUNDARIES_COMPLETE: + output << " if ~jacobian_eval" << endl; + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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]; + if (k==0) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")"; + else if (k==1) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_)*y(it_+1, " << var+1 << ")"; + else if (k>0) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")"; + else if (k<0) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; + if (k==0) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; + else if (k==1) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; + else if (k>0) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; + else if (k<0) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; + writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << var+1 + << ", equation=" << 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 (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + output << " " << Uf[ModelBlock->Block_List[j].Equation[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[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; + int eqr=ModelBlock->Block_List[j].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[j].Size;i++) + output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; +#endif + + output << " else" << endl; + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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-ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lag_Endo)*ModelBlock->Block_List[j].Size << ") = "; + writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size; + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1_x(" << eqr+1 << ", " + << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = "; + writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable (exogenous)=" << symbol_table.getName(var) + << "(" << k << ") " << var+1 << " " << varr+1 + << ", equation=" << eq+1 << endl; + } + } + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=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++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i]; + output << " g1_o(" << eqr+1 << ", " + << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = "; + writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(var) + << "(" << k << ") " << 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; + } + prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; + output.close(); + } +} + +void +DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const +{ + struct Uff_l + { + int u, var, lag; + Uff_l *pNext; + }; + + struct Uff + { + Uff_l *Ufl, *Ufl_First; + int eqr; + }; + + int i,j,k,m, v, ModelBlock_Aggregated_Count, k0, k1; + string tmp_s; + ostringstream tmp_output; + ofstream code_file; + NodeID lhs=NULL, rhs=NULL; + BinaryOpNode *eq_node; + bool lhs_rhs_done; + Uff Uf[symbol_table.endo_nbr()]; + map<NodeID, int> reference_count; + map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number; + int prev_Simulation_Type=-1; + //SymbolicGaussElimination SGE; + bool file_open=false; + temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); + //---------------------------------------------------------------------- + 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 + code_file.write(&FDIMT, sizeof(FDIMT)); + k=temporary_terms.size(); + code_file.write(reinterpret_cast<char *>(&k),sizeof(k)); + //search for successive and identical blocks + i=k=k0=0; + ModelBlock_Aggregated_Count=-1; + for (j = 0;j < ModelBlock->Size;j++) + { + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) + && (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 )) + { + } + else + { + k=k0=0; + ModelBlock_Aggregated_Count++; + } + k0+=ModelBlock->Block_List[j].Size; + ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0; + ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k; + prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; + } + ModelBlock_Aggregated_Count++; + //For each block + j=0; + for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++) + { + k1=j; + if (k0>0) + code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); + code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK)); + v=ModelBlock_Aggregated_Number[k0]; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + v=ModelBlock->Block_List[j].Simulation_Type; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + for (k=0; k<ModelBlock_Aggregated_Size[k0]; k++) + { + for (i=0; i < ModelBlock->Block_List[j].Size;i++) + { + code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i])); + code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i])); + code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i])); + } + j++; + } + j=k1; + if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || + ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE) + { + code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].is_linear),sizeof(ModelBlock->Block_List[j].is_linear)); + v=block_triangular.ModelBlock->Block_List[j].IM_lead_lag[block_triangular.ModelBlock->Block_List[j].Max_Lag + block_triangular.ModelBlock->Block_List[j].Max_Lead].u_finish + 1; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + v=symbol_table.endo_nbr(); + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + v=block_triangular.ModelBlock->Block_List[j].Max_Lag; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + v=block_triangular.ModelBlock->Block_List[j].Max_Lead; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) + //{ + int u_count_int=0; + Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open, + ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE); + v=u_count_int; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + file_open=true; + //} + } + for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++) + { + //For a block composed of a single equation determines whether we have to evaluate or to solve the equation + if (ModelBlock->Block_List[j].Size==1) + { + lhs_rhs_done=true; + eq_node = equations[ModelBlock->Block_List[j].Equation[0]]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + } + else + lhs_rhs_done=false; + // The equations + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); + //The Temporary terms + temporary_terms_type tt2; +#ifdef DEBUGC + k=0; +#endif + 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++) + { + (*it)->compile(code_file,false, output_type, tt2, map_idx); + code_file.write(&FSTPT, sizeof(FSTPT)); + map_idx_type::const_iterator ii=map_idx.find((*it)->idx); + v=(int)ii->second; + code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); + // 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 = ModelBlock->Block_List[j].Temporary_terms->begin(); + it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) + { + map_idx_type::const_iterator ii=map_idx.find((*it)->idx); + cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n"; + } +#endif + if (!lhs_rhs_done) + { + eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + } + switch (ModelBlock->Block_List[j].Simulation_Type) + { + case EVALUATE_BACKWARD: + case EVALUATE_FORWARD: + rhs->compile(code_file,false, output_type, temporary_terms, map_idx); + lhs->compile(code_file,true, output_type, temporary_terms, map_idx); + break; + case EVALUATE_BACKWARD_R: + case EVALUATE_FORWARD_R: + lhs->compile(code_file,false, output_type, temporary_terms, map_idx); + rhs->compile(code_file,true, output_type, temporary_terms, map_idx); + break; + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + v=ModelBlock->Block_List[j].Equation[i]; + Uf[v].eqr=i; + Uf[v].Ufl=NULL; + goto end; + case SOLVE_TWO_BOUNDARIES_COMPLETE: + case SOLVE_TWO_BOUNDARIES_SIMPLE: + v=ModelBlock->Block_List[j].Equation[i]; + Uf[v].eqr=i; + Uf[v].Ufl=NULL; + goto end; + default: + end: + lhs->compile(code_file,false, output_type, temporary_terms, map_idx); + rhs->compile(code_file,false, output_type, temporary_terms, map_idx); + code_file.write(&FBINARY, sizeof(FBINARY)); + int v=oMinus; + code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); + code_file.write(&FSTPR, sizeof(FSTPR)); + code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); +#ifdef CONDITION + if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE) + output << " condition[" << i << "]=0;\n"; +#endif + } + } + code_file.write(&FENDEQU, sizeof(FENDEQU)); + // The Jacobian if we have to solve the block + 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) + { + switch (ModelBlock->Block_List[j].Simulation_Type) + { + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: + compileDerivative(code_file, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, output_type, map_idx); + code_file.write(&FSTPG, sizeof(FSTPG)); + v=0; + code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); + break; + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + int v=ModelBlock->Block_List[j].Equation[eqr]; + if (!Uf[v].Ufl) + { + Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl_First=Uf[v].Ufl; + } + else + { + Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl=Uf[v].Ufl->pNext; + } + Uf[v].Ufl->pNext=NULL; + Uf[v].Ufl->u=u; + Uf[v].Ufl->var=var; + compileDerivative(code_file, eq, var, 0, output_type, map_idx); + code_file.write(&FSTPU, sizeof(FSTPU)); + code_file.write(reinterpret_cast<char *>(&u), sizeof(u)); + } + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + code_file.write(&FLDR, sizeof(FLDR)); + code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); + code_file.write(&FLDZ, sizeof(FLDZ)); + int v=ModelBlock->Block_List[j].Equation[i]; + for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext) + { + code_file.write(&FLDU, sizeof(FLDU)); + code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); + code_file.write(&FLDV, sizeof(FLDV)); + char vc=eEndogenous; + code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc)); + code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var)); + int v1=0; + code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); + code_file.write(&FBINARY, sizeof(FBINARY)); + v1=oTimes; + code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); + code_file.write(&FCUML, sizeof(FCUML)); + } + 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; + } + code_file.write(&FBINARY, sizeof(FBINARY)); + v=oMinus; + code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); + code_file.write(&FSTPU, sizeof(FSTPU)); + code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); + } + break; + case SOLVE_TWO_BOUNDARIES_COMPLETE: + case SOLVE_TWO_BOUNDARIES_SIMPLE: + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + int v=ModelBlock->Block_List[j].Equation[eqr]; + if (!Uf[v].Ufl) + { + Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl_First=Uf[v].Ufl; + } + else + { + Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl=Uf[v].Ufl->pNext; + } + Uf[v].Ufl->pNext=NULL; + Uf[v].Ufl->u=u; + Uf[v].Ufl->var=var; + Uf[v].Ufl->lag=k; + compileDerivative(code_file, eq, var, k, output_type, map_idx); + code_file.write(&FSTPU, sizeof(FSTPU)); + code_file.write(reinterpret_cast<char *>(&u), sizeof(u)); +#ifdef CONDITION + output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n"; + output << " condition[" << eqr << "]=u[" << u << "+Per_u_];\n"; +#endif + } + } + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + code_file.write(&FLDR, sizeof(FLDR)); + code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); + code_file.write(&FLDZ, sizeof(FLDZ)); + int v=ModelBlock->Block_List[j].Equation[i]; + for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext) + { + code_file.write(&FLDU, sizeof(FLDU)); + code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); + code_file.write(&FLDV, sizeof(FLDV)); + char vc=eEndogenous; + code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc)); + int v1=Uf[v].Ufl->var; + code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); + v1=Uf[v].Ufl->lag; + code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); + code_file.write(&FBINARY, sizeof(FBINARY)); + v1=oTimes; + code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); + code_file.write(&FCUML, sizeof(FCUML)); + } + 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; + } + code_file.write(&FBINARY, sizeof(FBINARY)); + v=oMinus; + code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); + code_file.write(&FSTPU, sizeof(FSTPU)); + code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); +#ifdef CONDITION + output << " if (fabs(condition[" << i << "])<fabs(u[" << i << "+Per_u_]))\n"; + output << " condition[" << i << "]=u[" << i << "+Per_u_];\n"; +#endif + } +#ifdef CONDITION + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + output << " u[" << u << "+Per_u_] /= condition[" << eqr << "];\n"; + } + } + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n"; +#endif + break; + default: + break; + } + + prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; + } + j++; + } + } + code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); + code_file.write(&FEND, sizeof(FEND)); + 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; + + writeDynamicModel(mDynamicModelFile); + + 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; + + // Writing the function body + writeDynamicModel(mDynamicModelFile); + + // 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, *g2;" << 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() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl + << " /* Create a C pointer to a copy of the output matrix g1. */" << endl + << " g1 = mxGetPr(plhs[1]);" << endl + << " }" << endl + << endl + << " g2 = NULL;" << endl + << " if (nlhs >= 3)" << endl + << " {" << endl + << " /* Set the output pointer to the output matrix g2. */" << endl; + int g2_ncols = variable_table.getDynJacobianColsNbr(computeJacobianExo)*variable_table.getDynJacobianColsNbr(computeJacobianExo); + mDynamicModelFile << " plhs[2] = mxCreateSparse(" << equations.size() << ", " << g2_ncols << ", " + << 5*g2_ncols << ", mxREAL);" << endl + << " /* Create a C pointer to a copy of the output matrix g1. */" << endl + << " g2 = mxGetPr(plhs[2]);" << endl + << " }" << endl + << endl + << " /* Call the C subroutines. */" << endl + << " Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);" << endl + << "}" << endl; + mDynamicModelFile.close(); +} + +string +DynamicModel::reform(const string name1) const +{ + string name=name1; + int pos = name.find("\\", 0); + while (pos >= 0) + { + if (name.substr(pos + 1, 1) != "\\") + { + name = name.insert(pos, "\\"); + pos++; + } + pos++; + pos = name.find("\\", pos); + } + return (name); +} + +void +DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num, + int &u_count_int, bool &file_open, bool is_two_boundaries) const +{ + int j; + std::ofstream SaveCode; + if (file_open) + SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate ); + else + SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::binary); + if (!SaveCode.is_open()) + { + cout << "Error : Can't open file \"" << bin_basename << ".bin\" for writing\n"; + exit(EXIT_FAILURE); + } + u_count_int=0; + for (int m=0;m<=block_triangular.ModelBlock->Block_List[num].Max_Lead+block_triangular.ModelBlock->Block_List[num].Max_Lag;m++) + { + int k1=m-block_triangular.ModelBlock->Block_List[num].Max_Lag; + for (j=0;j<block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].size;j++) + { + int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size; + int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j]; + int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j]; + SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); + SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr)); + SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1)); + SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u)); + u_count_int++; + } + } + if(is_two_boundaries) + { + for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++) + { + int eqr1=j; + int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods + +block_triangular.incidencematrix.Model_Max_Lead_Endo); + int k1=0; + SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); + SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr)); + SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1)); + SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); + u_count_int++; + } + } + //cout << "u_count_int=" << u_count_int << "\n"; + for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++) + { + int varr=block_triangular.ModelBlock->Block_List[num].Variable[j]; + SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr)); + } + for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++) + { + int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j]; + SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); + } + SaveCode.close(); +} + +void +DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const +{ + string sp; + ofstream mDynamicModelFile; + ostringstream tmp, tmp1, tmp_eq; + int prev_Simulation_Type, tmp_i; + //SymbolicGaussElimination SGE; + bool OK; + chdir(basename.c_str()); + string filename = dynamic_basename + ".m"; + 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 << "%\n"; + mDynamicModelFile << "% " << filename << " : Computes dynamic model for Dynare\n"; + mDynamicModelFile << "%\n"; + mDynamicModelFile << "% Warning : this file is generated automatically by Dynare\n"; + mDynamicModelFile << "% from model file (.mod)\n\n"; + mDynamicModelFile << "%/\n"; + + int i, k, Nb_SGE=0; + bool skip_head, open_par=false; + if (computeJacobian || computeJacobianExo || computeHessian) + { + mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n"; + mDynamicModelFile << " global oo_ options_ M_ ;\n"; + mDynamicModelFile << " g2=[];g3=[];\n"; + //Temporary variables declaration + OK=true; + ostringstream tmp_output; + for (temporary_terms_type::const_iterator it = temporary_terms.begin(); + it != temporary_terms.end(); it++) + { + if (OK) + OK=false; + else + tmp_output << " "; + (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms); + } + if (tmp_output.str().length()>0) + mDynamicModelFile << " global " << tmp_output.str() << " M_ ;\n"; + + mDynamicModelFile << " T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n"; + tmp_output.str(""); + for (temporary_terms_type::const_iterator it = temporary_terms.begin(); + it != temporary_terms.end(); it++) + { + tmp_output << " "; + (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); + tmp_output << "=T_init;\n"; + } + if (tmp_output.str().length()>0) + mDynamicModelFile << tmp_output.str(); + + mDynamicModelFile << " y_kmin=M_.maximum_lag;\n"; + mDynamicModelFile << " y_kmax=M_.maximum_lead;\n"; + mDynamicModelFile << " y_size=M_.endo_nbr;\n"; + mDynamicModelFile << " if(length(varargin)>0)\n"; + mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n"; + mDynamicModelFile << " params=varargin{3};\n"; + mDynamicModelFile << " it_=varargin{4};\n"; + /*i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+ + symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1); + mDynamicModelFile << " g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n";*/ + mDynamicModelFile << " Per_u_=0;\n"; + mDynamicModelFile << " Per_y_=it_*y_size;\n"; + mDynamicModelFile << " y=varargin{1};\n"; + mDynamicModelFile << " ys=y(it_,:);\n"; + mDynamicModelFile << " x=varargin{2};\n"; + prev_Simulation_Type=-1; + tmp.str(""); + tmp_eq.str(""); + for (int count_call=1, i = 0;i < block_triangular.ModelBlock->Size;i++, count_call++) + { + k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; + if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) && + ((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R) + || (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))) + { + mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; + tmp_eq.str(""); + mDynamicModelFile << " y_index=[" << tmp.str() << "];\n"; + tmp.str(""); + mDynamicModelFile << tmp1.str(); + tmp1.str(""); + } + for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) + { + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; + } + if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R) + { + if (i==block_triangular.ModelBlock->Size-1) + { + mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; + tmp_eq.str(""); + mDynamicModelFile << " y_index=[" << tmp.str() << "];\n"; + tmp.str(""); + mDynamicModelFile << tmp1.str(); + tmp1.str(""); + } + } + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && + (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) + skip_head=true; + else + skip_head=false; + switch (k) + { + case EVALUATE_FORWARD: + case EVALUATE_BACKWARD: + case EVALUATE_FORWARD_R: + case EVALUATE_BACKWARD_R: + if (!skip_head) + { + tmp1 << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n"; + tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n"; + } + break; + case SOLVE_FORWARD_SIMPLE: + case SOLVE_BACKWARD_SIMPLE: + mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n"; + mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n"; + mDynamicModelFile << " residual(y_index_eq)=r;\n"; + tmp_eq.str(""); + tmp.str(""); + break; + case SOLVE_FORWARD_COMPLETE: + case SOLVE_BACKWARD_COMPLETE: + mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; + mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n"; + mDynamicModelFile << " residual(y_index_eq)=r;\n"; + break; + case SOLVE_TWO_BOUNDARIES_COMPLETE: + case SOLVE_TWO_BOUNDARIES_SIMPLE: + int j; + mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; + mDynamicModelFile << " y_index = ["; + for (j=0;j<tmp_i;j++) + for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) + { + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr(); + } + int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1; + for (j=0;j<tmp_ix;j++) + for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].nb_exo;ik++) + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr()+symbol_table.endo_nbr()*tmp_i; + mDynamicModelFile << " ];\n"; + tmp.str(""); + tmp_eq.str(""); + //mDynamicModelFile << " ga = [];\n"; + j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1) + + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1); + /*mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " << + block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";*/ + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; + mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << variable_table.max_lag << ", 1, " << variable_table.max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n"; + /*if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead) + mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n"; + else + mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";*/ + mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n"; + break; + } + prev_Simulation_Type=k; + } + if (tmp1.str().length()) + { + mDynamicModelFile << tmp1.str(); + tmp1.str(""); + } + mDynamicModelFile << " varargout{1}=residual;\n"; + mDynamicModelFile << " varargout{2}=dr;\n"; + mDynamicModelFile << " return;\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n"; + mDynamicModelFile << " if(options_.simulation_method==0)\n"; + mDynamicModelFile << " mthd='Sparse LU';\n"; + mDynamicModelFile << " elseif(options_.simulation_method==2)\n"; + mDynamicModelFile << " mthd='GMRES';\n"; + mDynamicModelFile << " elseif(options_.simulation_method==3)\n"; + mDynamicModelFile << " mthd='BICGSTAB';\n"; + mDynamicModelFile << " else\n"; + mDynamicModelFile << " mthd='UNKNOWN';\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n"; + mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n"; + mDynamicModelFile << " fprintf('\\n') ;\n"; + mDynamicModelFile << " periods=options_.periods;\n"; + mDynamicModelFile << " maxit_=options_.maxit_;\n"; + mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n"; + mDynamicModelFile << " y=oo_.endo_simul';\n"; + mDynamicModelFile << " x=oo_.exo_simul;\n"; + } + prev_Simulation_Type=-1; + mDynamicModelFile << " params=M_.params;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.status = 0;\n"; + for (i = 0;i < block_triangular.ModelBlock->Size;i++) + { + k = block_triangular.ModelBlock->Block_List[i].Simulation_Type; + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && + (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) + skip_head=true; + else + skip_head=false; + if ((k == EVALUATE_FORWARD || k == EVALUATE_FORWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + if (!skip_head) + { + if (open_par) + { + mDynamicModelFile << " end\n"; + } + mDynamicModelFile << " oo_.deterministic_simulation.status = 1;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.error = 0;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.iterations = 0;\n"; + mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mDynamicModelFile << " else\n"; + mDynamicModelFile << " blck_num = 1;\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).status = 1;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n"; + mDynamicModelFile << " g1=[];g2=[];g3=[];\n"; + //mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n"; + mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n"; + } + //open_par=true; + } + else if ((k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + if (!skip_head) + { + if (open_par) + { + mDynamicModelFile << " end\n"; + } + mDynamicModelFile << " oo_.deterministic_simulation.status = 1;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.error = 0;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.iterations = 0;\n"; + mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mDynamicModelFile << " else\n"; + mDynamicModelFile << " blck_num = 1;\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).status = 1;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n"; + mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n"; + mDynamicModelFile << " g1=[];g2=[];g3=[];\n"; + mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n"; + } + } + else if ((k == SOLVE_FORWARD_COMPLETE || k == SOLVE_FORWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + if (open_par) + mDynamicModelFile << " end\n"; + open_par=false; + mDynamicModelFile << " g1=0;\n"; + mDynamicModelFile << " r=0;\n"; + tmp.str(""); + for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) + { + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + } + mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; + int nze, m; + for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) + nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; + mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mDynamicModelFile << " else\n"; + mDynamicModelFile << " blck_num = 1;\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" << + ", y, x, params, y_index, " << nze << + ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear << + ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n"; + + } + else if ((k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + if (open_par) + mDynamicModelFile << " end\n"; + open_par=false; + mDynamicModelFile << " g1=0;\n"; + mDynamicModelFile << " r=0;\n"; + tmp.str(""); + for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) + { + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + } + mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; + int nze, m; + for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) + nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; + mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mDynamicModelFile << " else\n"; + mDynamicModelFile << " blck_num = 1;\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" << + ", y, x, params, y_index, " << nze << + ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear << + ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n"; + } + else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE || k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + if (open_par) + mDynamicModelFile << " end\n"; + open_par=false; + Nb_SGE++; + int nze, m; + for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) + nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; + mDynamicModelFile << " y_index=["; + for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) + { + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + } + mDynamicModelFile << " ];\n"; + mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mDynamicModelFile << " else\n"; + mDynamicModelFile << " blck_num = 1;\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " y = solve_two_boundaries('" << dynamic_basename << "_" << i + 1 << "'" << + ", y, x, params, y_index, " << nze << + ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].Max_Lag << + ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << + ", " << block_triangular.ModelBlock->Block_List[i].is_linear << + ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method);\n"; + + } + prev_Simulation_Type=k; + } + if (open_par) + mDynamicModelFile << " end;\n"; + open_par=false; + mDynamicModelFile << " oo_.endo_simul = y';\n"; + mDynamicModelFile << "return;\n"; + + mDynamicModelFile.close(); + + writeModelEquationsOrdered_M( block_triangular.ModelBlock, dynamic_basename); + + chdir(".."); +} + +void +DynamicModel::writeDynamicModel(ostream &DynamicOutput) const +{ + ostringstream lsymetric; // Used when writing symetric elements in Hessian + ostringstream model_output; // Used for storing model equations + ostringstream jacobian_output; // Used for storing jacobian equations + ostringstream hessian_output; // Used for storing Hessian equations + ostringstream third_derivatives_output; + + ExprNodeOutputType output_type = (mode == eStandardMode || mode==eSparseMode ? oMatlabDynamicModel : oCDynamicModel); + + writeModelLocalVariables(model_output, output_type); + + writeTemporaryTerms(model_output, output_type); + + writeModelEquations(model_output, output_type); + + int nrows = equations.size(); + int nvars = variable_table.getDynJacobianColsNbr(computeJacobianExo); + int nvars_sq = nvars * nvars; + + // Writing Jacobian + if (computeJacobian || computeJacobianExo) + for (first_derivatives_type::const_iterator it = first_derivatives.begin(); + it != first_derivatives.end(); it++) + { + int eq = it->first.first; + int var = it->first.second; + NodeID d1 = it->second; + + if (computeJacobianExo || variable_table.getType(var) == eEndogenous) + { + ostringstream g1; + g1 << " g1"; + matrixHelper(g1, eq, variable_table.getDynJacobianCol(var), output_type); + + jacobian_output << g1.str() << "=" << g1.str() << "+"; + d1->writeOutput(jacobian_output, output_type, temporary_terms); + jacobian_output << ";" << endl; + } + } + + // Writing Hessian + if (computeHessian) + for (second_derivatives_type::const_iterator it = second_derivatives.begin(); + it != second_derivatives.end(); it++) + { + int eq = it->first.first; + int var1 = it->first.second.first; + int var2 = it->first.second.second; + NodeID d2 = it->second; + + int id1 = variable_table.getDynJacobianCol(var1); + int id2 = variable_table.getDynJacobianCol(var2); + + int col_nb = id1*nvars+id2; + int col_nb_sym = id2*nvars+id1; + + hessian_output << " g2"; + matrixHelper(hessian_output, eq, col_nb, output_type); + hessian_output << " = "; + d2->writeOutput(hessian_output, output_type, temporary_terms); + hessian_output << ";" << endl; + + // Treating symetric elements + if (id1 != id2) + { + lsymetric << " g2"; + matrixHelper(lsymetric, eq, col_nb_sym, output_type); + lsymetric << " = " << "g2"; + matrixHelper(lsymetric, eq, col_nb, output_type); + lsymetric << ";" << endl; + } + } + + // Writing third derivatives + if (computeThirdDerivatives) + for (third_derivatives_type::const_iterator it = third_derivatives.begin(); + it != third_derivatives.end(); it++) + { + int eq = it->first.first; + int var1 = it->first.second.first; + int var2 = it->first.second.second.first; + int var3 = it->first.second.second.second; + NodeID d3 = it->second; + + int id1 = variable_table.getDynJacobianCol(var1); + int id2 = variable_table.getDynJacobianCol(var2); + int id3 = variable_table.getDynJacobianCol(var3); + + // Reference column number for the g3 matrix + int ref_col = id1 * nvars_sq + id2 * nvars + id3; + + third_derivatives_output << " g3"; + matrixHelper(third_derivatives_output, eq, ref_col, output_type); + third_derivatives_output << " = "; + d3->writeOutput(third_derivatives_output, output_type, temporary_terms); + third_derivatives_output << ";" << endl; + + // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal) + set<int> cols; + cols.insert(id1 * nvars_sq + id3 * nvars + id2); + cols.insert(id2 * nvars_sq + id1 * nvars + id3); + cols.insert(id2 * nvars_sq + id3 * nvars + id1); + cols.insert(id3 * nvars_sq + id1 * nvars + id2); + cols.insert(id3 * nvars_sq + id2 * nvars + id1); + + for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++) + if (*it2 != ref_col) + { + third_derivatives_output << " g3"; + matrixHelper(third_derivatives_output, eq, *it2, output_type); + third_derivatives_output << " = " << "g3"; + matrixHelper(third_derivatives_output, eq, ref_col, output_type); + third_derivatives_output << ";" << endl; + } + } + + if (mode == eStandardMode) + { + DynamicOutput << "%" << endl + << "% Model equations" << endl + << "%" << endl + << endl + << "residual = zeros(" << nrows << ", 1);" << endl + << model_output.str(); + + if (computeJacobian || computeJacobianExo) + { + // Writing initialization instruction for matrix g1 + DynamicOutput << "if nargout >= 2," << endl + << " g1 = zeros(" << nrows << ", " << nvars << ");" << endl + << endl + << "%" << endl + << "% Jacobian matrix" << endl + << "%" << endl + << endl + << jacobian_output.str() + << "end" << endl; + } + if (computeHessian) + { + // Writing initialization instruction for matrix g2 + int ncols = nvars_sq; + DynamicOutput << "if nargout >= 3," << endl + << " g2 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl + << endl + << "%" << endl + << "% Hessian matrix" << endl + << "%" << endl + << endl + << hessian_output.str() + << lsymetric.str() + << "end;" << endl; + } + if (computeThirdDerivatives) + { + int ncols = nvars_sq * nvars; + DynamicOutput << "if nargout >= 4," << endl + << " g3 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl + << endl + << "%" << endl + << "% Third order derivatives" << endl + << "%" << endl + << endl + << third_derivatives_output.str() + << "end;" << endl; + } + } + else + { + DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *g2)" << endl + << "{" << endl + << " double lhs, rhs;" << endl + << endl + << " /* Residual equations */" << endl + << model_output.str(); + + if (computeJacobian || computeJacobianExo) + { + DynamicOutput << " /* Jacobian */" << endl + << " if (g1 == NULL)" << endl + << " return;" << endl + << " else" << endl + << " {" << endl + << jacobian_output.str() + << " }" << endl; + } + if (computeHessian) + { + DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl + << " if (g2 == NULL)" << endl + << " return;" << endl + << " else" << endl + << " {" << endl + << hessian_output.str() + << lsymetric.str() + << " }" << endl; + } + DynamicOutput << "}" << endl << endl; + } +} + +void +DynamicModel::writeOutput(ostream &output) const +{ + /* Writing initialisation for M_.lead_lag_incidence matrix + M_.lead_lag_incidence is a matrix with as many columns as there are + endogenous variables and as many rows as there are periods in the + models (nbr of rows = M_.max_lag+M_.max_lead+1) + + The matrix elements are equal to zero if a variable isn't present in the + model at a given period. + */ + output << "M_.lead_lag_incidence = ["; + // Loop on endogenous variables + int lag = 0; + for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++) + { + output << "\n\t"; + // Loop on periods + for (lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++) + { + // Print variableID if exists with current period, otherwise print 0 + try + { + int varID = variable_table.getID(symbol_table.getID(eEndogenous, endoID), lag); + output << " " << variable_table.getDynJacobianCol(varID) + 1; + } + catch (VariableTable::UnknownVariableKeyException &e) + { + output << " 0"; + } + } + output << ";"; + } + output << "]';\n"; + //In case of sparse model, writes the block structure of the model + + if (mode==eSparseMode || mode==eSparseDLLMode) + { + //int prev_Simulation_Type=-1; + //bool skip_the_head; + int k=0; + int count_lead_lag_incidence = 0; + int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo; + for (int j = 0;j < block_triangular.ModelBlock->Size;j++) + { + //For a block composed of a single equation determines wether we have to evaluate or to solve the equation + //skip_the_head=false; + k++; + count_lead_lag_incidence = 0; + int Block_size=block_triangular.ModelBlock->Block_List[j].Size; + max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ; + max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead; + max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ; + max_lead_endo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Endo; + max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ; + max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo; + bool evaluate=false; + vector<int> exogenous; + vector<int>::iterator it_exogenous; + exogenous.clear(); + ostringstream tmp_s, tmp_s_eq; + tmp_s.str(""); + tmp_s_eq.str(""); + for (int i=0;i<block_triangular.ModelBlock->Block_List[j].Size;i++) + { + tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1; + tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1; + } + for (int i=0;i<block_triangular.ModelBlock->Block_List[j].nb_exo;i++) + { + int ii=block_triangular.ModelBlock->Block_List[j].Exogenous[i]; + for (it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) /*cout << "*it_exogenous=" << *it_exogenous << "\n"*/; + if (it_exogenous==exogenous.end() || exogenous.begin()==exogenous.end()) + exogenous.push_back(ii); + } + output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n"; + output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_lead = " << max_lead << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag_endo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead_endo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_exo_lag = " << max_lag_exo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_exo_lead = " << max_lead_exo << ";\n"; + output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n"; + output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n"; + output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n"; + output << "M_.block_structure.block(" << k << ").exogenous = ["; + int i=0; + for (it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) + if (*it_exogenous>=0) + { + output << " " << *it_exogenous+1; + i++; + } + output << "];\n"; + output << "M_.block_structure.block(" << k << ").exo_nbr = " << i << ";\n"; + + output << "M_.block_structure.block(" << k << ").exo_det_nbr = " << block_triangular.ModelBlock->Block_List[j].nb_exo_det << ";\n"; + + tmp_s.str(""); + + bool done_IM=false; + if (!evaluate) + { + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n"; + for (int l=-max_lag_endo;l<max_lead_endo+1;l++) + { + bool *tmp_IM; + tmp_IM=block_triangular.incidencematrix.Get_IM(l, eEndogenous); + if (tmp_IM) + { + for (int l_var=0;l_var<block_triangular.ModelBlock->Block_List[j].Size;l_var++) + { + for (int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[j].Size;l_equ++) + if (tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr()+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) + { + count_lead_lag_incidence++; + if (tmp_s.str().length()) + tmp_s << " "; + tmp_s << count_lead_lag_incidence; + done_IM=true; + break; + } + if (!done_IM) + tmp_s << " 0"; + done_IM=false; + } + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; + tmp_s.str(""); + } + } + } + else + { + bool done_some_where; + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; + for (int l=-max_lag_endo;l<max_lead_endo+1;l++) + { + bool not_increm=true; + bool *tmp_IM; + tmp_IM=block_triangular.incidencematrix.Get_IM(l, eEndogenous); + int ii=j; + if (tmp_IM) + { + done_some_where = false; + while (ii-j<Block_size) + { + for (int l_var=0;l_var<block_triangular.ModelBlock->Block_List[ii].Size;l_var++) + { + for (int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[ii].Size;l_equ++) + if (tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr()+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) + { + //if(not_increm && l==-max_lag) + count_lead_lag_incidence++; + not_increm=false; + if (tmp_s.str().length()) + tmp_s << " "; + //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; + tmp_s << count_lead_lag_incidence; + done_IM=true; + break; + } + if (!done_IM) + tmp_s << " 0"; + else + done_some_where = true; + done_IM=false; + } + ii++; + } + output << tmp_s.str() << "\n"; + tmp_s.str(""); + } + } + output << "];\n"; + } + + } + for (int j=-block_triangular.incidencematrix.Model_Max_Lag_Endo;j<=block_triangular.incidencematrix.Model_Max_Lead_Endo;j++) + { + bool* IM = block_triangular.incidencematrix.Get_IM(j, eEndogenous); + if (IM) + { + bool new_entry=true; + output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").lead_lag = " << j << ";\n"; + output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").sparse_IM = ["; + for (int i=0;i<symbol_table.endo_nbr()*symbol_table.endo_nbr();i++) + { + if (IM[i]) + { + if (!new_entry) + output << " ; "; + else + output << " "; + output << i/symbol_table.endo_nbr()+1 << " " << i % symbol_table.endo_nbr()+1; + new_entry=false; + } + } + output << "];\n"; + } + } + } + // Writing initialization for some other variables + output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];\n"; + output << "M_.maximum_lag = " << variable_table.max_lag << ";\n"; + output << "M_.maximum_lead = " << variable_table.max_lead << ";\n"; + if (symbol_table.endo_nbr()) + { + output << "M_.maximum_endo_lag = " << variable_table.max_endo_lag << ";\n"; + output << "M_.maximum_endo_lead = " << variable_table.max_endo_lead << ";\n"; + output << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);\n"; + } + if (symbol_table.exo_nbr()) + { + output << "M_.maximum_exo_lag = " << variable_table.max_exo_lag << ";\n"; + output << "M_.maximum_exo_lead = " << variable_table.max_exo_lead << ";\n"; + output << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);\n"; + } + if (symbol_table.exo_det_nbr()) + { + output << "M_.maximum_exo_det_lag = " << variable_table.max_exo_det_lag << ";\n"; + output << "M_.maximum_exo_det_lead = " << variable_table.max_exo_det_lead << ";\n"; + output << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);\n"; + } + if (symbol_table.param_nbr()) + output << "M_.params = repmat(NaN," << symbol_table.param_nbr() << ", 1);\n"; +} + +void +DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m) +{ + int i=0; + int j=0; + bool *IM=NULL; + int a_variable_lag=-9999; + for (first_derivatives_type::iterator it = first_derivatives.begin(); + it != first_derivatives.end(); it++) + { + //cout << "it->first.second=" << it->first.second << " variable_table.getSymbolID(it->first.second)=" << variable_table.getSymbolID(it->first.second) << " Type=" << variable_table.getType(it->first.second) << " eEndogenous=" << eEndogenous << " eExogenous=" << eExogenous << " variable_table.getLag(it->first.second)=" << variable_table.getLag(it->first.second) << "\n"; + if (variable_table.getType(it->first.second) == eEndogenous) + { + NodeID Id = it->second; + double val = 0; + try + { + val = Id->eval(eval_context); + } + catch (ExprNode::EvalException &e) + { + cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ") [" << variable_table.getSymbolID(it->first.second) << "] !" << endl; + Id->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms); + cout << "\n"; + cerr << "DynamicModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ")!" << endl; + } + int eq=it->first.first; + int var=symbol_table.getTypeSpecificID(variable_table.getSymbolID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second); + int k1=variable_table.getLag(it->first.second); + if (a_variable_lag!=k1) + { + IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous); + a_variable_lag=k1; + } + if (k1==0) + { + j++; + (*j_m)[make_pair(eq,var)]=val; + } + if (IM[eq*symbol_table.endo_nbr()+var] && (fabs(val) < cutoff)) + { + if (block_triangular.bt_verbose) + cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")\n"; + block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous); + i++; + } + } + } + if (i>0) + { + cout << i << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded\n"; + cout << "the contemporaneous incidence matrix has " << j << " elements\n"; + } +} + +void +DynamicModel::BlockLinear(Model_Block *ModelBlock) +{ + int i,j,l,m,ll; + for (j = 0;j < ModelBlock->Size;j++) + { + if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || + ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE) + { + ll=ModelBlock->Block_List[j].Max_Lag; + for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[ll].size;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[ll].Equ_Index[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[ll].Var_Index[i]; + //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,0))); + first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var),0))); + if (it!= first_derivatives.end()) + { + NodeID Id = it->second; + set<pair<int, int> > endogenous; + Id->collectEndogenous(endogenous); + if (endogenous.size() > 0) + { + for (l=0;l<ModelBlock->Block_List[j].Size;l++) + { + if (endogenous.find(make_pair(ModelBlock->Block_List[j].Variable[l], 0)) != endogenous.end()) + { + ModelBlock->Block_List[j].is_linear=false; + goto follow; + } + } + } + } + } + } + else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) + { + for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + int k1=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]; + //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,k1))); + first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var),k1))); + NodeID Id = it->second; + if (it!= first_derivatives.end()) + { + set<pair<int, int> > endogenous; + Id->collectEndogenous(endogenous); + if (endogenous.size() > 0) + { + for (l=0;l<ModelBlock->Block_List[j].Size;l++) + { + if (endogenous.find(make_pair(ModelBlock->Block_List[j].Variable[l], k1)) != endogenous.end()) + { + ModelBlock->Block_List[j].is_linear=false; + goto follow; + } + } + } + } + } + } + } + follow: + i=0; + } +} + +void +DynamicModel::computingPass(const eval_context_type &eval_context, bool no_tmp_terms) +{ + // Computes dynamic jacobian columns + variable_table.computeDynJacobianCols(); + + // Determine derivation order + int order = 1; + if (computeThirdDerivatives) + order = 3; + else if (computeHessian) + order = 2; + + // Launch computations + cout << "Computing dynamic model derivatives at order " << order << "..."; + derive(order); + cout << " done" << endl; + + if (mode == eSparseDLLMode || mode == eSparseMode) + { + BuildIncidenceMatrix(); + + jacob_map j_m; + + evaluateJacobian(eval_context, &j_m); + + + if (block_triangular.bt_verbose) + { + cout << "The gross incidence matrix \n"; + block_triangular.incidencematrix.Print_IM(eEndogenous); + } + block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations); + BlockLinear(block_triangular.ModelBlock); + if (!no_tmp_terms) + computeTemporaryTermsOrdered(order, block_triangular.ModelBlock); + } + else + if (!no_tmp_terms) + computeTemporaryTerms(order); +} + +void +DynamicModel::writeDynamicFile(const string &basename) const +{ + switch (mode) + { + case eStandardMode: + writeDynamicMFile(basename + "_dynamic"); + break; + case eSparseMode: + writeSparseDynamicMFile(basename + "_dynamic", basename, mode); + block_triangular.Free_Block(block_triangular.ModelBlock); + block_triangular.incidencematrix.Free_IM(); + //block_triangular.Free_IM_X(block_triangular.First_IM_X); + break; + case eDLLMode: + writeDynamicCFile(basename + "_dynamic"); + break; + case eSparseDLLMode: + // create a directory to store all the files +#ifdef _WIN32 + mkdir(basename.c_str()); +#else + mkdir(basename.c_str(), 0777); +#endif + writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL, map_idx); + block_triangular.Free_Block(block_triangular.ModelBlock); + block_triangular.incidencematrix.Free_IM(); + //block_triangular.Free_IM_X(block_triangular.First_IM_X); + break; + } +} + +void +DynamicModel::toStatic(StaticModel &static_model) const +{ + static_model.mode = mode; + + // Convert model local variables (need to be done first) + for (map<int, NodeID>::const_iterator it = local_variables_table.begin(); + it != local_variables_table.end(); it++) + static_model.AddLocalParameter(symbol_table.getName(it->first), it->second->toStatic(static_model)); + + // Convert equations + for (vector<BinaryOpNode *>::const_iterator it = equations.begin(); + it != equations.end(); it++) + static_model.addEquation((*it)->toStatic(static_model)); +} + diff --git a/DynamicModel.hh b/DynamicModel.hh new file mode 100644 index 0000000000000000000000000000000000000000..10b434ed7bb70408c93267eab8dd791cfe671bbb --- /dev/null +++ b/DynamicModel.hh @@ -0,0 +1,96 @@ +/* + * 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/>. + */ + +#ifndef _DYNAMICMODEL_HH +#define _DYNAMICMODEL_HH + +using namespace std; + +#include "StaticModel.hh" +#include "BlockTriangular.hh" + +//! Stores a dynamic model +class DynamicModel : public ModelTree +{ +private: + //! Writes dynamic model file (Matlab version) + void writeDynamicMFile(const string &dynamic_basename) const; + //! Writes dynamic model file (C version) + /*! \todo add third derivatives handling */ + void writeDynamicCFile(const string &dynamic_basename) const; + //! Writes dynamic model file when SparseDLL option is on + void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const; + //! Writes the dynamic model equations and its derivatives + /*! \todo add third derivatives handling in C output */ + void writeDynamicModel(ostream &DynamicOutput) const; + //! Writes the Block reordred structure of the model in M output + void writeModelEquationsOrdered_M(Model_Block *ModelBlock, const string &dynamic_basename) const; + //! Writes the code of the Block reordred structure of the model in virtual machine bytecode + void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const; + //! Computes jacobian and prepares for equation normalization + /*! Using values from initval/endval blocks and parameter initializations: + - computes the jacobian for the model w.r. to contemporaneous variables + - removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff) + */ + void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m); + void BlockLinear(Model_Block *ModelBlock); + string reform(string name) const; + map_idx_type map_idx; + //! Build The incidence matrix form the modeltree + void BuildIncidenceMatrix(); + + void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock); + //! Write derivative code of an equation w.r. to a variable + void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const; + +public: + DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); + //! Absolute value under which a number is considered to be zero + double cutoff; + //! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc) + double markowitz; + //! the file containing the model and the derivatives code + ofstream code_file; + //! Whether dynamic Jacobian (w.r. to endogenous) should be written + bool computeJacobian; + //! Whether dynamic Jacobian (w.r. to endogenous and exogenous) should be written + bool computeJacobianExo; + //! Whether dynamic Hessian (w.r. to endogenous and exogenous) should be written + bool computeHessian; + //! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written + bool computeThirdDerivatives; + //! Execute computations (variable sorting + derivation) + /*! You must set computeJacobian, computeJacobianExo, computeHessian and computeThirdDerivatives to correct values before calling this function + \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */ + void computingPass(const eval_context_type &eval_context, bool no_tmp_terms); + //! Writes model initialization and lead/lag incidence matrix to output + void writeOutput(ostream &output) const; + //! Complete set to block decompose the model + BlockTriangular block_triangular; + //! Adds informations for simulation in a binary file + void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, + const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const; + //! Writes dynamic model file + void writeDynamicFile(const string &basename) const; + //! Converts to static model (only the equations) + /*! It assumes that the static model given in argument has just been allocated */ + void toStatic(StaticModel &static_model) const; +}; + +#endif diff --git a/ExprNode.cc b/ExprNode.cc index 57ada61a54279c9d4479c15459de262f36846a8b..b4918adbcbff2e39e367d30ae5cc3d59a4ac278b 100644 --- a/ExprNode.cc +++ b/ExprNode.cc @@ -99,6 +99,7 @@ ExprNode::writeOutput(ostream &output) writeOutput(output, oMatlabOutsideModel, temporary_terms_type()); } + NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) : ExprNode(datatree_arg), id(id_arg) @@ -115,7 +116,6 @@ NumConstNode::computeDerivative(int varID) return datatree.Zero; } - void NumConstNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const { @@ -158,7 +158,6 @@ NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou CompileCode.write(reinterpret_cast<char *>(&vard),sizeof(vard)); } - void NumConstNode::collectEndogenous(set<pair<int, int> > &result) const { @@ -169,6 +168,12 @@ NumConstNode::collectExogenous(set<pair<int, int> > &result) const { } +NodeID +NumConstNode::toStatic(DataTree &static_datatree) const +{ + return static_datatree.AddNumConstant(datatree.num_constants.get(id)); +} + VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg) : ExprNode(datatree_arg), @@ -479,7 +484,6 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou } } - void VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, @@ -493,9 +497,6 @@ VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count, datatree.local_variables_table[symb_id]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); } - - - void VariableNode::collectEndogenous(set<pair<int, int> > &result) const { @@ -514,6 +515,12 @@ VariableNode::collectExogenous(set<pair<int, int> > &result) const datatree.local_variables_table[symb_id]->collectExogenous(result); } +NodeID +VariableNode::toStatic(DataTree &static_datatree) const +{ + return static_datatree.AddVariable(datatree.symbol_table.getName(symb_id)); +} + UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) : ExprNode(datatree_arg), @@ -923,6 +930,49 @@ UnaryOpNode::collectExogenous(set<pair<int, int> > &result) const arg->collectExogenous(result); } +NodeID +UnaryOpNode::toStatic(DataTree &static_datatree) const +{ + NodeID sarg = arg->toStatic(static_datatree); + switch(op_code) + { + case oUminus: + return static_datatree.AddUMinus(sarg); + case oExp: + return static_datatree.AddExp(sarg); + case oLog: + return static_datatree.AddLog(sarg); + case oLog10: + return static_datatree.AddLog10(sarg); + case oCos: + return static_datatree.AddCos(sarg); + case oSin: + return static_datatree.AddSin(sarg); + case oTan: + return static_datatree.AddTan(sarg); + case oAcos: + return static_datatree.AddACos(sarg); + case oAsin: + return static_datatree.AddASin(sarg); + case oAtan: + return static_datatree.AddATan(sarg); + case oCosh: + return static_datatree.AddCosH(sarg); + case oSinh: + return static_datatree.AddSinH(sarg); + case oTanh: + return static_datatree.AddTanH(sarg); + case oAcosh: + return static_datatree.AddACosH(sarg); + case oAsinh: + return static_datatree.AddASinH(sarg); + case oAtanh: + return static_datatree.AddATanH(sarg); + case oSqrt: + return static_datatree.AddSqRt(sarg); + } +} + BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, BinaryOpcode op_code_arg, const NodeID arg2_arg) : @@ -1405,6 +1455,45 @@ BinaryOpNode::collectExogenous(set<pair<int, int> > &result) const arg2->collectExogenous(result); } +NodeID +BinaryOpNode::toStatic(DataTree &static_datatree) const +{ + NodeID sarg1 = arg1->toStatic(static_datatree); + NodeID sarg2 = arg2->toStatic(static_datatree); + switch(op_code) + { + case oPlus: + return static_datatree.AddPlus(sarg1, sarg2); + case oMinus: + return static_datatree.AddMinus(sarg1, sarg2); + case oTimes: + return static_datatree.AddTimes(sarg1, sarg2); + case oDivide: + return static_datatree.AddDivide(sarg1, sarg2); + case oPower: + return static_datatree.AddPower(sarg1, sarg2); + case oEqual: + return static_datatree.AddEqual(sarg1, sarg2); + case oMax: + return static_datatree.AddMaX(sarg1, sarg2); + case oMin: + return static_datatree.AddMin(sarg1, sarg2); + case oLess: + return static_datatree.AddLess(sarg1, sarg2); + case oGreater: + return static_datatree.AddGreater(sarg1, sarg2); + case oLessEqual: + return static_datatree.AddLessEqual(sarg1, sarg2); + case oGreaterEqual: + return static_datatree.AddGreaterEqual(sarg1, sarg2); + case oEqualEqual: + return static_datatree.AddEqualEqual(sarg1, sarg2); + case oDifferent: + return static_datatree.AddDifferent(sarg1, sarg2); + } +} + + TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) : ExprNode(datatree_arg), @@ -1697,6 +1786,19 @@ TrinaryOpNode::collectExogenous(set<pair<int, int> > &result) const arg3->collectExogenous(result); } +NodeID +TrinaryOpNode::toStatic(DataTree &static_datatree) const +{ + NodeID sarg1 = arg1->toStatic(static_datatree); + NodeID sarg2 = arg2->toStatic(static_datatree); + NodeID sarg3 = arg2->toStatic(static_datatree); + switch(op_code) + { + case oNormcdf: + return static_datatree.AddNormcdf(sarg1, sarg2, sarg3); + } +} + UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg, int symb_id_arg, @@ -1792,3 +1894,13 @@ UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutput cerr << "UnknownFunctionNode::compile: operation impossible!" << endl; exit(EXIT_FAILURE); } + +NodeID +UnknownFunctionNode::toStatic(DataTree &static_datatree) const +{ + vector<NodeID> static_arguments; + for(vector<NodeID>::const_iterator it = arguments.begin(); + it != arguments.end(); it++) + static_arguments.push_back((*it)->toStatic(static_datatree)); + return static_datatree.AddUnknownFunction(datatree.symbol_table.getName(symb_id), static_arguments); +} diff --git a/ExprNode.hh b/ExprNode.hh index e3134e874229a85af67c5a4eb6db7a4206fe98a0..ab7f3a06d30ac855ecdd2aa2fd0e344f3240e320 100644 --- a/ExprNode.hh +++ b/ExprNode.hh @@ -86,7 +86,7 @@ typedef map<int, double> eval_context_type; class ExprNode { friend class DataTree; - friend class ModelTree; + friend class DynamicModel; friend class ExprNodeLess; friend class NumConstNode; friend class VariableNode; @@ -159,6 +159,12 @@ public: virtual double eval(const eval_context_type &eval_context) const throw (EvalException) = 0; virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const = 0; + //! Creates a static version of this node + /*! + This method duplicates the current node by creating a similar node from which all leads/lags have been stripped, + adds the result in the static_datatree argument (and not in the original datatree), and returns it. + */ + virtual NodeID toStatic(DataTree &static_datatree) const = 0; }; //! Object used to compare two nodes (using their indexes) @@ -186,6 +192,7 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; + virtual NodeID toStatic(DataTree &static_datatree) const; }; //! Symbol or variable node @@ -214,12 +221,12 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; + virtual NodeID toStatic(DataTree &static_datatree) const; }; //! Unary operator node class UnaryOpNode : public ExprNode { - friend class DataTree; private: const NodeID arg; const UnaryOpcode op_code; @@ -243,12 +250,16 @@ public: static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; + //! Returns operand + NodeID get_arg() const { return(arg); }; + //! Returns op code + UnaryOpcode get_op_code() const { return(op_code); }; + virtual NodeID toStatic(DataTree &static_datatree) const; }; //! Binary operator node class BinaryOpNode : public ExprNode { - friend class ModelTree; private: const NodeID arg1, arg2; const BinaryOpcode op_code; @@ -273,15 +284,15 @@ public: static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; - virtual NodeID get_arg1() { return(arg1);}; - virtual NodeID get_arg2() { return(arg2);}; + //! Returns first operand + NodeID get_arg1() const { return(arg1); }; + //! Returns second operand + NodeID get_arg2() const { return(arg2); }; + //! Returns op code + BinaryOpcode get_op_code() const { return(op_code); }; + virtual NodeID toStatic(DataTree &static_datatree) const; }; -enum TrinaryOpcode - { - oNormcdf - }; - //! Trinary operator node class TrinaryOpNode : public ExprNode { @@ -310,6 +321,7 @@ public: static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; + virtual NodeID toStatic(DataTree &static_datatree) const; }; //! Unknown function node @@ -337,6 +349,7 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; + virtual NodeID toStatic(DataTree &static_datatree) const; }; //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model diff --git a/Makefile b/Makefile index 9abcaca511a3e9bf5563fa4d5451e9664de969d8..9c0153d286f38604626d44a35770323fcbea3395 100644 --- a/Makefile +++ b/Makefile @@ -38,6 +38,8 @@ MAIN_OBJS = \ DynareBison.o \ ComputingTasks.o \ ModelTree.o \ + StaticModel.o \ + DynamicModel.o \ NumericalConstants.o \ NumericalInitialization.o \ Shocks.o \ diff --git a/ModFile.cc b/ModFile.cc index 864cc3e4c0fc4396a7fa6fc830187d8d6715f155..b3ae75732c16c7b591cc744cced6351e9af843ae 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -24,7 +24,8 @@ #include "ModFile.hh" ModFile::ModFile() : expressions_tree(symbol_table, num_constants), - model_tree(symbol_table, num_constants), + static_model(symbol_table, num_constants), + dynamic_model(symbol_table, num_constants), linear(false) { } @@ -58,7 +59,7 @@ ModFile::evalAllExpressions() } // Evaluate model local variables - model_tree.fillEvalContext(global_eval_context); + dynamic_model.fillEvalContext(global_eval_context); cout << "done" << endl; @@ -100,7 +101,7 @@ ModFile::checkPass() || mod_file_struct.ramsey_policy_present; // Allow empty model only when doing a standalone BVAR estimation - if (model_tree.equation_number() == 0 + if (dynamic_model.equation_number() == 0 && (mod_file_struct.check_present || mod_file_struct.simul_present || stochastic_statement_present)) @@ -109,7 +110,7 @@ ModFile::checkPass() exit(EXIT_FAILURE); } - if (mod_file_struct.simul_present && stochastic_statement_present && model_tree.mode==0) + if (mod_file_struct.simul_present && stochastic_statement_present && dynamic_model.mode == 0) { cerr << "ERROR: A .mod file cannot contain both a simul command and one of {stoch_simul, estimation, forecast, osr, ramsey_policy}" << endl; exit(EXIT_FAILURE); @@ -125,23 +126,29 @@ ModFile::checkPass() */ if (!mod_file_struct.ramsey_policy_present && !((mod_file_struct.bvar_density_present || mod_file_struct.bvar_forecast_present) - && model_tree.equation_number() == 0) - && (model_tree.equation_number() != symbol_table.endo_nbr())) + && dynamic_model.equation_number() == 0) + && (dynamic_model.equation_number() != symbol_table.endo_nbr())) { - cerr << "ERROR: There are " << model_tree.equation_number() << " equations but " << symbol_table.endo_nbr() << " endogenous variables!" << endl; + cerr << "ERROR: There are " << dynamic_model.equation_number() << " equations but " << symbol_table.endo_nbr() << " endogenous variables!" << endl; exit(EXIT_FAILURE); } + + cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl; } void ModFile::computingPass(bool no_tmp_terms) { // Mod file may have no equation (for example in a standalone BVAR estimation) - if (model_tree.equation_number() > 0) + if (dynamic_model.equation_number() > 0) { - // Set things to compute + // Compute static model and its derivatives + dynamic_model.toStatic(static_model); + static_model.computingPass(no_tmp_terms); + + // Set things to compute for dynamic model if (mod_file_struct.simul_present) - model_tree.computeJacobian = true; + dynamic_model.computeJacobian = true; else { if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3) @@ -149,13 +156,13 @@ ModFile::computingPass(bool no_tmp_terms) cerr << "Incorrect order option..." << endl; exit(EXIT_FAILURE); } - model_tree.computeJacobianExo = true; + dynamic_model.computeJacobianExo = true; if (mod_file_struct.order_option >= 2) - model_tree.computeHessian = true; + dynamic_model.computeHessian = true; if (mod_file_struct.order_option == 3) - model_tree.computeThirdDerivatives = true; + dynamic_model.computeThirdDerivatives = true; } - model_tree.computingPass(global_eval_context, no_tmp_terms); + dynamic_model.computingPass(global_eval_context, no_tmp_terms); } for(vector<Statement *>::iterator it = statements.begin(); it != statements.end(); it++) @@ -209,18 +216,17 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const << "warning warning_old_state" << endl; mOutputFile << "logname_ = '" << basename << ".log';" << endl; mOutputFile << "diary " << basename << ".log" << endl; - mOutputFile << "options_.model_mode = " << model_tree.mode << ";\n"; - if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode) + mOutputFile << "options_.model_mode = " << dynamic_model.mode << ";\n"; + if (dynamic_model.mode == eSparseMode || dynamic_model.mode == eSparseDLLMode) { mOutputFile << "addpath " << basename << ";\n"; - mOutputFile << "delete('" << basename << "_static.m');\n"; mOutputFile << "delete('" << basename << "_dynamic.m');\n"; } - if (model_tree.equation_number() > 0) + if (dynamic_model.equation_number() > 0) { - if (model_tree.mode == eDLLMode || model_tree.mode == eSparseDLLMode) + if (dynamic_model.mode == eDLLMode || dynamic_model.mode == eSparseDLLMode) { mOutputFile << "if exist('" << basename << "_static.c')" << endl; mOutputFile << " clear " << basename << "_static" << endl; @@ -244,11 +250,11 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const if (linear == 1) mOutputFile << "options_.linear = 1;" << endl; - if (model_tree.equation_number() > 0) + if (dynamic_model.equation_number() > 0) { - model_tree.writeOutput(mOutputFile); - model_tree.writeStaticFile(basename); - model_tree.writeDynamicFile(basename); + dynamic_model.writeOutput(mOutputFile); + static_model.writeStaticFile(basename); + dynamic_model.writeDynamicFile(basename); } // Print statements @@ -259,7 +265,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl; mOutputFile << "diary off" << endl; - if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode) + if (dynamic_model.mode == eSparseMode || dynamic_model.mode == eSparseDLLMode) mOutputFile << "rmpath " << basename << ";\n"; mOutputFile << endl << "disp(['Total computing time : ' dynsec2hms(toc) ]);" << endl; diff --git a/ModFile.hh b/ModFile.hh index 1ea1a26b54e738500d245029b2df3349ca71af20..8ee32c1271345d975613ec7461b6b368f8a40223 100644 --- a/ModFile.hh +++ b/ModFile.hh @@ -28,7 +28,8 @@ using namespace std; #include "SymbolTable.hh" #include "NumericalConstants.hh" #include "NumericalInitialization.hh" -#include "ModelTree.hh" +#include "StaticModel.hh" +#include "DynamicModel.hh" #include "VariableTable.hh" #include "Statement.hh" @@ -44,8 +45,10 @@ public: NumericalConstants num_constants; //! Expressions outside model block DataTree expressions_tree; - //! Model equations and their derivatives - ModelTree model_tree; + //! Static model + StaticModel static_model; + //! Dynamic model + DynamicModel dynamic_model; //! Option linear bool linear; //! Global evaluation context @@ -69,7 +72,7 @@ public: //! Execute computations /*! \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */ void computingPass(bool no_tmp_terms); - //! Writes Matlab/Scilab output files + //! Writes Matlab/Octave output files /*! \param basename The base name used for writing output files. Should be the name of the mod file without its extension \param clear_all Should a "clear all" instruction be written to output ? diff --git a/ModelTree.cc b/ModelTree.cc index 86aeec4fc3ab0a8fc75bb44c7ade2ba84797fa57..830199e5c97ed703d77209bbc384c236ee998581 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -19,37 +19,13 @@ #include <cstdlib> #include <iostream> -#include <fstream> -#include <sstream> -#include <cstring> -#include <cmath> - -// For mkdir() and chdir() -#ifdef _WIN32 -# include <direct.h> -#else -# include <unistd.h> -# include <sys/stat.h> -# include <sys/types.h> -#endif #include "ModelTree.hh" -#include "ModelGraph.hh" - ModelTree::ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg) : DataTree(symbol_table_arg, num_constants_arg), - mode(eStandardMode), - cutoff(1e-15), - markowitz(0.7), - new_SGE(true), - computeJacobian(false), - computeJacobianExo(false), - computeHessian(false), - computeStaticHessian(false), - computeThirdDerivatives(false), - block_triangular(symbol_table_arg) + mode(eStandardMode) { } @@ -71,23 +47,10 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag, output << 0; } -void -ModelTree::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 ModelTree::derive(int order) { - cout << "Processing derivation ..." << endl; - cout << " Processing Order 1... "; for (int var = 0; var < variable_table.size(); var++) for (int eq = 0; eq < (int) equations.size(); eq++) { @@ -96,11 +59,9 @@ ModelTree::derive(int order) continue; first_derivatives[make_pair(eq, var)] = d1; } - cout << "done" << endl; if (order >= 2) { - cout << " Processing Order 2... "; for (first_derivatives_type::const_iterator it = first_derivatives.begin(); it != first_derivatives.end(); it++) { @@ -117,12 +78,10 @@ ModelTree::derive(int order) second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2; } } - cout << "done" << endl; } if (order >= 3) { - cout << " Processing Order 3... "; for (second_derivatives_type::const_iterator it = second_derivatives.begin(); it != second_derivatives.end(); it++) { @@ -143,7 +102,6 @@ ModelTree::derive(int order) third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3; } } - cout << "done" << endl; } } @@ -223,35 +181,6 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t } } - -void -ModelTree::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->arg1; - Id->collectEndogenous(endogenous); - Id = eq_node->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->arg1; - Id->collectExogenous(exogenous); - Id = eq_node->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 ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const { @@ -259,12 +188,12 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) { BinaryOpNode *eq_node = equations[eq]; - NodeID lhs = eq_node->arg1; + NodeID lhs = eq_node->get_arg1(); output << "lhs ="; lhs->writeOutput(output, output_type, temporary_terms); output << ";" << endl; - NodeID rhs = eq_node->arg2; + NodeID rhs = eq_node->get_arg2(); output << "rhs ="; rhs->writeOutput(output, output_type, temporary_terms); output << ";" << endl; @@ -274,2873 +203,19 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) } void -ModelTree::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 -ModelTree::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->arg1; - rhs = eq_node->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+1 - << ", equation=" << eq+1 << endl; - } - } - //jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr; - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; - output << " g1_x(" << eqr+1 << ", " - << varr+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr() << ") = "; - writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - } - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=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++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i]; - output << " g1_o(" << eqr+1 << ", " - << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << 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: - 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 - << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - } - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; - output << " g1_x(" << eqr+1 << ", " << varr+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = "; - writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - } - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=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++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i]; - output << " g1_o(" << eqr+1 << ", " - << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - } - } - output << " varargout{1}=g1_x;\n"; - output << " varargout{2}=g1_o;\n"; - output << " else" << endl; - - 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 << ") = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - output << " end;\n"; - break; - case SOLVE_TWO_BOUNDARIES_SIMPLE: - case SOLVE_TWO_BOUNDARIES_COMPLETE: - output << " if ~jacobian_eval" << endl; - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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]; - if (k==0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")"; - else if (k==1) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_)*y(it_+1, " << var+1 << ")"; - else if (k>0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")"; - else if (k<0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; - if (k==0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; - else if (k==1) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; - else if (k>0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; - else if (k<0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << var+1 - << ", equation=" << 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 (i = 0;i < ModelBlock->Block_List[j].Size;i++) - { - output << " " << Uf[ModelBlock->Block_List[j].Equation[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[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].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[j].Size;i++) - output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; -#endif - - output << " else" << endl; - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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-ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lag_Endo)*ModelBlock->Block_List[j].Size << ") = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - } - jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size; - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; - output << " g1_x(" << eqr+1 << ", " - << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = "; - writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable (exogenous)=" << symbol_table.getName(var) - << "(" << k << ") " << var+1 << " " << varr+1 - << ", equation=" << eq+1 << endl; - } - } - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=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++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i]; - output << " g1_o(" << eqr+1 << ", " - << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(" << k << ") " << 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; - } - prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; - output.close(); - } -} - -void -ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const string &static_basename) const -{ - int i,j,k,m, var, eq, g1_index = 1; - string tmp_s, sps; - ostringstream tmp_output, tmp1_output, global_output; - NodeID lhs=NULL, rhs=NULL; - BinaryOpNode *eq_node; - map<NodeID, int> reference_count; - int prev_Simulation_Type=-1; - int nze=0; - bool *IM, *IMl; - ofstream output; - temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); - //---------------------------------------------------------------------- - 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 - tmp1_output.str(""); - tmp1_output << static_basename << "_" << j+1 << ".m"; - output.open(tmp1_output.str().c_str(), ios::out | ios::binary); - output << "%\n"; - output << "% " << tmp1_output.str() << " : Computes static 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] = " << static_basename << "_" << j+1 << "(y, x, params, jacobian_eval)\n"; - else - output << "function [residual, g1, g2, g3] = " << static_basename << "_" << j+1 << "(y, x, params, jacobian_eval)\n"; - output << " % ////////////////////////////////////////////////////////////////////////" << endl - << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " - << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl - << " % // Simulation type "; - output << BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //" << endl - << " % ////////////////////////////////////////////////////////////////////////" << endl; - //The Temporary terms - //output << global_output.str(); - 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"; - } - - int n=ModelBlock->Block_List[j].Size; - int n1=symbol_table.endo_nbr(); - IM=(bool*)malloc(n*n*sizeof(bool)); - memset(IM, 0, n*n*sizeof(bool)); - for (m=-ModelBlock->Block_List[j].Max_Lag;m<=ModelBlock->Block_List[j].Max_Lead;m++) - { - IMl=block_triangular.incidencematrix.Get_IM(m, eEndogenous); - if (IMl) - { - for (i=0;i<n;i++) - { - eq=ModelBlock->Block_List[j].Equation[i]; - for (k=0;k<n;k++) - { - var=ModelBlock->Block_List[j].Variable[k]; - IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var]; - } - } - } - } - for (nze=0, i=0;i<n*n;i++) - { - nze+=IM[i]; - } - memset(IM, 0, n*n*sizeof(bool)); - 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 << " g1=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n"; - output << " residual=zeros(" << ModelBlock->Block_List[j].Size << ",1);\n"; - } - 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, oMatlabStaticModelSparse, temporary_terms); - output << " = "; - (*it)->writeOutput(output, oMatlabStaticModelSparse, tt2); - // Insert current node into tt2 - tt2.insert(*it); - output << ";" << endl; - } - //cout << "variable_table.getSymbolID(variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]))=" << variable_table.getSymbolID(variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i])) << "\n"; - string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i])); - output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " - << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; - eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; - lhs = eq_node->arg1; - rhs = eq_node->arg2; - tmp_output.str(""); - lhs->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms); - output << " "; - switch (ModelBlock->Block_List[j].Simulation_Type) - { - case EVALUATE_BACKWARD: - case EVALUATE_FORWARD: - output << tmp_output.str(); - output << " = "; - rhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); - output << ";\n"; - break; - case EVALUATE_BACKWARD_R: - case EVALUATE_FORWARD_R: - rhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); - output << " = "; - lhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); - output << ";\n"; - break; - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FORWARD_SIMPLE: - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_COMPLETE: - case SOLVE_TWO_BOUNDARIES_COMPLETE: - case SOLVE_TWO_BOUNDARIES_SIMPLE: - goto end; - default: - end: - output << sps << "residual(" << i+1 << ") = ("; - output << tmp_output.str(); - output << ") - ("; - rhs->writeOutput(output, oMatlabStaticModelSparse, 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 - output << " " << sps << "% Jacobian " << endl; - switch (ModelBlock->Block_List[j].Simulation_Type) - { - case EVALUATE_BACKWARD: - case EVALUATE_FORWARD: - case EVALUATE_BACKWARD_R: - case EVALUATE_FORWARD_R: - output << " if(jacobian_eval)\n"; - output << " g1( " << g1_index << ", " << g1_index << ")="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[0]), 0, oMatlabStaticModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[0])) - << "(" << variable_table.getLag(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[0])) - << ") " << ModelBlock->Block_List[j].Variable[0]+1 - << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; - output << " end\n"; - break; - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FORWARD_SIMPLE: - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_COMPLETE: - case SOLVE_TWO_BOUNDARIES_COMPLETE: - case SOLVE_TWO_BOUNDARIES_SIMPLE: - output << " g2=0;g3=0;\n"; - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabStaticModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) - << "(" << k << ") " << var+1 - << ", equation=" << 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 - } - } -#ifdef CONDITION - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].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[j].Size;i++) - output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; -#endif - break; - default: - break; - } - prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; - free(IM); - output << "return;\n"; - output.close(); - } - //output << "return;\n\n\n"; -} - - -void -ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const -{ - struct Uff_l - { - int u, var, lag; - Uff_l *pNext; - }; - - struct Uff - { - Uff_l *Ufl, *Ufl_First; - int eqr; - }; - - int i,j,k,m, v, ModelBlock_Aggregated_Count, k0, k1; - string tmp_s; - ostringstream tmp_output; - ofstream code_file; - NodeID lhs=NULL, rhs=NULL; - BinaryOpNode *eq_node; - bool lhs_rhs_done; - Uff Uf[symbol_table.endo_nbr()]; - map<NodeID, int> reference_count; - map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number; - int prev_Simulation_Type=-1; - //SymbolicGaussElimination SGE; - bool file_open=false; - temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); - //---------------------------------------------------------------------- - 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 - code_file.write(&FDIMT, sizeof(FDIMT)); - k=temporary_terms.size(); - code_file.write(reinterpret_cast<char *>(&k),sizeof(k)); - //search for successive and identical blocks - i=k=k0=0; - ModelBlock_Aggregated_Count=-1; - for (j = 0;j < ModelBlock->Size;j++) - { - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) - && (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 )) - { - } - else - { - k=k0=0; - ModelBlock_Aggregated_Count++; - } - k0+=ModelBlock->Block_List[j].Size; - ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0; - ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k; - prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; - } - ModelBlock_Aggregated_Count++; - //For each block - j=0; - for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++) - { - k1=j; - if (k0>0) - code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); - code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK)); - v=ModelBlock_Aggregated_Number[k0]; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - v=ModelBlock->Block_List[j].Simulation_Type; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - for (k=0; k<ModelBlock_Aggregated_Size[k0]; k++) - { - for (i=0; i < ModelBlock->Block_List[j].Size;i++) - { - code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i])); - code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i])); - code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i])); - } - j++; - } - j=k1; - if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || - ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE) - { - code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].is_linear),sizeof(ModelBlock->Block_List[j].is_linear)); - v=block_triangular.ModelBlock->Block_List[j].IM_lead_lag[block_triangular.ModelBlock->Block_List[j].Max_Lag + block_triangular.ModelBlock->Block_List[j].Max_Lead].u_finish + 1; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - v=symbol_table.endo_nbr(); - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - v=block_triangular.ModelBlock->Block_List[j].Max_Lag; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - v=block_triangular.ModelBlock->Block_List[j].Max_Lead; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) - //{ - int u_count_int=0; - Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open, - ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE); - v=u_count_int; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - file_open=true; - //} - } - for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++) - { - //For a block composed of a single equation determines whether we have to evaluate or to solve the equation - if (ModelBlock->Block_List[j].Size==1) - { - lhs_rhs_done=true; - eq_node = equations[ModelBlock->Block_List[j].Equation[0]]; - lhs = eq_node->arg1; - rhs = eq_node->arg2; - } - else - lhs_rhs_done=false; - // The equations - for (i = 0;i < ModelBlock->Block_List[j].Size;i++) - { - //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); - //The Temporary terms - temporary_terms_type tt2; -#ifdef DEBUGC - k=0; -#endif - 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++) - { - (*it)->compile(code_file,false, output_type, tt2, map_idx); - code_file.write(&FSTPT, sizeof(FSTPT)); - map_idx_type::const_iterator ii=map_idx.find((*it)->idx); - v=(int)ii->second; - code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); - // 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 = ModelBlock->Block_List[j].Temporary_terms->begin(); - it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) - { - map_idx_type::const_iterator ii=map_idx.find((*it)->idx); - cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n"; - } -#endif - if (!lhs_rhs_done) - { - eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; - lhs = eq_node->arg1; - rhs = eq_node->arg2; - } - switch (ModelBlock->Block_List[j].Simulation_Type) - { - case EVALUATE_BACKWARD: - case EVALUATE_FORWARD: - rhs->compile(code_file,false, output_type, temporary_terms, map_idx); - lhs->compile(code_file,true, output_type, temporary_terms, map_idx); - break; - case EVALUATE_BACKWARD_R: - case EVALUATE_FORWARD_R: - lhs->compile(code_file,false, output_type, temporary_terms, map_idx); - rhs->compile(code_file,true, output_type, temporary_terms, map_idx); - break; - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_COMPLETE: - v=ModelBlock->Block_List[j].Equation[i]; - Uf[v].eqr=i; - Uf[v].Ufl=NULL; - goto end; - case SOLVE_TWO_BOUNDARIES_COMPLETE: - case SOLVE_TWO_BOUNDARIES_SIMPLE: - v=ModelBlock->Block_List[j].Equation[i]; - Uf[v].eqr=i; - Uf[v].Ufl=NULL; - goto end; - default: - end: - lhs->compile(code_file,false, output_type, temporary_terms, map_idx); - rhs->compile(code_file,false, output_type, temporary_terms, map_idx); - code_file.write(&FBINARY, sizeof(FBINARY)); - int v=oMinus; - code_file.write(reinterpret_cast<char *>(&v),sizeof(v)); - code_file.write(&FSTPR, sizeof(FSTPR)); - code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); -#ifdef CONDITION - if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE) - output << " condition[" << i << "]=0;\n"; -#endif - } - } - code_file.write(&FENDEQU, sizeof(FENDEQU)); - // The Jacobian if we have to solve the block - 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) - { - switch (ModelBlock->Block_List[j].Simulation_Type) - { - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FORWARD_SIMPLE: - compileDerivative(code_file, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, output_type, map_idx); - code_file.write(&FSTPG, sizeof(FSTPG)); - v=0; - code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); - break; - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_COMPLETE: - 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - int v=ModelBlock->Block_List[j].Equation[eqr]; - if (!Uf[v].Ufl) - { - Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); - Uf[v].Ufl_First=Uf[v].Ufl; - } - else - { - Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); - Uf[v].Ufl=Uf[v].Ufl->pNext; - } - Uf[v].Ufl->pNext=NULL; - Uf[v].Ufl->u=u; - Uf[v].Ufl->var=var; - compileDerivative(code_file, eq, var, 0, output_type, map_idx); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast<char *>(&u), sizeof(u)); - } - for (i = 0;i < ModelBlock->Block_List[j].Size;i++) - { - code_file.write(&FLDR, sizeof(FLDR)); - code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); - code_file.write(&FLDZ, sizeof(FLDZ)); - int v=ModelBlock->Block_List[j].Equation[i]; - for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext) - { - code_file.write(&FLDU, sizeof(FLDU)); - code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); - code_file.write(&FLDV, sizeof(FLDV)); - char vc=eEndogenous; - code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc)); - code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var)); - int v1=0; - code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); - code_file.write(&FBINARY, sizeof(FBINARY)); - v1=oTimes; - code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); - code_file.write(&FCUML, sizeof(FCUML)); - } - 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; - } - code_file.write(&FBINARY, sizeof(FBINARY)); - v=oMinus; - code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); - } - break; - case SOLVE_TWO_BOUNDARIES_COMPLETE: - case SOLVE_TWO_BOUNDARIES_SIMPLE: - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - int v=ModelBlock->Block_List[j].Equation[eqr]; - if (!Uf[v].Ufl) - { - Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); - Uf[v].Ufl_First=Uf[v].Ufl; - } - else - { - Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); - Uf[v].Ufl=Uf[v].Ufl->pNext; - } - Uf[v].Ufl->pNext=NULL; - Uf[v].Ufl->u=u; - Uf[v].Ufl->var=var; - Uf[v].Ufl->lag=k; - compileDerivative(code_file, eq, var, k, output_type, map_idx); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast<char *>(&u), sizeof(u)); -#ifdef CONDITION - output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n"; - output << " condition[" << eqr << "]=u[" << u << "+Per_u_];\n"; -#endif - } - } - for (i = 0;i < ModelBlock->Block_List[j].Size;i++) - { - code_file.write(&FLDR, sizeof(FLDR)); - code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); - code_file.write(&FLDZ, sizeof(FLDZ)); - int v=ModelBlock->Block_List[j].Equation[i]; - for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext) - { - code_file.write(&FLDU, sizeof(FLDU)); - code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); - code_file.write(&FLDV, sizeof(FLDV)); - char vc=eEndogenous; - code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc)); - int v1=Uf[v].Ufl->var; - code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); - v1=Uf[v].Ufl->lag; - code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); - code_file.write(&FBINARY, sizeof(FBINARY)); - v1=oTimes; - code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1)); - code_file.write(&FCUML, sizeof(FCUML)); - } - 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; - } - code_file.write(&FBINARY, sizeof(FBINARY)); - v=oMinus; - code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); -#ifdef CONDITION - output << " if (fabs(condition[" << i << "])<fabs(u[" << i << "+Per_u_]))\n"; - output << " condition[" << i << "]=u[" << i << "+Per_u_];\n"; -#endif - } -#ifdef CONDITION - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - output << " u[" << u << "+Per_u_] /= condition[" << eqr << "];\n"; - } - } - for (i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n"; -#endif - break; - default: - break; - } - - prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; - } - j++; - } - } - code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); - code_file.write(&FEND, sizeof(FEND)); - code_file.close(); -} - - -void -ModelTree::writeStaticMFile(const string &static_basename) const -{ - string filename = static_basename + ".m"; - - ofstream mStaticModelFile; - mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); - if (!mStaticModelFile.is_open()) - { - cerr << "Error: Can't open file " << filename << " for writing" << endl; - exit(EXIT_FAILURE); - } - // Writing comments and function definition command - mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "(y, x, params)" << endl - << "%" << endl - << "% Status : Computes static model for Dynare" << endl - << "%" << endl - << "% Warning : this file is generated automatically by Dynare" << endl - << "% from model file (.mod)" << endl << endl; - - writeStaticModel(mStaticModelFile); - - mStaticModelFile.close(); -} - - -void -ModelTree::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; - - writeDynamicModel(mDynamicModelFile); - - mDynamicModelFile.close(); -} - -void -ModelTree::writeStaticCFile(const string &static_basename) const -{ - string filename = static_basename + ".c"; - - ofstream mStaticModelFile; - mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); - if (!mStaticModelFile.is_open()) - { - cerr << "Error: Can't open file " << filename << " for writing" << endl; - exit(EXIT_FAILURE); - } - mStaticModelFile << "/*" << endl - << " * " << filename << " : Computes static model for Dynare" << 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; - - // Writing the function Static - writeStaticModel(mStaticModelFile); - - // Writing the gateway routine - mStaticModelFile << "/* 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;" << 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 - << " 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() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix g1. */" << endl - << " g1 = mxGetPr(plhs[1]);" << endl - << " }" << endl - << endl - << " /* Call the C Static. */" << endl - << " Static(y, x, params, residual, g1);" << endl - << "}" << endl; - - mStaticModelFile.close(); -} - -void -ModelTree::writeDynamicCFile(const string &dynamic_basename) const +ModelTree::addEquation(NodeID eq) { - string filename = dynamic_basename + ".c"; - ofstream mDynamicModelFile; + BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq); - mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary); - if (!mDynamicModelFile.is_open()) + if (beq == NULL || beq->get_op_code() != oEqual) { - 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; - - // Writing the function body - writeDynamicModel(mDynamicModelFile); - - // 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, *g2;" << 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() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix g1. */" << endl - << " g1 = mxGetPr(plhs[1]);" << endl - << " }" << endl - << endl - << " g2 = NULL;" << endl - << " if (nlhs >= 3)" << endl - << " {" << endl - << " /* Set the output pointer to the output matrix g2. */" << endl; - int g2_ncols = variable_table.getDynJacobianColsNbr(computeJacobianExo)*variable_table.getDynJacobianColsNbr(computeJacobianExo); - mDynamicModelFile << " plhs[2] = mxCreateSparse(" << equations.size() << ", " << g2_ncols << ", " - << 5*g2_ncols << ", mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix g1. */" << endl - << " g2 = mxGetPr(plhs[2]);" << endl - << " }" << endl - << endl - << " /* Call the C subroutines. */" << endl - << " Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);" << endl - << "}" << endl; - mDynamicModelFile.close(); -} - -void -ModelTree::writeStaticModel(ostream &StaticOutput) const -{ - ostringstream model_output; // Used for storing model equations - ostringstream jacobian_output; // Used for storing jacobian equations - ostringstream hessian_output; - ostringstream lsymetric; // For symmetric elements in hessian - - ExprNodeOutputType output_type = (mode == eDLLMode ? oCStaticModel : oMatlabStaticModel); - - writeModelLocalVariables(model_output, output_type); - - writeTemporaryTerms(model_output, output_type); - - writeModelEquations(model_output, output_type); - - // Write Jacobian w.r. to endogenous only - for (first_derivatives_type::const_iterator it = first_derivatives.begin(); - it != first_derivatives.end(); it++) - { - int eq = it->first.first; - int var = it->first.second; - NodeID d1 = it->second; - - if (variable_table.getType(var) == eEndogenous) - { - ostringstream g1; - g1 << " g1"; - matrixHelper(g1, eq, symbol_table.getTypeSpecificID(variable_table.getSymbolID(var)), output_type); - - jacobian_output << g1.str() << "=" << g1.str() << "+"; - d1->writeOutput(jacobian_output, output_type, temporary_terms); - jacobian_output << ";" << endl; - } - } - - // Write Hessian w.r. to endogenous only - if (computeStaticHessian) - for (second_derivatives_type::const_iterator it = second_derivatives.begin(); - it != second_derivatives.end(); it++) - { - int eq = it->first.first; - int var1 = it->first.second.first; - int var2 = it->first.second.second; - NodeID d2 = it->second; - - // Keep only derivatives w.r. to endogenous variables - if (variable_table.getType(var1) == eEndogenous - && variable_table.getType(var2) == eEndogenous) - { - int id1 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var1)); - int id2 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var2)); - - int col_nb = id1*symbol_table.endo_nbr()+id2; - int col_nb_sym = id2*symbol_table.endo_nbr()+id1; - - hessian_output << " g2"; - matrixHelper(hessian_output, eq, col_nb, output_type); - hessian_output << " = "; - d2->writeOutput(hessian_output, output_type, temporary_terms); - hessian_output << ";" << endl; - - // Treating symetric elements - if (var1 != var2) - { - lsymetric << " g2"; - matrixHelper(lsymetric, eq, col_nb_sym, output_type); - lsymetric << " = " << "g2"; - matrixHelper(lsymetric, eq, col_nb, output_type); - lsymetric << ";" << endl; - } - } - } - - // Writing ouputs - if (mode != eDLLMode) - { - StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl - << "%" << endl - << "% Model equations" << endl - << "%" << endl - << endl - << model_output.str() - << "if ~isreal(residual)" << endl - << " residual = real(residual)+imag(residual).^2;" << endl - << "end" << endl - << "if nargout >= 2," << endl - << " g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl - << endl - << "%" << endl - << "% Jacobian matrix" << endl - << "%" << endl - << endl - << jacobian_output.str() - << " if ~isreal(g1)" << endl - << " g1 = real(g1)+2*imag(g1);" << endl - << " end" << endl - << "end" << endl; - if (computeStaticHessian) - { - StaticOutput << "if nargout >= 3,\n"; - // Writing initialization instruction for matrix g2 - int ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr(); - StaticOutput << " g2 = sparse([],[],[], " << equations.size() << ", " << ncols << ", " << 5*ncols << ");" << endl - << endl - << "%" << endl - << "% Hessian matrix" << endl - << "%" << endl - << endl - << hessian_output.str() - << lsymetric.str() - << "end;" << endl; - } - } - else - { - StaticOutput << "void Static(double *y, double *x, double *params, double *residual, double *g1)" << endl - << "{" << endl - << " double lhs, rhs;" << endl - // Writing residual equations - << " /* Residual equations */" << endl - << " if (residual == NULL)" << endl - << " return;" << endl - << " else" << endl - << " {" << endl - << model_output.str() - // Writing Jacobian - << " /* Jacobian for endogenous variables without lag */" << endl - << " if (g1 == NULL)" << endl - << " return;" << endl - << " else" << endl - << " {" << endl - << jacobian_output.str() - << " }" << endl - << " }" << endl - << "}" << endl << endl; - } -} - -string -ModelTree::reform(const string name1) const -{ - string name=name1; - int pos = name.find("\\", 0); - while (pos >= 0) - { - if (name.substr(pos + 1, 1) != "\\") - { - name = name.insert(pos, "\\"); - pos++; - } - pos++; - pos = name.find("\\", pos); - } - return (name); -} - -void -ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num, - int &u_count_int, bool &file_open, bool is_two_boundaries) const -{ - int j; - std::ofstream SaveCode; - if (file_open) - SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate ); - else - SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::binary); - if (!SaveCode.is_open()) - { - cout << "Error : Can't open file \"" << bin_basename << ".bin\" for writing\n"; - exit(EXIT_FAILURE); - } - u_count_int=0; - for (int m=0;m<=block_triangular.ModelBlock->Block_List[num].Max_Lead+block_triangular.ModelBlock->Block_List[num].Max_Lag;m++) - { - int k1=m-block_triangular.ModelBlock->Block_List[num].Max_Lag; - for (j=0;j<block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].size;j++) - { - int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size; - int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j]; - int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j]; - SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); - SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr)); - SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1)); - SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u)); - u_count_int++; - } - } - if(is_two_boundaries) - { - for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++) - { - int eqr1=j; - int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods - +block_triangular.incidencematrix.Model_Max_Lead_Endo); - int k1=0; - SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); - SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr)); - SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1)); - SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); - u_count_int++; - } - } - //cout << "u_count_int=" << u_count_int << "\n"; - for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++) - { - int varr=block_triangular.ModelBlock->Block_List[num].Variable[j]; - SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr)); - } - for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++) - { - int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j]; - SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1)); - } - SaveCode.close(); -} - -void -ModelTree::writeSparseStaticMFile(const string &static_basename, const string &basename, const int mode) const -{ - string filename; - ofstream mStaticModelFile; - ostringstream tmp, tmp1, tmp_eq; - int i, k, prev_Simulation_Type, ga_index = 1; - bool skip_head, open_par=false; - - chdir(basename.c_str()); - filename = static_basename + ".m"; - mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); - if (!mStaticModelFile.is_open()) - { - cerr << "Error: Can't open file " << filename << " for writing" << endl; - exit(EXIT_FAILURE); - } - mStaticModelFile << "%\n"; - mStaticModelFile << "% " << filename << " : Computes static model for Dynare\n"; - mStaticModelFile << "%\n"; - mStaticModelFile << "% Warning : this file is generated automatically by Dynare\n"; - mStaticModelFile << "% from model file (.mod)\n\n"; - mStaticModelFile << "%/\n"; - mStaticModelFile << "function [varargout] = " << static_basename << "(varargin)\n"; - mStaticModelFile << " global oo_ M_ options_ ys0_ ;\n"; - bool OK=true; - ostringstream tmp_output; - for (temporary_terms_type::const_iterator it = temporary_terms.begin(); - it != temporary_terms.end(); it++) - { - if (OK) - OK=false; - else - tmp_output << " "; - (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); - } - if (tmp_output.str().length()>0) - mStaticModelFile << " global " << tmp_output.str() << " M_ ;\n"; - mStaticModelFile << " T_init=0;\n"; - tmp_output.str(""); - for (temporary_terms_type::const_iterator it = temporary_terms.begin(); - it != temporary_terms.end(); it++) - { - tmp_output << " "; - (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); - tmp_output << "=T_init;\n"; - } - if (tmp_output.str().length()>0) - mStaticModelFile << tmp_output.str(); - - mStaticModelFile << " y_kmin=M_.maximum_lag;\n"; - mStaticModelFile << " y_kmax=M_.maximum_lead;\n"; - mStaticModelFile << " y_size=M_.endo_nbr;\n"; - - - mStaticModelFile << " if(length(varargin)>0)\n"; - mStaticModelFile << " %A simple evaluation of the static model\n"; - mStaticModelFile << " y=varargin{1}(:);\n"; - mStaticModelFile << " ys=y;\n"; - mStaticModelFile << " g1=[];\n"; - mStaticModelFile << " x=varargin{2}(:);\n"; - mStaticModelFile << " params=varargin{3}(:);\n"; - mStaticModelFile << " residual=zeros(1, " << symbol_table.endo_nbr() << ");\n"; - prev_Simulation_Type=-1; - tmp.str(""); - tmp_eq.str(""); - for (i=0;i<block_triangular.ModelBlock->Size;i++) - { - k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) && - ((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R) - || (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))) - { - mStaticModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; - tmp_eq.str(""); - mStaticModelFile << " y_index=[" << tmp.str() << "];\n"; - tmp.str(""); - mStaticModelFile << tmp1.str(); - tmp1.str(""); - } - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; - } - if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R) - { - if (i==block_triangular.ModelBlock->Size-1) - { - mStaticModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; - tmp_eq.str(""); - mStaticModelFile << " y_index=[" << tmp.str() << "];\n"; - tmp.str(""); - mStaticModelFile << tmp1.str(); - tmp1.str(""); - } - } - - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && - (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) - skip_head=true; - else - skip_head=false; - - switch (k) - { - case EVALUATE_FORWARD: - case EVALUATE_BACKWARD: - case EVALUATE_FORWARD_R: - case EVALUATE_BACKWARD_R: - if (!skip_head) - { - ga_index = 1; - tmp1 << " [y, ga]=" << static_basename << "_" << i + 1 << "(y, x, params, 1);\n"; - tmp1 << " residual(y_index)=ys(y_index)-y(y_index);\n"; - tmp1 << " g1(y_index_eq, y_index) = ga;\n"; - } - else - ga_index++; - break; - case SOLVE_FORWARD_COMPLETE: - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_SIMPLE: - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_TWO_BOUNDARIES_COMPLETE: - case SOLVE_TWO_BOUNDARIES_SIMPLE: - mStaticModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - mStaticModelFile << " y_index = ["; - mStaticModelFile << tmp.str(); - mStaticModelFile << " ];\n"; - tmp.str(""); - tmp_eq.str(""); - mStaticModelFile << " [r, ga]=" << static_basename << "_" << i + 1 << "(y, x, params, 1);\n"; - mStaticModelFile << " g1(y_index_eq, y_index) = ga;\n"; - mStaticModelFile << " residual(y_index)=r;\n"; - break; - } - prev_Simulation_Type=k; - } - mStaticModelFile << " varargout{1}=residual';\n"; - mStaticModelFile << " varargout{2}=g1;\n"; - mStaticModelFile << " return;\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " %The deterministic simulation of the block decomposed static model\n"; - mStaticModelFile << " periods=options_.periods;\n"; - mStaticModelFile << " maxit_=options_.maxit_;\n"; - mStaticModelFile << " solve_tolf=options_.solve_tolf;\n"; - mStaticModelFile << " y=oo_.steady_state;\n"; - mStaticModelFile << " x=oo_.exo_steady_state;\n"; - mStaticModelFile << " params=M_.params;\n"; - mStaticModelFile << " varargout{2}=1;\n"; - prev_Simulation_Type=-1; - int Blck_Num = 0; - for (i = 0;i < block_triangular.ModelBlock->Size;i++) - { - k = block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && - (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) - skip_head=true; - else - { - skip_head=false; - Blck_Num++; - } - if ((k == EVALUATE_FORWARD || k == EVALUATE_FORWARD_R || k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (!skip_head) - { - if (open_par) - { - mStaticModelFile << " end\n"; - } - mStaticModelFile << " y = " << static_basename << "_" << i + 1 << "(y, x, params, 0);\n"; - } - open_par=false; - } - else if ((k == SOLVE_FORWARD_SIMPLE || k == SOLVE_BACKWARD_SIMPLE || k == SOLVE_FORWARD_COMPLETE || k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_TWO_BOUNDARIES_COMPLETE || k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - { - mStaticModelFile << "end\n"; - } - open_par=false; - mStaticModelFile << " y_index=["; - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - } - mStaticModelFile << " ];\n"; - - - mStaticModelFile << " g1=0;g2=0;g3=0;\n"; - - int nze, m; - for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) - nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; - mStaticModelFile << " [y, info] = solve_one_boundary('" << static_basename << "_" << i + 1 << "'" << - ", y, x, params, y_index, " << nze << - ", 1, " << block_triangular.ModelBlock->Block_List[i].is_linear << - ", " << Blck_Num << ", y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 0, 0);\n"; - mStaticModelFile << " if(info<=0)\n" - << " varagout(2) = 0;\n" - << " varagout(1) = i+1;\n" - << " return;\n" - << " end;\n"; - - } - prev_Simulation_Type=k; - } - if (open_par) - mStaticModelFile << " end;\n"; - mStaticModelFile << " oo_.steady_state = y;\n"; - mStaticModelFile << " if isempty(ys0_)\n"; - mStaticModelFile << " oo_.endo_simul(:,1:M_.maximum_lag) = oo_.steady_state * ones(1,M_.maximum_lag);\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " if(~options_.homotopy_mode)\n"; - mStaticModelFile << " disp('Steady State value');\n"; - mStaticModelFile << " disp([strcat(M_.endo_names,' : ') num2str(oo_.steady_state,'%f')]);\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " varargout{2}=info;\n"; - mStaticModelFile << " varargout{1}=oo_.steady_state;\n"; - mStaticModelFile << "return;\n"; - writeModelStaticEquationsOrdered_M(block_triangular.ModelBlock, static_basename); - mStaticModelFile.close(); - chdir(".."); -} - -void -ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const -{ - string sp; - ofstream mDynamicModelFile; - ostringstream tmp, tmp1, tmp_eq; - int prev_Simulation_Type, tmp_i; - //SymbolicGaussElimination SGE; - bool OK; - chdir(basename.c_str()); - string filename = dynamic_basename + ".m"; - 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 << "%\n"; - mDynamicModelFile << "% " << filename << " : Computes dynamic model for Dynare\n"; - mDynamicModelFile << "%\n"; - mDynamicModelFile << "% Warning : this file is generated automatically by Dynare\n"; - mDynamicModelFile << "% from model file (.mod)\n\n"; - mDynamicModelFile << "%/\n"; - - int i, k, Nb_SGE=0; - bool skip_head, open_par=false; - if (computeJacobian || computeJacobianExo || computeHessian) - { - mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n"; - mDynamicModelFile << " global oo_ options_ M_ ;\n"; - mDynamicModelFile << " g2=[];g3=[];\n"; - //Temporary variables declaration - OK=true; - ostringstream tmp_output; - for (temporary_terms_type::const_iterator it = temporary_terms.begin(); - it != temporary_terms.end(); it++) - { - if (OK) - OK=false; - else - tmp_output << " "; - (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms); - } - if (tmp_output.str().length()>0) - mDynamicModelFile << " global " << tmp_output.str() << " M_ ;\n"; - - mDynamicModelFile << " T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n"; - tmp_output.str(""); - for (temporary_terms_type::const_iterator it = temporary_terms.begin(); - it != temporary_terms.end(); it++) - { - tmp_output << " "; - (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); - tmp_output << "=T_init;\n"; - } - if (tmp_output.str().length()>0) - mDynamicModelFile << tmp_output.str(); - - mDynamicModelFile << " y_kmin=M_.maximum_lag;\n"; - mDynamicModelFile << " y_kmax=M_.maximum_lead;\n"; - mDynamicModelFile << " y_size=M_.endo_nbr;\n"; - mDynamicModelFile << " if(length(varargin)>0)\n"; - mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n"; - mDynamicModelFile << " params=varargin{3};\n"; - mDynamicModelFile << " it_=varargin{4};\n"; - /*i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+ - symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1); - mDynamicModelFile << " g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n";*/ - mDynamicModelFile << " Per_u_=0;\n"; - mDynamicModelFile << " Per_y_=it_*y_size;\n"; - mDynamicModelFile << " y=varargin{1};\n"; - mDynamicModelFile << " ys=y(it_,:);\n"; - mDynamicModelFile << " x=varargin{2};\n"; - prev_Simulation_Type=-1; - tmp.str(""); - tmp_eq.str(""); - for (int count_call=1, i = 0;i < block_triangular.ModelBlock->Size;i++, count_call++) - { - k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) && - ((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R) - || (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))) - { - mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; - tmp_eq.str(""); - mDynamicModelFile << " y_index=[" << tmp.str() << "];\n"; - tmp.str(""); - mDynamicModelFile << tmp1.str(); - tmp1.str(""); - } - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; - } - if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R) - { - if (i==block_triangular.ModelBlock->Size-1) - { - mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; - tmp_eq.str(""); - mDynamicModelFile << " y_index=[" << tmp.str() << "];\n"; - tmp.str(""); - mDynamicModelFile << tmp1.str(); - tmp1.str(""); - } - } - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && - (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) - skip_head=true; - else - skip_head=false; - switch (k) - { - case EVALUATE_FORWARD: - case EVALUATE_BACKWARD: - case EVALUATE_FORWARD_R: - case EVALUATE_BACKWARD_R: - if (!skip_head) - { - tmp1 << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n"; - tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n"; - } - break; - case SOLVE_FORWARD_SIMPLE: - case SOLVE_BACKWARD_SIMPLE: - mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n"; - mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n"; - mDynamicModelFile << " residual(y_index_eq)=r;\n"; - tmp_eq.str(""); - tmp.str(""); - break; - case SOLVE_FORWARD_COMPLETE: - case SOLVE_BACKWARD_COMPLETE: - mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n"; - mDynamicModelFile << " residual(y_index_eq)=r;\n"; - break; - case SOLVE_TWO_BOUNDARIES_COMPLETE: - case SOLVE_TWO_BOUNDARIES_SIMPLE: - int j; - mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; - mDynamicModelFile << " y_index = ["; - for (j=0;j<tmp_i;j++) - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr(); - } - int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1; - for (j=0;j<tmp_ix;j++) - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].nb_exo;ik++) - mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr()+symbol_table.endo_nbr()*tmp_i; - mDynamicModelFile << " ];\n"; - tmp.str(""); - tmp_eq.str(""); - //mDynamicModelFile << " ga = [];\n"; - j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1) - + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1); - /*mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " << - block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";*/ - tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; - mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << variable_table.max_lag << ", 1, " << variable_table.max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n"; - /*if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead) - mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n"; - else - mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";*/ - mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n"; - break; - } - prev_Simulation_Type=k; - } - if (tmp1.str().length()) - { - mDynamicModelFile << tmp1.str(); - tmp1.str(""); - } - mDynamicModelFile << " varargout{1}=residual;\n"; - mDynamicModelFile << " varargout{2}=dr;\n"; - mDynamicModelFile << " return;\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n"; - mDynamicModelFile << " if(options_.simulation_method==0)\n"; - mDynamicModelFile << " mthd='Sparse LU';\n"; - mDynamicModelFile << " elseif(options_.simulation_method==2)\n"; - mDynamicModelFile << " mthd='GMRES';\n"; - mDynamicModelFile << " elseif(options_.simulation_method==3)\n"; - mDynamicModelFile << " mthd='BICGSTAB';\n"; - mDynamicModelFile << " else\n"; - mDynamicModelFile << " mthd='UNKNOWN';\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n"; - mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n"; - mDynamicModelFile << " fprintf('\\n') ;\n"; - mDynamicModelFile << " periods=options_.periods;\n"; - mDynamicModelFile << " maxit_=options_.maxit_;\n"; - mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n"; - mDynamicModelFile << " y=oo_.endo_simul';\n"; - mDynamicModelFile << " x=oo_.exo_simul;\n"; - } - prev_Simulation_Type=-1; - mDynamicModelFile << " params=M_.params;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.status = 0;\n"; - for (i = 0;i < block_triangular.ModelBlock->Size;i++) - { - k = block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && - (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) - skip_head=true; - else - skip_head=false; - if ((k == EVALUATE_FORWARD || k == EVALUATE_FORWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (!skip_head) - { - if (open_par) - { - mDynamicModelFile << " end\n"; - } - mDynamicModelFile << " oo_.deterministic_simulation.status = 1;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.error = 0;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.iterations = 0;\n"; - mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; - mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; - mDynamicModelFile << " else\n"; - mDynamicModelFile << " blck_num = 1;\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).status = 1;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n"; - mDynamicModelFile << " g1=[];g2=[];g3=[];\n"; - //mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n"; - mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n"; - } - //open_par=true; - } - else if ((k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (!skip_head) - { - if (open_par) - { - mDynamicModelFile << " end\n"; - } - mDynamicModelFile << " oo_.deterministic_simulation.status = 1;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.error = 0;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.iterations = 0;\n"; - mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; - mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; - mDynamicModelFile << " else\n"; - mDynamicModelFile << " blck_num = 1;\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).status = 1;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n"; - mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n"; - mDynamicModelFile << " g1=[];g2=[];g3=[];\n"; - mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n"; - } - } - else if ((k == SOLVE_FORWARD_COMPLETE || k == SOLVE_FORWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - mDynamicModelFile << " g1=0;\n"; - mDynamicModelFile << " r=0;\n"; - tmp.str(""); - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - } - mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; - int nze, m; - for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) - nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; - mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; - mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; - mDynamicModelFile << " else\n"; - mDynamicModelFile << " blck_num = 1;\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" << - ", y, x, params, y_index, " << nze << - ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear << - ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n"; - - } - else if ((k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - mDynamicModelFile << " g1=0;\n"; - mDynamicModelFile << " r=0;\n"; - tmp.str(""); - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - } - mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; - int nze, m; - for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) - nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; - mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; - mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; - mDynamicModelFile << " else\n"; - mDynamicModelFile << " blck_num = 1;\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" << - ", y, x, params, y_index, " << nze << - ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear << - ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n"; - } - else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE || k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - Nb_SGE++; - int nze, m; - for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) - nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; - mDynamicModelFile << " y_index=["; - for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++) - { - mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - } - mDynamicModelFile << " ];\n"; - mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; - mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; - mDynamicModelFile << " else\n"; - mDynamicModelFile << " blck_num = 1;\n"; - mDynamicModelFile << " end;\n"; - mDynamicModelFile << " y = solve_two_boundaries('" << dynamic_basename << "_" << i + 1 << "'" << - ", y, x, params, y_index, " << nze << - ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].Max_Lag << - ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << - ", " << block_triangular.ModelBlock->Block_List[i].is_linear << - ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method);\n"; - - } - prev_Simulation_Type=k; - } - if (open_par) - mDynamicModelFile << " end;\n"; - open_par=false; - mDynamicModelFile << " oo_.endo_simul = y';\n"; - mDynamicModelFile << "return;\n"; - - mDynamicModelFile.close(); - - writeModelEquationsOrdered_M( block_triangular.ModelBlock, dynamic_basename); - - chdir(".."); -} - -void -ModelTree::writeDynamicModel(ostream &DynamicOutput) const -{ - ostringstream lsymetric; // Used when writing symetric elements in Hessian - ostringstream model_output; // Used for storing model equations - ostringstream jacobian_output; // Used for storing jacobian equations - ostringstream hessian_output; // Used for storing Hessian equations - ostringstream third_derivatives_output; - - ExprNodeOutputType output_type = (mode == eStandardMode || mode==eSparseMode ? oMatlabDynamicModel : oCDynamicModel); - - writeModelLocalVariables(model_output, output_type); - - writeTemporaryTerms(model_output, output_type); - - writeModelEquations(model_output, output_type); - - int nrows = equations.size(); - int nvars = variable_table.getDynJacobianColsNbr(computeJacobianExo); - int nvars_sq = nvars * nvars; - - // Writing Jacobian - if (computeJacobian || computeJacobianExo) - for (first_derivatives_type::const_iterator it = first_derivatives.begin(); - it != first_derivatives.end(); it++) - { - int eq = it->first.first; - int var = it->first.second; - NodeID d1 = it->second; - - if (computeJacobianExo || variable_table.getType(var) == eEndogenous) - { - ostringstream g1; - g1 << " g1"; - matrixHelper(g1, eq, variable_table.getDynJacobianCol(var), output_type); - - jacobian_output << g1.str() << "=" << g1.str() << "+"; - d1->writeOutput(jacobian_output, output_type, temporary_terms); - jacobian_output << ";" << endl; - } - } - - // Writing Hessian - if (computeHessian) - for (second_derivatives_type::const_iterator it = second_derivatives.begin(); - it != second_derivatives.end(); it++) - { - int eq = it->first.first; - int var1 = it->first.second.first; - int var2 = it->first.second.second; - NodeID d2 = it->second; - - int id1 = variable_table.getDynJacobianCol(var1); - int id2 = variable_table.getDynJacobianCol(var2); - - int col_nb = id1*nvars+id2; - int col_nb_sym = id2*nvars+id1; - - hessian_output << " g2"; - matrixHelper(hessian_output, eq, col_nb, output_type); - hessian_output << " = "; - d2->writeOutput(hessian_output, output_type, temporary_terms); - hessian_output << ";" << endl; - - // Treating symetric elements - if (id1 != id2) - { - lsymetric << " g2"; - matrixHelper(lsymetric, eq, col_nb_sym, output_type); - lsymetric << " = " << "g2"; - matrixHelper(lsymetric, eq, col_nb, output_type); - lsymetric << ";" << endl; - } - } - - // Writing third derivatives - if (computeThirdDerivatives) - for (third_derivatives_type::const_iterator it = third_derivatives.begin(); - it != third_derivatives.end(); it++) - { - int eq = it->first.first; - int var1 = it->first.second.first; - int var2 = it->first.second.second.first; - int var3 = it->first.second.second.second; - NodeID d3 = it->second; - - int id1 = variable_table.getDynJacobianCol(var1); - int id2 = variable_table.getDynJacobianCol(var2); - int id3 = variable_table.getDynJacobianCol(var3); - - // Reference column number for the g3 matrix - int ref_col = id1 * nvars_sq + id2 * nvars + id3; - - third_derivatives_output << " g3"; - matrixHelper(third_derivatives_output, eq, ref_col, output_type); - third_derivatives_output << " = "; - d3->writeOutput(third_derivatives_output, output_type, temporary_terms); - third_derivatives_output << ";" << endl; - - // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal) - set<int> cols; - cols.insert(id1 * nvars_sq + id3 * nvars + id2); - cols.insert(id2 * nvars_sq + id1 * nvars + id3); - cols.insert(id2 * nvars_sq + id3 * nvars + id1); - cols.insert(id3 * nvars_sq + id1 * nvars + id2); - cols.insert(id3 * nvars_sq + id2 * nvars + id1); - - for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++) - if (*it2 != ref_col) - { - third_derivatives_output << " g3"; - matrixHelper(third_derivatives_output, eq, *it2, output_type); - third_derivatives_output << " = " << "g3"; - matrixHelper(third_derivatives_output, eq, ref_col, output_type); - third_derivatives_output << ";" << endl; - } - } - - if (mode == eStandardMode) - { - DynamicOutput << "%" << endl - << "% Model equations" << endl - << "%" << endl - << endl - << "residual = zeros(" << nrows << ", 1);" << endl - << model_output.str(); - - if (computeJacobian || computeJacobianExo) - { - // Writing initialization instruction for matrix g1 - DynamicOutput << "if nargout >= 2," << endl - << " g1 = zeros(" << nrows << ", " << nvars << ");" << endl - << endl - << "%" << endl - << "% Jacobian matrix" << endl - << "%" << endl - << endl - << jacobian_output.str() - << "end" << endl; - } - if (computeHessian) - { - // Writing initialization instruction for matrix g2 - int ncols = nvars_sq; - DynamicOutput << "if nargout >= 3," << endl - << " g2 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl - << endl - << "%" << endl - << "% Hessian matrix" << endl - << "%" << endl - << endl - << hessian_output.str() - << lsymetric.str() - << "end;" << endl; - } - if (computeThirdDerivatives) - { - int ncols = nvars_sq * nvars; - DynamicOutput << "if nargout >= 4," << endl - << " g3 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl - << endl - << "%" << endl - << "% Third order derivatives" << endl - << "%" << endl - << endl - << third_derivatives_output.str() - << "end;" << endl; - } - } - else - { - DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *g2)" << endl - << "{" << endl - << " double lhs, rhs;" << endl - << endl - << " /* Residual equations */" << endl - << model_output.str(); - - if (computeJacobian || computeJacobianExo) - { - DynamicOutput << " /* Jacobian */" << endl - << " if (g1 == NULL)" << endl - << " return;" << endl - << " else" << endl - << " {" << endl - << jacobian_output.str() - << " }" << endl; - } - if (computeHessian) - { - DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl - << " if (g2 == NULL)" << endl - << " return;" << endl - << " else" << endl - << " {" << endl - << hessian_output.str() - << lsymetric.str() - << " }" << endl; - } - DynamicOutput << "}" << endl << endl; - } -} - -void -ModelTree::writeOutput(ostream &output) const -{ - /* Writing initialisation for M_.lead_lag_incidence matrix - M_.lead_lag_incidence is a matrix with as many columns as there are - endogenous variables and as many rows as there are periods in the - models (nbr of rows = M_.max_lag+M_.max_lead+1) - - The matrix elements are equal to zero if a variable isn't present in the - model at a given period. - */ - output << "M_.lead_lag_incidence = ["; - // Loop on endogenous variables - int lag = 0; - for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++) - { - output << "\n\t"; - // Loop on periods - for (lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++) - { - // Print variableID if exists with current period, otherwise print 0 - try - { - int varID = variable_table.getID(symbol_table.getID(eEndogenous, endoID), lag); - output << " " << variable_table.getDynJacobianCol(varID) + 1; - } - catch (VariableTable::UnknownVariableKeyException &e) - { - output << " 0"; - } - } - output << ";"; - } - output << "]';\n"; - //In case of sparse model, writes the block structure of the model - - if (mode==eSparseMode || mode==eSparseDLLMode) - { - //int prev_Simulation_Type=-1; - //bool skip_the_head; - int k=0; - int count_lead_lag_incidence = 0; - int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo; - for (int j = 0;j < block_triangular.ModelBlock->Size;j++) - { - //For a block composed of a single equation determines wether we have to evaluate or to solve the equation - //skip_the_head=false; - k++; - count_lead_lag_incidence = 0; - int Block_size=block_triangular.ModelBlock->Block_List[j].Size; - max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ; - max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead; - max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ; - max_lead_endo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Endo; - max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ; - max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo; - bool evaluate=false; - vector<int> exogenous; - vector<int>::iterator it_exogenous; - exogenous.clear(); - ostringstream tmp_s, tmp_s_eq; - tmp_s.str(""); - tmp_s_eq.str(""); - for (int i=0;i<block_triangular.ModelBlock->Block_List[j].Size;i++) - { - tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1; - tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1; - } - for (int i=0;i<block_triangular.ModelBlock->Block_List[j].nb_exo;i++) - { - int ii=block_triangular.ModelBlock->Block_List[j].Exogenous[i]; - for (it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) /*cout << "*it_exogenous=" << *it_exogenous << "\n"*/; - if (it_exogenous==exogenous.end() || exogenous.begin()==exogenous.end()) - exogenous.push_back(ii); - } - output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n"; - output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_lead = " << max_lead << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag_endo << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead_endo << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_exo_lag = " << max_lag_exo << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_exo_lead = " << max_lead_exo << ";\n"; - output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n"; - output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n"; - output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n"; - output << "M_.block_structure.block(" << k << ").exogenous = ["; - int i=0; - for (it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) - if (*it_exogenous>=0) - { - output << " " << *it_exogenous+1; - i++; - } - output << "];\n"; - output << "M_.block_structure.block(" << k << ").exo_nbr = " << i << ";\n"; - - output << "M_.block_structure.block(" << k << ").exo_det_nbr = " << block_triangular.ModelBlock->Block_List[j].nb_exo_det << ";\n"; - - tmp_s.str(""); - - bool done_IM=false; - if (!evaluate) - { - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n"; - for (int l=-max_lag_endo;l<max_lead_endo+1;l++) - { - bool *tmp_IM; - tmp_IM=block_triangular.incidencematrix.Get_IM(l, eEndogenous); - if (tmp_IM) - { - for (int l_var=0;l_var<block_triangular.ModelBlock->Block_List[j].Size;l_var++) - { - for (int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[j].Size;l_equ++) - if (tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr()+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) - { - count_lead_lag_incidence++; - if (tmp_s.str().length()) - tmp_s << " "; - tmp_s << count_lead_lag_incidence; - done_IM=true; - break; - } - if (!done_IM) - tmp_s << " 0"; - done_IM=false; - } - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; - tmp_s.str(""); - } - } - } - else - { - bool done_some_where; - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; - for (int l=-max_lag_endo;l<max_lead_endo+1;l++) - { - bool not_increm=true; - bool *tmp_IM; - tmp_IM=block_triangular.incidencematrix.Get_IM(l, eEndogenous); - int ii=j; - if (tmp_IM) - { - done_some_where = false; - while (ii-j<Block_size) - { - for (int l_var=0;l_var<block_triangular.ModelBlock->Block_List[ii].Size;l_var++) - { - for (int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[ii].Size;l_equ++) - if (tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr()+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) - { - //if(not_increm && l==-max_lag) - count_lead_lag_incidence++; - not_increm=false; - if (tmp_s.str().length()) - tmp_s << " "; - //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; - tmp_s << count_lead_lag_incidence; - done_IM=true; - break; - } - if (!done_IM) - tmp_s << " 0"; - else - done_some_where = true; - done_IM=false; - } - ii++; - } - output << tmp_s.str() << "\n"; - tmp_s.str(""); - } - } - output << "];\n"; - } - - } - for (int j=-block_triangular.incidencematrix.Model_Max_Lag_Endo;j<=block_triangular.incidencematrix.Model_Max_Lead_Endo;j++) - { - bool* IM = block_triangular.incidencematrix.Get_IM(j, eEndogenous); - if (IM) - { - bool new_entry=true; - output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").lead_lag = " << j << ";\n"; - output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").sparse_IM = ["; - for (int i=0;i<symbol_table.endo_nbr()*symbol_table.endo_nbr();i++) - { - if (IM[i]) - { - if (!new_entry) - output << " ; "; - else - output << " "; - output << i/symbol_table.endo_nbr()+1 << " " << i % symbol_table.endo_nbr()+1; - new_entry=false; - } - } - output << "];\n"; - } - } - } - // Writing initialization for some other variables - output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];\n"; - output << "M_.maximum_lag = " << variable_table.max_lag << ";\n"; - output << "M_.maximum_lead = " << variable_table.max_lead << ";\n"; - if (symbol_table.endo_nbr()) - { - output << "M_.maximum_endo_lag = " << variable_table.max_endo_lag << ";\n"; - output << "M_.maximum_endo_lead = " << variable_table.max_endo_lead << ";\n"; - output << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);\n"; - } - if (symbol_table.exo_nbr()) - { - output << "M_.maximum_exo_lag = " << variable_table.max_exo_lag << ";\n"; - output << "M_.maximum_exo_lead = " << variable_table.max_exo_lead << ";\n"; - output << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);\n"; - } - if (symbol_table.exo_det_nbr()) - { - output << "M_.maximum_exo_det_lag = " << variable_table.max_exo_det_lag << ";\n"; - output << "M_.maximum_exo_det_lead = " << variable_table.max_exo_det_lead << ";\n"; - output << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);\n"; - } - if (symbol_table.param_nbr()) - output << "M_.params = repmat(NaN," << symbol_table.param_nbr() << ", 1);\n"; -} - -void -ModelTree::addEquation(NodeID eq) -{ - BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq); - - if (beq == NULL || beq->op_code != oEqual) - { - cerr << "ModelTree::addEquation: you didn't provide an equal node!" << endl; + cerr << "ModelTree::addEquation: you didn't provide an equal node!" << endl; exit(EXIT_FAILURE); } equations.push_back(beq); } -void -ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m) -{ - int i=0; - int j=0; - bool *IM=NULL; - int a_variable_lag=-9999; - for (first_derivatives_type::iterator it = first_derivatives.begin(); - it != first_derivatives.end(); it++) - { - //cout << "it->first.second=" << it->first.second << " variable_table.getSymbolID(it->first.second)=" << variable_table.getSymbolID(it->first.second) << " Type=" << variable_table.getType(it->first.second) << " eEndogenous=" << eEndogenous << " eExogenous=" << eExogenous << " variable_table.getLag(it->first.second)=" << variable_table.getLag(it->first.second) << "\n"; - if (variable_table.getType(it->first.second) == eEndogenous) - { - NodeID Id = it->second; - double val = 0; - try - { - val = Id->eval(eval_context); - } - catch (ExprNode::EvalException &e) - { - cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ") [" << variable_table.getSymbolID(it->first.second) << "] !" << endl; - Id->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms); - cout << "\n"; - cerr << "ModelTree::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ")!" << endl; - } - int eq=it->first.first; - int var=symbol_table.getTypeSpecificID(variable_table.getSymbolID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second); - int k1=variable_table.getLag(it->first.second); - if (a_variable_lag!=k1) - { - IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous); - a_variable_lag=k1; - } - if (k1==0) - { - j++; - (*j_m)[make_pair(eq,var)]=val; - } - if (IM[eq*symbol_table.endo_nbr()+var] && (fabs(val) < cutoff)) - { - if (block_triangular.bt_verbose) - cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")\n"; - block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous); - i++; - } - } - } - if (i>0) - { - cout << i << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded\n"; - cout << "the contemporaneous incidence matrix has " << j << " elements\n"; - } -} - -void -ModelTree::BlockLinear(Model_Block *ModelBlock) -{ - int i,j,l,m,ll; - for (j = 0;j < ModelBlock->Size;j++) - { - if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || - ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE) - { - ll=ModelBlock->Block_List[j].Max_Lag; - for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[ll].size;i++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[ll].Equ_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[ll].Var_Index[i]; - //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,0))); - first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var),0))); - if (it!= first_derivatives.end()) - { - NodeID Id = it->second; - set<pair<int, int> > endogenous; - Id->collectEndogenous(endogenous); - if (endogenous.size() > 0) - { - for (l=0;l<ModelBlock->Block_List[j].Size;l++) - { - if (endogenous.find(make_pair(ModelBlock->Block_List[j].Variable[l], 0)) != endogenous.end()) - { - ModelBlock->Block_List[j].is_linear=false; - goto follow; - } - } - } - } - } - } - else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) - { - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - int k1=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]; - //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,k1))); - first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var),k1))); - NodeID Id = it->second; - if (it!= first_derivatives.end()) - { - set<pair<int, int> > endogenous; - Id->collectEndogenous(endogenous); - if (endogenous.size() > 0) - { - for (l=0;l<ModelBlock->Block_List[j].Size;l++) - { - if (endogenous.find(make_pair(ModelBlock->Block_List[j].Variable[l], k1)) != endogenous.end()) - { - ModelBlock->Block_List[j].is_linear=false; - goto follow; - } - } - } - } - } - } - } - follow: - i=0; - } -} - -void -ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_terms) -{ - cout << equations.size() << " equation(s) found" << endl; - - // Computes dynamic jacobian columns - variable_table.computeDynJacobianCols(); - - // Determine derivation order - int order = 1; - if (computeThirdDerivatives) - order = 3; - else if (computeHessian || computeStaticHessian) - order = 2; - - // Launch computations - derive(order); - - if (mode == eSparseDLLMode || mode == eSparseMode) - { - BuildIncidenceMatrix(); - - jacob_map j_m; - - evaluateJacobian(eval_context, &j_m); - - - if (block_triangular.bt_verbose) - { - cout << "The gross incidence matrix \n"; - block_triangular.incidencematrix.Print_IM(eEndogenous); - } - block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations); - BlockLinear(block_triangular.ModelBlock); - if (!no_tmp_terms) - computeTemporaryTermsOrdered(order, block_triangular.ModelBlock); - } - else - if (!no_tmp_terms) - computeTemporaryTerms(order); -} - -void -ModelTree::writeStaticFile(const string &basename) const -{ - switch (mode) - { - case eStandardMode: - /*case eSparseDLLMode:*/ - writeStaticMFile(basename + "_static"); - break; - case eSparseDLLMode: - case eSparseMode: - // create a directory to store all files -#ifdef _WIN32 - mkdir(basename.c_str()); -#else - mkdir(basename.c_str(), 0777); -#endif - - writeSparseStaticMFile(basename + "_static", basename, mode); - break; - case eDLLMode: - writeStaticCFile(basename + "_static"); - break; - } -} - -void -ModelTree::writeDynamicFile(const string &basename) const -{ - switch (mode) - { - case eStandardMode: - writeDynamicMFile(basename + "_dynamic"); - break; - case eSparseMode: - writeSparseDynamicMFile(basename + "_dynamic", basename, mode); - block_triangular.Free_Block(block_triangular.ModelBlock); - block_triangular.incidencematrix.Free_IM(); - //block_triangular.Free_IM_X(block_triangular.First_IM_X); - break; - case eDLLMode: - writeDynamicCFile(basename + "_dynamic"); - break; - case eSparseDLLMode: - // create a directory to store all the files -#ifdef _WIN32 - mkdir(basename.c_str()); -#else - mkdir(basename.c_str(), 0777); -#endif - writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL, map_idx); - block_triangular.Free_Block(block_triangular.ModelBlock); - block_triangular.incidencematrix.Free_IM(); - //block_triangular.Free_IM_X(block_triangular.First_IM_X); - break; - } -} - void ModelTree::matrixHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const { diff --git a/ModelTree.hh b/ModelTree.hh index b5ff7fb4bef67eb535002c464e66836548dceff9..330bc8ead3f6e22f794d52530f314d23b2bcb17c 100644 --- a/ModelTree.hh +++ b/ModelTree.hh @@ -26,12 +26,8 @@ using namespace std; #include <vector> #include <map> #include <ostream> -#include <algorithm> -#include "SymbolTable.hh" -#include "NumericalConstants.hh" #include "DataTree.hh" -#include "BlockTriangular.hh" //! The three in which ModelTree can work enum ModelTreeMode @@ -42,10 +38,10 @@ enum ModelTreeMode eSparseDLLMode //!< Sparse DLL mode (static file in Matlab, dynamic file in C with block decomposition plus a binary file) }; -//! Stores a model's equations and derivatives +//! Shared code for static and dynamic models class ModelTree : public DataTree { -private: +protected: //! Stores declared equations vector<BinaryOpNode *> equations; @@ -77,19 +73,13 @@ private: //! Temporary terms (those which will be noted Txxxx) temporary_terms_type temporary_terms; - map_idx_type map_idx; //! Computes derivatives of ModelTree void derive(int order); //! Write derivative of an equation w.r. to a variable void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; - //! Write derivative code of an equation w.r. to a variable - void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const; //! Computes temporary terms void computeTemporaryTerms(int order); - void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock); - //! Build The incidence matrix form the modeltree - void BuildIncidenceMatrix(); //! Writes temporary terms void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const; //! Writes model local variables @@ -97,39 +87,6 @@ private: void writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const; //! Writes model equations void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const; - //! Writes the static model equations and its derivatives - /*! \todo handle hessian in C output */ - void writeStaticModel(ostream &StaticOutput) const; - //! Writes the dynamic model equations and its derivatives - /*! \todo add third derivatives handling in C output */ - void writeDynamicModel(ostream &DynamicOutput) const; - //! Writes the Block reordred structure of the model in M output - void writeModelEquationsOrdered_M(Model_Block *ModelBlock, const string &dynamic_basename) const; - //! Writes the Block reordred structure of the static model in M output - void writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const string &static_basename) const; - //! Writes the code of the Block reordred structure of the model in virtual machine bytecode - void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const; - //! Writes static model file (Matlab version) - void writeStaticMFile(const string &static_basename) const; - //! Writes static model file (C version) - void writeStaticCFile(const string &static_basename) const; - //! Writes static model file when Sparse option is on (Matlab version) - void writeSparseStaticMFile(const string &static_basename, const string &basename, const int mode) const; - //! Writes dynamic model file (Matlab version) - void writeDynamicMFile(const string &dynamic_basename) const; - //! Writes dynamic model file (C version) - /*! \todo add third derivatives handling */ - void writeDynamicCFile(const string &dynamic_basename) const; - //! Writes dynamic model file when SparseDLL option is on - void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const; - //! Computes jacobian and prepares for equation normalization - /*! Using values from initval/endval blocks and parameter initializations: - - computes the jacobian for the model w.r. to contemporaneous variables - - removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff) - */ - void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m); - void BlockLinear(Model_Block *ModelBlock); - string reform(string name) const; //! Writes either (i+1,j+1) or [i+j*n_i] whether we are in Matlab or C mode void matrixHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const; @@ -138,41 +95,8 @@ public: ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); //! Mode in which the ModelTree is supposed to work (Matlab, DLL or SparseDLL) ModelTreeMode mode; - //! Absolute value under which a number is considered to be zero - double cutoff; - //! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc) - double markowitz; - //! Use a graphical and symbolic version of the symbolic gaussian elimination (new_SGE = false) or use direct gaussian elimination (new_SGE = true) - bool new_SGE; - //! the file containing the model and the derivatives code - ofstream code_file; //! Declare a node as an equation of the model void addEquation(NodeID eq); - //! Whether dynamic Jacobian (w.r. to endogenous) should be written - bool computeJacobian; - //! Whether dynamic Jacobian (w.r. to endogenous and exogenous) should be written - bool computeJacobianExo; - //! Whether dynamic Hessian (w.r. to endogenous and exogenous) should be written - bool computeHessian; - //! Whether static Hessian (w.r. to endogenous only) should be written - bool computeStaticHessian; - //! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written - bool computeThirdDerivatives; - //! Execute computations (variable sorting + derivation) - /*! You must set computeJacobian, computeJacobianExo, computeHessian, computeStaticHessian and computeThirdDerivatives to correct values before calling this function - \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */ - void computingPass(const eval_context_type &eval_context, bool no_tmp_terms); - //! Writes model initialization and lead/lag incidence matrix to output - void writeOutput(ostream &output) const; - //! Writes static model file - void writeStaticFile(const string &basename) const; - //! Writes dynamic model file - void writeDynamicFile(const string &basename) const; - //! Complete set to block decompose the model - BlockTriangular block_triangular; - //! Adds informations for simulation in a binary file - void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, - const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const; //! Returns the number of equations in the model int equation_number() const; }; diff --git a/ParsingDriver.cc b/ParsingDriver.cc index 159e9efcea186e4c0b85fa0a4ee8b8d759d22ad2..26fe3942e9ee1d4b04f45ef65fa05321f29e3a4d 100644 --- a/ParsingDriver.cc +++ b/ParsingDriver.cc @@ -392,7 +392,7 @@ ParsingDriver::end_homotopy() void ParsingDriver::begin_model() { - set_current_data_tree(&mod_file->model_tree); + set_current_data_tree(&mod_file->dynamic_model); } void @@ -585,10 +585,7 @@ ParsingDriver::add_to_row(NodeID v) void ParsingDriver::steady() { - if (mod_file->model_tree.mode == eSparseMode || mod_file->model_tree.mode == eSparseDLLMode) - mod_file->addStatement(new SteadySparseStatement(options_list)); - else - mod_file->addStatement(new SteadyStatement(options_list)); + mod_file->addStatement(new SteadyStatement(options_list)); options_list.clear(); } @@ -619,10 +616,10 @@ ParsingDriver::option_num(const string &name_option, const string &opt) && (options_list.num_options.find(name_option) != options_list.num_options.end())) error("option " + name_option + " declared twice"); - if ((name_option == "periods") && (mod_file->model_tree.mode == eSparseDLLMode || mod_file->model_tree.mode == eSparseMode)) - mod_file->model_tree.block_triangular.periods = atoi(opt.c_str()); + if ((name_option == "periods") && (mod_file->dynamic_model.mode == eSparseDLLMode || mod_file->dynamic_model.mode == eSparseMode)) + mod_file->dynamic_model.block_triangular.periods = atoi(opt.c_str()); else if (name_option == "cutoff") - mod_file->model_tree.cutoff = atof(opt.c_str()); + mod_file->dynamic_model.cutoff = atof(opt.c_str()); options_list.num_options[name_option] = opt; } @@ -686,7 +683,7 @@ void ParsingDriver::stoch_simul() void ParsingDriver::simulate() { - if (mod_file->model_tree.mode == eSparseDLLMode || mod_file->model_tree.mode == eSparseMode) + if (mod_file->dynamic_model.mode == eSparseDLLMode || mod_file->dynamic_model.mode == eSparseMode) simul_sparse(); else simul(); @@ -695,7 +692,7 @@ void ParsingDriver::simulate() void ParsingDriver::simul_sparse() { - mod_file->addStatement(new SimulSparseStatement(options_list, mod_file->model_tree.mode)); + mod_file->addStatement(new SimulSparseStatement(options_list, mod_file->dynamic_model.mode)); options_list.clear(); } @@ -1033,7 +1030,7 @@ ParsingDriver::run_model_comparison() void ParsingDriver::begin_planner_objective() { - set_current_data_tree(new ModelTree(mod_file->symbol_table, mod_file->num_constants)); + set_current_data_tree(new StaticModel(mod_file->symbol_table, mod_file->num_constants)); } void @@ -1043,7 +1040,7 @@ ParsingDriver::end_planner_objective(NodeID expr) NodeID eq = model_tree->AddEqual(expr, model_tree->Zero); model_tree->addEquation(eq); - mod_file->addStatement(new PlannerObjectiveStatement(model_tree)); + mod_file->addStatement(new PlannerObjectiveStatement(dynamic_cast<StaticModel *>(model_tree))); reset_data_tree(); } @@ -1120,7 +1117,7 @@ ParsingDriver::change_type(SymbolType new_type, vector<string *> *var_list) // Check if symbol already used in a VariableNode if (mod_file->expressions_tree.isSymbolUsed(id) - || mod_file->model_tree.isSymbolUsed(id)) + || mod_file->dynamic_model.isSymbolUsed(id)) error("You cannot modify the type of symbol " + **it + " after having used it in an expression"); mod_file->symbol_table.changeType(id, new_type); diff --git a/StaticModel.cc b/StaticModel.cc new file mode 100644 index 0000000000000000000000000000000000000000..24fb837adec0d7fa5417d3fe29c94c8525b763a4 --- /dev/null +++ b/StaticModel.cc @@ -0,0 +1,289 @@ +/* + * 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 "StaticModel.hh" + +StaticModel::StaticModel(SymbolTable &symbol_table_arg, + NumericalConstants &num_constants_arg) : + ModelTree(symbol_table_arg, num_constants_arg), + computeStaticHessian(false) +{ +} + +void +StaticModel::writeStaticMFile(const string &static_basename) const +{ + string filename = static_basename + ".m"; + + ofstream mStaticModelFile; + mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); + if (!mStaticModelFile.is_open()) + { + cerr << "Error: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + // Writing comments and function definition command + mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "(y, x, params)" << endl + << "%" << endl + << "% Status : Computes static model for Dynare" << endl + << "%" << endl + << "% Warning : this file is generated automatically by Dynare" << endl + << "% from model file (.mod)" << endl << endl; + + writeStaticModel(mStaticModelFile); + + mStaticModelFile.close(); +} + +void +StaticModel::writeStaticCFile(const string &static_basename) const +{ + string filename = static_basename + ".c"; + + ofstream mStaticModelFile; + mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); + if (!mStaticModelFile.is_open()) + { + cerr << "Error: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + mStaticModelFile << "/*" << endl + << " * " << filename << " : Computes static model for Dynare" << 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; + + // Writing the function Static + writeStaticModel(mStaticModelFile); + + // Writing the gateway routine + mStaticModelFile << "/* 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;" << 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 + << " 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() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl + << " /* Create a C pointer to a copy of the output matrix g1. */" << endl + << " g1 = mxGetPr(plhs[1]);" << endl + << " }" << endl + << endl + << " /* Call the C Static. */" << endl + << " Static(y, x, params, residual, g1);" << endl + << "}" << endl; + + mStaticModelFile.close(); +} + +void +StaticModel::writeStaticModel(ostream &StaticOutput) const +{ + ostringstream model_output; // Used for storing model equations + ostringstream jacobian_output; // Used for storing jacobian equations + ostringstream hessian_output; + ostringstream lsymetric; // For symmetric elements in hessian + + ExprNodeOutputType output_type = (mode == eDLLMode ? oCStaticModel : oMatlabStaticModel); + + writeModelLocalVariables(model_output, output_type); + + writeTemporaryTerms(model_output, output_type); + + writeModelEquations(model_output, output_type); + + // Write Jacobian w.r. to endogenous only + for (first_derivatives_type::const_iterator it = first_derivatives.begin(); + it != first_derivatives.end(); it++) + { + int eq = it->first.first; + int var = it->first.second; + NodeID d1 = it->second; + + if (variable_table.getType(var) == eEndogenous) + { + ostringstream g1; + g1 << " g1"; + matrixHelper(g1, eq, symbol_table.getTypeSpecificID(variable_table.getSymbolID(var)), output_type); + + jacobian_output << g1.str() << "=" << g1.str() << "+"; + d1->writeOutput(jacobian_output, output_type, temporary_terms); + jacobian_output << ";" << endl; + } + } + + // Write Hessian w.r. to endogenous only + if (computeStaticHessian) + for (second_derivatives_type::const_iterator it = second_derivatives.begin(); + it != second_derivatives.end(); it++) + { + int eq = it->first.first; + int var1 = it->first.second.first; + int var2 = it->first.second.second; + NodeID d2 = it->second; + + // Keep only derivatives w.r. to endogenous variables + if (variable_table.getType(var1) == eEndogenous + && variable_table.getType(var2) == eEndogenous) + { + int id1 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var1)); + int id2 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var2)); + + int col_nb = id1*symbol_table.endo_nbr()+id2; + int col_nb_sym = id2*symbol_table.endo_nbr()+id1; + + hessian_output << " g2"; + matrixHelper(hessian_output, eq, col_nb, output_type); + hessian_output << " = "; + d2->writeOutput(hessian_output, output_type, temporary_terms); + hessian_output << ";" << endl; + + // Treating symetric elements + if (var1 != var2) + { + lsymetric << " g2"; + matrixHelper(lsymetric, eq, col_nb_sym, output_type); + lsymetric << " = " << "g2"; + matrixHelper(lsymetric, eq, col_nb, output_type); + lsymetric << ";" << endl; + } + } + } + + // Writing ouputs + if (mode != eDLLMode) + { + StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl + << "%" << endl + << "% Model equations" << endl + << "%" << endl + << endl + << model_output.str() + << "if ~isreal(residual)" << endl + << " residual = real(residual)+imag(residual).^2;" << endl + << "end" << endl + << "if nargout >= 2," << endl + << " g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl + << endl + << "%" << endl + << "% Jacobian matrix" << endl + << "%" << endl + << endl + << jacobian_output.str() + << " if ~isreal(g1)" << endl + << " g1 = real(g1)+2*imag(g1);" << endl + << " end" << endl + << "end" << endl; + if (computeStaticHessian) + { + StaticOutput << "if nargout >= 3,\n"; + // Writing initialization instruction for matrix g2 + int ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr(); + StaticOutput << " g2 = sparse([],[],[], " << equations.size() << ", " << ncols << ", " << 5*ncols << ");" << endl + << endl + << "%" << endl + << "% Hessian matrix" << endl + << "%" << endl + << endl + << hessian_output.str() + << lsymetric.str() + << "end;" << endl; + } + } + else + { + StaticOutput << "void Static(double *y, double *x, double *params, double *residual, double *g1)" << endl + << "{" << endl + << " double lhs, rhs;" << endl + // Writing residual equations + << " /* Residual equations */" << endl + << " if (residual == NULL)" << endl + << " return;" << endl + << " else" << endl + << " {" << endl + << model_output.str() + // Writing Jacobian + << " /* Jacobian for endogenous variables without lag */" << endl + << " if (g1 == NULL)" << endl + << " return;" << endl + << " else" << endl + << " {" << endl + << jacobian_output.str() + << " }" << endl + << " }" << endl + << "}" << endl << endl; + } +} + +void +StaticModel::writeStaticFile(const string &basename) const +{ + switch (mode) + { + case eStandardMode: + case eSparseDLLMode: + case eSparseMode: + writeStaticMFile(basename + "_static"); + break; + case eDLLMode: + writeStaticCFile(basename + "_static"); + break; + } +} + +void +StaticModel::computingPass(bool no_tmp_terms) +{ + // Determine derivation order + int order = 1; + if (computeStaticHessian) + order = 2; + + // Launch computations + cout << "Computing static model derivatives at order " << order << "..."; + derive(order); + cout << " done" << endl; + + if (!no_tmp_terms) + computeTemporaryTerms(order); +} diff --git a/StaticModel.hh b/StaticModel.hh new file mode 100644 index 0000000000000000000000000000000000000000..0a2c62ccedc362fa9cf15baef4d3438c537baf39 --- /dev/null +++ b/StaticModel.hh @@ -0,0 +1,51 @@ +/* + * 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/>. + */ + +#ifndef _STATICMODEL_HH +#define _STATICMODEL_HH + +using namespace std; + +#include "ModelTree.hh" + +//! Stores a static model +class StaticModel : public ModelTree +{ +private: + //! Writes the static model equations and its derivatives + /*! \todo handle hessian in C output */ + void writeStaticModel(ostream &StaticOutput) const; + //! Writes static model file (Matlab version) + void writeStaticMFile(const string &static_basename) const; + //! Writes static model file (C version) + void writeStaticCFile(const string &static_basename) const; + +public: + StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); + //! Whether static Hessian (w.r. to endogenous only) should be written + bool computeStaticHessian; + //! Execute computations (derivation) + /*! You must set computeStaticHessian before calling this function + \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */ + void computingPass(bool no_tmp_terms); + //! Writes static model file + void writeStaticFile(const string &basename) const; +}; + +#endif