diff --git a/StaticDllModel.cc b/StaticDllModel.cc new file mode 100644 index 00000000..da6cf823 --- /dev/null +++ b/StaticDllModel.cc @@ -0,0 +1,1434 @@ +/* + * 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 . + */ + +#include +#include +#include +#include +#include +#include +#include "StaticDllModel.hh" + +// For mkdir() and chdir() +#ifdef _WIN32 +# include +#else +# include +# include +# include +#endif + +StaticDllModel::StaticDllModel(SymbolTable &symbol_table_arg, + NumericalConstants &num_constants_arg) : + ModelTree(symbol_table_arg, num_constants_arg), + max_lag(0), max_lead(0), + max_endo_lag(0), max_endo_lead(0), + max_exo_lag(0), max_exo_lead(0), + max_exo_det_lag(0), max_exo_det_lead(0), + dynJacobianColsNbr(0), + mode(eStandardMode), + cutoff(1e-15), + markowitz(0.7), + mfs(0), + block_triangular(symbol_table_arg, num_constants_arg) +{ +} + +NodeID +StaticDllModel::AddVariable(const string &name, int lag) +{ + return AddVariableInternal(name, lag); +} + +void +StaticDllModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const + { + //first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag))); + first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag))); + if (it != first_derivatives.end()) + (it->second)->compile(code_file, false, temporary_terms, map_idx, false); + else + code_file.write(&FLDZ, sizeof(FLDZ)); + } + + +void +StaticDllModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, map_idx_type &map_idx) const +{ + map >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag))); + if (it != first_chain_rule_derivatives.end()) + (it->second)->compile(code_file, false, temporary_terms, map_idx, false); + else + code_file.write(&FLDZ, sizeof(FLDZ)); +} + + +void +StaticDllModel::BuildIncidenceMatrix() +{ + set > 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 >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++) + { + block_triangular.incidencematrix.fill_IM(eq, it_endogenous->first, 0, eEndogenous); + } + exogenous.clear(); + Id = eq_node->get_arg1(); + Id->collectExogenous(exogenous); + Id = eq_node->get_arg2(); + Id->collectExogenous(exogenous); + for (set >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) + { + block_triangular.incidencematrix.fill_IM(eq, it_exogenous->first, 0, eExogenous); + } + } +} + +void +StaticDllModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock) +{ + map > first_occurence; + map reference_count; + int i, j, eqr, varr, lag; + temporary_terms_type vect; + ostringstream tmp_output; + BinaryOpNode *eq_node; + first_derivatives_type::const_iterator it; + first_chain_rule_derivatives_type::const_iterator it_chr; + 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++) + { + if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S && iBlock_List[j].Nb_Recursives && ModelBlock->Block_List[j].Equation_Normalized[i]) + ModelBlock->Block_List[j].Equation_Normalized[i]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx); + else + { + eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; + eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx); + } + } + for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++) + { + pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i); + lag=it.first.first; + int eqr=it.second.first; + int varr=it.second.second; + it_chr=first_chain_rule_derivatives.find(make_pair(eqr, make_pair( varr, lag))); + it_chr->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++) + { + // Collecte the temporary terms reordered + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S && iBlock_List[j].Nb_Recursives && ModelBlock->Block_List[j].Equation_Normalized[i]) + ModelBlock->Block_List[j].Equation_Normalized[i]->collectTemporary_terms(temporary_terms, ModelBlock, j); + else + { + eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; + eq_node->collectTemporary_terms(temporary_terms, ModelBlock, j); + } + + } + for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++) + { + pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i); + lag=it.first.first; + eqr=it.second.first; + varr=it.second.second; + it_chr=first_chain_rule_derivatives.find(make_pair(eqr, make_pair( varr, lag))); + it_chr->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 +StaticDllModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &static_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; + map reference_count; + ofstream output; + int nze, nze_exo, nze_other_endo; + vector feedback_variables; + //For each block + for (j = 0;j < ModelBlock->Size;j++) + { + //recursive_variables.clear(); + feedback_variables.clear(); + //For a block composed of a single equation determines wether we have to evaluate or to solve the equation + nze = 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; + 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) + output << "function y = " << static_basename << "_" << j+1 << "(y, x, params)\n"; + else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE + || ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE + || ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE + || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) + output << "function [residual, y, g1] = " << static_basename << "_" << j+1 << "(y, x, params)\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) + output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives + << ", " << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives << ", " << nze << ");\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"; + } + if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) + output << " residual=zeros(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives << ",1);\n"; + // 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; + } + string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, 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(""); + /*if((ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) and (iBlock_List[j].Nb_Recursives)) + lhs->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms); + else*/ + lhs->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms); + switch (ModelBlock->Block_List[j].Simulation_Type) + { + case EVALUATE_BACKWARD: + case EVALUATE_FORWARD: +evaluation: if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) + output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel + << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl; + output << " "; + if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE) + { + output << tmp_output.str(); + output << " = "; + rhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); + } + else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S) + { + output << "%" << tmp_output.str(); + output << " = "; + if (ModelBlock->Block_List[j].Equation_Normalized[i]) + { + rhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); + output << "\n "; + tmp_output.str(""); + eq_node = (BinaryOpNode *)ModelBlock->Block_List[j].Equation_Normalized[i]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + lhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); + output << " = "; + rhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); + } + } + else + { + cerr << "Type missmatch for equation " << ModelBlock->Block_List[j].Equation[i]+1 << "\n"; + exit(EXIT_FAILURE); + } + output << ";\n"; + break; + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + if (iBlock_List[j].Nb_Recursives) + goto evaluation; + feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]); + output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel + << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl; + output << " " << "residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ("; + goto end; + default: +end: + output << tmp_output.str(); + output << ") - ("; + rhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); + output << ");\n"; + } + } + // 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: + break; + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++) + { + pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i); + k=it.first.first; + int eq=it.first.second.first; + int var=it.first.second.second; + int eqr=it.second.first; + int varr=it.second.second; + output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << ", " + << var+1-ModelBlock->Block_List[j].Nb_Recursives << ") = "; + writeChainRuleDerivative(output, eqr, varr, k, oMatlabStaticModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr)) + << " " << varr+1 << ", equation=" << eqr+1 << endl; + } + break; + default: + break; + } + output.close(); + } + } + +void +StaticDllModel::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, map_idx_type map_idx) const + { + struct Uff_l + { + int u, var, lag; + Uff_l *pNext; + }; + + struct Uff + { + Uff_l *Ufl, *Ufl_First; + }; + + int i,j,k,v; + string tmp_s; + ostringstream tmp_output; + ofstream code_file; + NodeID lhs=NULL, rhs=NULL; + BinaryOpNode *eq_node; + Uff Uf[symbol_table.endo_nbr()]; + map reference_count; + vector feedback_variables; + bool file_open=false; + string main_name=file_name; + main_name+=".cod"; + code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate ); + if (!code_file.is_open()) + { + cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + exit(EXIT_FAILURE); + } + //Temporary variables declaration + code_file.write(&FDIMST, sizeof(FDIMST)); + k=temporary_terms.size(); + code_file.write(reinterpret_cast(&k),sizeof(k)); + + for (j = 0; j < ModelBlock->Size ;j++) + { + feedback_variables.clear(); + if (j>0) + code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); + code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK)); + v=ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives; + code_file.write(reinterpret_cast(&v),sizeof(v)); + v=ModelBlock->Block_List[j].Simulation_Type; + code_file.write(reinterpret_cast(&v),sizeof(v)); + int count_u; + for (i=ModelBlock->Block_List[j].Nb_Recursives; i < ModelBlock->Block_List[j].Size;i++) + { + code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i])); + code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i])); + code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i])); + } + if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE) + { + int u_count_int=0; + Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open); + code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].is_linear),sizeof(ModelBlock->Block_List[j].is_linear)); + v = u_count_int ; + code_file.write(reinterpret_cast(&v),sizeof(v)); + v=symbol_table.endo_nbr(); + code_file.write(reinterpret_cast(&v),sizeof(v)); + v=block_triangular.ModelBlock->Block_List[j].Max_Lag; + code_file.write(reinterpret_cast(&v),sizeof(v)); + v=block_triangular.ModelBlock->Block_List[j].Max_Lead; + code_file.write(reinterpret_cast(&v),sizeof(v)); + + v=u_count_int; + code_file.write(reinterpret_cast(&v),sizeof(v)); + file_open=true; + } + // The equations + //cout << block_triangular.BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " j=" << j << endl; + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + //The Temporary terms + //cout << "equation = " << ModelBlock->Block_List[j].Equation[i] << " variable = " << ModelBlock->Block_List[j].Variable[i] << " r[" << i << "] " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl; + temporary_terms_type tt2; + tt2.clear(); + 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, tt2, map_idx, false); + code_file.write(&FSTPST, sizeof(FSTPST)); + map_idx_type::const_iterator ii=map_idx.find((*it)->idx); + v=(int)ii->second; + code_file.write(reinterpret_cast(&v), sizeof(v)); + // Insert current node into tt2 + tt2.insert(*it); + } + switch (ModelBlock->Block_List[j].Simulation_Type) + { +evaluation: + case EVALUATE_BACKWARD: + case EVALUATE_FORWARD: + if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE) + { + eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + rhs->compile(code_file, false, temporary_terms, map_idx, false); + lhs->compile(code_file, true, temporary_terms, map_idx, false); + } + else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S) + { + eq_node = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[i]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + rhs->compile(code_file, false, temporary_terms, map_idx, false); + lhs->compile(code_file, true, temporary_terms, map_idx, false); + } + break; + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + if (iBlock_List[j].Nb_Recursives) + goto evaluation; + feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]); + v=ModelBlock->Block_List[j].Equation[i]; + Uf[v].Ufl=NULL; + goto end; + default: +end: + eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + lhs->compile(code_file, false, temporary_terms, map_idx, false); + rhs->compile(code_file, false, temporary_terms, map_idx, false); + code_file.write(&FBINARY, sizeof(FBINARY)); + int v=oMinus; + code_file.write(reinterpret_cast(&v),sizeof(v)); + code_file.write(&FSTPR, sizeof(FSTPR)); + v = i - ModelBlock->Block_List[j].Nb_Recursives; + code_file.write(reinterpret_cast(&v), sizeof(v)); + } + } + 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) + { + 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, map_idx); + code_file.write(&FSTPG, sizeof(FSTPG)); + v=0; + code_file.write(reinterpret_cast(&v), sizeof(v)); + break; + + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + count_u = feedback_variables.size(); + for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++) + { + pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i); + k=it.first.first; + int eq=it.first.second.first; + int var=it.first.second.second; + int eqr=it.second.first; + int varr=it.second.second; + int v=ModelBlock->Block_List[j].Equation[eq]; + if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives) + { + 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=count_u; + Uf[v].Ufl->var=varr; + Uf[v].Ufl->lag=k; + compileChainRuleDerivative(code_file, eqr, varr, k, map_idx); + code_file.write(&FSTPSU, sizeof(FSTPSU)); + code_file.write(reinterpret_cast(&count_u), sizeof(count_u)); + count_u++; + } + } + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + if(i>=ModelBlock->Block_List[j].Nb_Recursives) + { + code_file.write(&FLDR, sizeof(FLDR)); + v = i-ModelBlock->Block_List[j].Nb_Recursives; + code_file.write(reinterpret_cast(&v), sizeof(v)); + code_file.write(&FLDZ, sizeof(FLDZ)); + 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(&FLDSU, sizeof(FLDSU)); + code_file.write(reinterpret_cast(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); + code_file.write(&FLDSV, sizeof(FLDSV)); + char vc=eEndogenous; + code_file.write(reinterpret_cast(&vc), sizeof(vc)); + int v1=Uf[v].Ufl->var; + code_file.write(reinterpret_cast(&v1), sizeof(v1)); + code_file.write(&FBINARY, sizeof(FBINARY)); + v1=oTimes; + code_file.write(reinterpret_cast(&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(&v), sizeof(v)); + code_file.write(&FSTPSU, sizeof(FSTPSU)); + v = i - ModelBlock->Block_List[j].Nb_Recursives; + code_file.write(reinterpret_cast(&v), sizeof(v)); + } + } + break; + default: + break; + } + } + } + code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); + code_file.write(&FEND, sizeof(FEND)); + code_file.close(); + } + + + +void +StaticDllModel::Write_Inf_To_Bin_File(const string &static_basename, const string &bin_basename, const int &num, + int &u_count_int, bool &file_open) const + { + int j; + std::ofstream SaveCode; + if (file_open) + SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate ); + else + SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::binary); + if (!SaveCode.is_open()) + { + cout << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing\n"; + exit(EXIT_FAILURE); + } + u_count_int=0; + int Size = block_triangular.ModelBlock->Block_List[num].Size - block_triangular.ModelBlock->Block_List[num].Nb_Recursives; + for(int i=0; i<(int)block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->size();i++) + { + //Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag)))); + pair< pair >, pair > it = block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->at(i); + int k=it.first.first; + int eq=it.first.second.first; + + int var_init=it.first.second.second; + /*int eqr=it.second.first; + int varr=it.second.second;*/ + if(eq>=block_triangular.ModelBlock->Block_List[num].Nb_Recursives and var_init>=block_triangular.ModelBlock->Block_List[num].Nb_Recursives) + { + int v=eq-block_triangular.ModelBlock->Block_List[num].Nb_Recursives; + SaveCode.write(reinterpret_cast(&v), sizeof(v)); + int var=it.first.second.second-block_triangular.ModelBlock->Block_List[num].Nb_Recursives + k * Size; + SaveCode.write(reinterpret_cast(&var), sizeof(var)); + SaveCode.write(reinterpret_cast(&k), sizeof(k)); + int u = u_count_int + Size; + SaveCode.write(reinterpret_cast(&u), sizeof(u)); + //cout << "eq=" << v << ", var=" << var << ", lag=" << k << " u=" << u << "\n"; + u_count_int++; + } + } + /*cout << "u_count_int=" << u_count_int << endl; + cout << "block_triangular.ModelBlock->Block_List[" << num << "].Nb_Recursives=" << block_triangular.ModelBlock->Block_List[num].Nb_Recursives << " block_triangular.ModelBlock->Block_List[" << num << "].Size=" << block_triangular.ModelBlock->Block_List[num].Size << endl;*/ + for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;jBlock_List[num].Size;j++) + { + int varr=block_triangular.ModelBlock->Block_List[num].Variable[j]; + //cout << "j=" << j << " varr=" << varr << "\n"; + SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); + } + for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;jBlock_List[num].Size;j++) + { + int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j]; + SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); + } + SaveCode.close(); + } + +void +StaticDllModel::writeSparseStaticMFile(const string &static_basename, const string &basename, const int mode) const + { + string sp; + ofstream mStaticModelFile; + ostringstream tmp, tmp1, tmp_eq; + int prev_Simulation_Type; + //SymbolicGaussElimination SGE; + bool OK; + chdir(basename.c_str()); + string 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"; + + int i, k; + + mStaticModelFile << "function [varargout] = " << static_basename << "(varargin)\n"; + mStaticModelFile << " global oo_ options_ M_ ;\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) + 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, oMatlabStaticModel, temporary_terms); + tmp_output << "=T_init;\n"; + } + if (tmp_output.str().length()>0) + mStaticModelFile << tmp_output.str(); + + mStaticModelFile << " y_size=M_.endo_nbr;\n"; + mStaticModelFile << " %it is the deterministic simulation of the block decomposed static model\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 oo_.exo_det_steady_state];\n"; + + prev_Simulation_Type=-1; + mStaticModelFile << " params=M_.params;\n"; + mStaticModelFile << " 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 ((k == EVALUATE_FORWARD) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + mStaticModelFile << " oo_.deterministic_simulation.status = 1;\n"; + mStaticModelFile << " oo_.deterministic_simulation.error = 0;\n"; + mStaticModelFile << " oo_.deterministic_simulation.iterations = 0;\n"; + mStaticModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mStaticModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mStaticModelFile << " else\n"; + mStaticModelFile << " blck_num = 1;\n"; + mStaticModelFile << " end;\n"; + mStaticModelFile << " oo_.deterministic_simulation.block(blck_num).status = 1;\n"; + mStaticModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n"; + mStaticModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n"; + //mStaticModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n"; + mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x, params);\n"; + mStaticModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n"; + mStaticModelFile << " if(isnan(tmp) | isinf(tmp))\n"; + mStaticModelFile << " disp(['Inf or Nan value during the evaluation of block " << i <<"']);\n"; + mStaticModelFile << " return;\n"; + mStaticModelFile << " end;\n"; + } + else if ((k == EVALUATE_BACKWARD) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + mStaticModelFile << " oo_.deterministic_simulation.status = 1;\n"; + mStaticModelFile << " oo_.deterministic_simulation.error = 0;\n"; + mStaticModelFile << " oo_.deterministic_simulation.iterations = 0;\n"; + mStaticModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mStaticModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mStaticModelFile << " else\n"; + mStaticModelFile << " blck_num = 1;\n"; + mStaticModelFile << " end;\n"; + mStaticModelFile << " oo_.deterministic_simulation.block(blck_num).status = 1;\n"; + mStaticModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n"; + mStaticModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n"; + mStaticModelFile << " " << static_basename << "_" << i + 1 << "(y, x, params);\n"; + mStaticModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n"; + mStaticModelFile << " if(isnan(tmp) | isinf(tmp))\n"; + mStaticModelFile << " disp(['Inf or Nan value during the evaluation of block " << i <<"']);\n"; + mStaticModelFile << " return;\n"; + mStaticModelFile << " end;\n"; + } + else if ((k == SOLVE_FORWARD_COMPLETE || k == SOLVE_FORWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + /*mStaticModelFile << " g1=0;\n"; + mStaticModelFile << " r=0;\n";*/ + tmp.str(""); + for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_List[i].Size;ik++) + { + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + } + mStaticModelFile << " 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++)*/ + int nze=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[0].size; + mStaticModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mStaticModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mStaticModelFile << " else\n"; + mStaticModelFile << " blck_num = 1;\n"; + mStaticModelFile << " end;\n"; + mStaticModelFile << " y = solve_one_boundary('" << static_basename << "_" << i + 1 << "'" << + ", y, x, params, y_index, " << nze << + ", 1, " << block_triangular.ModelBlock->Block_List[i].is_linear << + ", blck_num, 0, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 0, 0);\n"; + mStaticModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n"; + mStaticModelFile << " if(isnan(tmp) | isinf(tmp))\n"; + mStaticModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n"; + mStaticModelFile << " return;\n"; + mStaticModelFile << " end;\n"; + } + else if ((k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + /*mStaticModelFile << " g1=0;\n"; + mStaticModelFile << " r=0;\n";*/ + tmp.str(""); + for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_List[i].Size;ik++) + { + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + } + mStaticModelFile << " 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++)*/ + int nze=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[0].size; + mStaticModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; + mStaticModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; + mStaticModelFile << " else\n"; + mStaticModelFile << " blck_num = 1;\n"; + mStaticModelFile << " end;\n"; + mStaticModelFile << " y = solve_one_boundary('" << static_basename << "_" << i + 1 << "'" << + ", y, x, params, y_index, " << nze << + ", 0, " << 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 << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n"; + mStaticModelFile << " if(isnan(tmp) | isinf(tmp))\n"; + mStaticModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n"; + mStaticModelFile << " return;\n"; + mStaticModelFile << " end;\n"; + } + } + mStaticModelFile << " oo_.steady_state = y';\n"; + mStaticModelFile << "return;\n"; + + mStaticModelFile.close(); + + writeModelEquationsOrdered_M( block_triangular.ModelBlock, static_basename); + + chdir(".."); + } + + +void +StaticDllModel::writeOutput(ostream &output, const string &basename) const +{ + output << "options_.model_mode = " << mode << ";" << endl; + + // Erase possible remnants of previous runs + if (mode != eStandardMode) + output << "delete('" << basename << "_static.m');" << endl; + if (mode != eDLLMode) + output << "erase_compiled_function('" + basename + "_static');" << endl; + + // Special setup for DLL or Sparse modes + if (mode == eDLLMode) + output << "mex -O LDFLAGS='-pthread -shared -Wl,--no-undefined' " << basename << "_static.c" << endl; + if (mode == eSparseMode || mode == eSparseDLLMode) + output << "addpath " << basename << ";" << endl; +} + +void +StaticDllModel::writeOutputPostComputing(ostream &output, const string &basename) const +{ + if (mode == eSparseMode || mode == eSparseDLLMode) + output << "rmpath " << basename << ";" << endl; +} + +void +StaticDllModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic) +{ + 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 (getTypeByDerivID(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(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ") [" << getSymbIDByDerivID(it->first.second) << "] !" << endl; + Id->writeOutput(cout, oMatlabStaticModelSparse, temporary_terms); + cout << "\n"; + cerr << "StaticDllModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ")!" << endl; + } + int eq=it->first.first; + int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second); + int k1 = getLagByDerivID(it->first.second); + if (a_variable_lag!=k1) + { + IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous); + a_variable_lag=k1; + } + if (k1==0 or !dynamic) + { + 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++; + } + } + } + //Get ride of the elements of the incidence matrix equal to Zero + IM=block_triangular.incidencematrix.Get_IM(0, eEndogenous); + for (int i=0;i0) + { + 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 +StaticDllModel::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;iBlock_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,getDerivID(symbol_table.getID(eEndogenous, var),0))); + if (it!= first_derivatives.end()) + { + NodeID Id = it->second; + set > endogenous; + Id->collectEndogenous(endogenous); + if (endogenous.size() > 0) + { + for (l=0;lBlock_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;iBlock_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,getDerivID(symbol_table.getID(eEndogenous, var),k1))); + NodeID Id = it->second; + if (it!= first_derivatives.end()) + { + set > endogenous; + Id->collectEndogenous(endogenous); + if (endogenous.size() > 0) + { + for (l=0;lBlock_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; + } +} + + +map >, NodeID> +StaticDllModel::collect_first_order_derivatives_endogenous() +{ + map >, NodeID> endo_derivatives; + for (first_derivatives_type::iterator it2 = first_derivatives.begin(); + it2 != first_derivatives.end(); it2++) + { + if (getTypeByDerivID(it2->first.second)==eEndogenous) + { + int eq = it2->first.first; + int var=symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second)); + int lag=getLagByDerivID(it2->first.second); + //if (lag==0) + endo_derivatives[make_pair(eq, make_pair(var, lag))] = it2->second; + } + } + return endo_derivatives; +} + + + +void +StaticDllModel::computingPass(const eval_context_type &eval_context, bool no_tmp_terms) +{ + + // Computes static jacobian columns + computeStatJacobianCols(); + + // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic + set vars; + for (deriv_id_table_t::const_iterator it = deriv_id_table.begin(); + it != deriv_id_table.end(); it++) + { + SymbolType type = symbol_table.getType(it->first.first); + if (type == eEndogenous) + vars.insert(it->second); + } + + // Launch computations + cout << "Computing static model derivatives:" << endl + << " - order 1" << endl; + computeJacobian(vars); + //cout << "mode=" << mode << " eSparseDLLMode=" << eSparseDLLMode << " eSparseMode=" << eSparseMode << "\n"; + if (mode == eSparseDLLMode || mode == eSparseMode) + { + BuildIncidenceMatrix(); + + jacob_map j_m; + evaluateJacobian(eval_context, &j_m, true); + + + if (block_triangular.bt_verbose) + { + cout << "The gross incidence matrix \n"; + block_triangular.incidencematrix.Print_IM(eEndogenous); + } + t_etype equation_simulation_type; + map >, NodeID> first_order_endo_derivatives = collect_first_order_derivatives_endogenous(); + + block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations, equation_simulation_type, first_order_endo_derivatives, mfs, cutoff); + /*for (int j = 0;j < block_triangular.ModelBlock->Size;j++) + { + for (int i = 0;i < block_triangular.ModelBlock->Block_List[j].Size;i++) + { + if (iBlock_List[j].Nb_Recursives ) + cout << "block=" << j << " R i=" << i << " equation=" << block_triangular.ModelBlock->Block_List[j].Equation[i]+1 << " variable=" << block_triangular.ModelBlock->Block_List[j].Variable[i]+1 << endl; + else + cout << "block=" << j << " S i=" << i << " equation=" << block_triangular.ModelBlock->Block_List[j].Equation[i]+1 << " variable=" << block_triangular.ModelBlock->Block_List[j].Variable[i]+1 << endl; + } + }*/ + + BlockLinear(block_triangular.ModelBlock); + + computeChainRuleJacobian(block_triangular.ModelBlock); + + if (!no_tmp_terms) + computeTemporaryTermsOrdered(block_triangular.ModelBlock); + + } + else + if (!no_tmp_terms) + computeTemporaryTerms(mode == eStandardMode); +} + +void +StaticDllModel::writeStaticFile(const string &basename) const + { + int r; + switch (mode) + { + case eStandardMode: + break; + case eSparseMode: +#ifdef _WIN32 + r = mkdir(basename.c_str()); +#else + r = mkdir(basename.c_str(), 0777); +#endif + if (r < 0 && errno != EEXIST) + { + perror("ERROR"); + exit(EXIT_FAILURE); + } + writeSparseStaticMFile(basename + "_static", 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: + break; + case eSparseDLLMode: + // create a directory to store all the files +#ifdef _WIN32 + r = mkdir(basename.c_str()); +#else + r = mkdir(basename.c_str(), 0777); +#endif + if (r < 0 && errno != EEXIST) + { + perror("ERROR"); + exit(EXIT_FAILURE); + } + writeModelEquationsCodeOrdered(basename + "_static", block_triangular.ModelBlock, basename, 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; + } + } + + +int +StaticDllModel::computeDerivID(int symb_id, int lag) +{ + // Check if static variable already has a deriv_id + pair key = make_pair(symb_id, lag); + deriv_id_table_t::const_iterator it = deriv_id_table.find(key); + if (it != deriv_id_table.end()) + return it->second; + + // Create a new deriv_id + int deriv_id = deriv_id_table.size(); + + deriv_id_table[key] = deriv_id; + inv_deriv_id_table.push_back(key); + + SymbolType type = symbol_table.getType(symb_id); + + if (type == eEndogenous) + dynJacobianColsNbr++; + + return deriv_id; +} + +SymbolType +StaticDllModel::getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException) +{ + return symbol_table.getType(getSymbIDByDerivID(deriv_id)); +} + +int +StaticDllModel::getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException) +{ + if (deriv_id < 0 || deriv_id >= (int) inv_deriv_id_table.size()) + throw UnknownDerivIDException(); + + return inv_deriv_id_table[deriv_id].second; +} + +int +StaticDllModel::getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException) +{ + if (deriv_id < 0 || deriv_id >= (int) inv_deriv_id_table.size()) + throw UnknownDerivIDException(); + + return inv_deriv_id_table[deriv_id].first; +} + +int +StaticDllModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException) +{ + deriv_id_table_t::const_iterator it = deriv_id_table.find(make_pair(symb_id, lag)); + if (it == deriv_id_table.end()) + throw UnknownDerivIDException(); + else + return it->second; +} + +void +StaticDllModel::computeStatJacobianCols() +{ + /* Sort the static endogenous variables by lexicographic order over (lag, type_specific_symbol_id) + and fill the static columns for exogenous and exogenous deterministic */ + map, int> ordered_dyn_endo; + + for (deriv_id_table_t::const_iterator it = deriv_id_table.begin(); + it != deriv_id_table.end(); it++) + { + const int &symb_id = it->first.first; + const int &lag = it->first.second; + const int &deriv_id = it->second; + SymbolType type = symbol_table.getType(symb_id); + int tsid = symbol_table.getTypeSpecificID(symb_id); + + switch (type) + { + case eEndogenous: + ordered_dyn_endo[make_pair(lag, tsid)] = deriv_id; + break; + case eExogenous: + // At this point, dynJacobianColsNbr contains the number of static endogenous + break; + case eExogenousDet: + // At this point, dynJacobianColsNbr contains the number of static endogenous + break; + case eParameter: + // We don't assign a static jacobian column to parameters + break; + default: + // Shut up GCC + cerr << "StaticDllModel::computeStatJacobianCols: impossible case" << endl; + exit(EXIT_FAILURE); + } + } + + // Fill in static jacobian columns for endogenous + int sorted_id = 0; + for (map, int>::const_iterator it = ordered_dyn_endo.begin(); + it != ordered_dyn_endo.end(); it++) + dyn_jacobian_cols_table[it->second] = sorted_id++; + +} + +int +StaticDllModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException) +{ + map::const_iterator it = dyn_jacobian_cols_table.find(deriv_id); + if (it == dyn_jacobian_cols_table.end()) + throw UnknownDerivIDException(); + else + return it->second; +} + + +void +StaticDllModel::computeChainRuleJacobian(Model_Block *ModelBlock) +{ + map recursive_variables; + first_chain_rule_derivatives.clear(); + for(int blck = 0; blckSize; blck++) + { + recursive_variables.clear(); + if (ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) + { + ModelBlock->Block_List[blck].Chain_Rule_Derivatives->clear(); + for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++) + { + if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S) + recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i]; + else + recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]]; + } + map >, pair >, int> Derivatives = block_triangular.get_Derivatives(ModelBlock, blck); + + map >, pair >, int>::const_iterator it = Derivatives.begin(); + //#pragma omp parallel for shared(it, blck) + for(int i=0; i<(int)Derivatives.size(); i++) + { + int Deriv_type = it->second; + pair >, pair > it_l(it->first); + it++; + int lag = it_l.first.first; + int eq = it_l.first.second.first; + int var = it_l.first.second.second; + int eqr = it_l.second.first; + int varr = it_l.second.second; + if(Deriv_type == 0) + { + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))]; + } + else if (Deriv_type == 1) + { + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables); + } + else if (Deriv_type == 2) + { + if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S && eqBlock_List[blck].Nb_Recursives) + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables); + else + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables); + } + ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr))); + } + } + else if( ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_SIMPLE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_SIMPLE + or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_COMPLETE) + { + ModelBlock->Block_List[blck].Chain_Rule_Derivatives->clear(); + for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++) + { + if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S) + recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i]; + else + recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]]; + } + for(int eq = ModelBlock->Block_List[blck].Nb_Recursives; eq < ModelBlock->Block_List[blck].Size; eq++) + { + int eqr = ModelBlock->Block_List[blck].Equation[eq]; + for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++) + { + int varr = ModelBlock->Block_List[blck].Variable[var]; + NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables); + if (d1 == Zero) + continue; + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1; + ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(0, make_pair(eq, var)), make_pair(eqr, varr))); + } + } + } + } +} + + + +void +StaticDllModel::computeParamsDerivatives() +{ + for (deriv_id_table_t::const_iterator it = deriv_id_table.begin(); + it != deriv_id_table.end(); it++) + { + if (symbol_table.getType(it->first.first) != eParameter) + continue; + + int param = it->second; + + for (first_derivatives_type::const_iterator it2 = first_derivatives.begin(); + it2 != first_derivatives.end(); it2++) + { + int eq = it2->first.first; + int var = it2->first.second; + NodeID d1 = it2->second; + + NodeID d2 = d1->getDerivative(param); + if (d2 == Zero) + continue; + params_derivatives[make_pair(eq, make_pair(var, param))] = d2; + } + } +} + +void +StaticDllModel::computeParamsDerivativesTemporaryTerms() +{ + map reference_count; + params_derivs_temporary_terms.clear(); + + for (second_derivatives_type::iterator it = params_derivatives.begin(); + it != params_derivatives.end(); it++) + it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true); +} + +void +StaticDllModel::writeParamsDerivativesFile(const string &basename) const + { + if (!params_derivatives.size()) + return; + + string filename = basename + "_params_derivs.m"; + + ofstream paramsDerivsFile; + paramsDerivsFile.open(filename.c_str(), ios::out | ios::binary); + if (!paramsDerivsFile.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + paramsDerivsFile << "function gp = " << basename << "_params_derivs(y, x, params, it_)" << endl + << "%" << endl + << "% Warning : this file is generated automatically by Dynare" << endl + << "% from model file (.mod)" << endl << endl; + + + writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, oMatlabStaticModel); + + paramsDerivsFile << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " + << symbol_table.param_nbr() << ");" << endl; + + for (second_derivatives_type::const_iterator it = params_derivatives.begin(); + it != params_derivatives.end(); it++) + { + int eq = it->first.first; + int var = it->first.second.first; + int param = it->first.second.second; + NodeID d2 = it->second; + + int var_col = getDynJacobianCol(var) + 1; + int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; + + paramsDerivsFile << "gp(" << eq+1 << ", " << var_col << ", " << param_col << ") = "; + d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms); + paramsDerivsFile << ";" << endl; + } + + paramsDerivsFile.close(); + } + + + +void +StaticDllModel::writeChainRuleDerivative(ostream &output, int eqr, int varr, int lag, + ExprNodeOutputType output_type, + const temporary_terms_type &temporary_terms) const +{ + map >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag))); + if (it != first_chain_rule_derivatives.end()) + (it->second)->writeOutput(output, output_type, temporary_terms); + else + output << 0; +} + + +void +StaticDllModel::writeLatexFile(const string &basename) const + { + writeLatexModelFile(basename + "_static.tex", oLatexStaticModel); + } + +void +StaticDllModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const +{ + output << LEFT_ARRAY_SUBSCRIPT(output_type); + if (IS_MATLAB(output_type)) + output << eq_nb + 1 << ", " << col_nb + 1; + else + output << eq_nb + col_nb * equations.size(); + output << RIGHT_ARRAY_SUBSCRIPT(output_type); +} + +void +StaticDllModel::hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const +{ + output << LEFT_ARRAY_SUBSCRIPT(output_type); + if (IS_MATLAB(output_type)) + output << row_nb + 1 << ", " << col_nb + 1; + else + output << row_nb + col_nb * NNZDerivatives[1]; + output << RIGHT_ARRAY_SUBSCRIPT(output_type); +} + + diff --git a/StaticDllModel.hh b/StaticDllModel.hh new file mode 100644 index 00000000..dd067710 --- /dev/null +++ b/StaticDllModel.hh @@ -0,0 +1,192 @@ +/* + * 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 . + */ + +#ifndef _STATICDLLMODEL_HH +#define _STATICDLLMODEL_HH + +using namespace std; + +#include + +#include "StaticModel.hh" +#include "BlockTriangular.hh" + +//! Stores a static model +class StaticDllModel : public ModelTree +{ +public: + //! The modes in which StaticDllModel can work + enum mode_t + { + eStandardMode, //!< Standard mode (static file in Matlab) + eSparseMode, //!< Sparse mode (static file in Matlab with block decomposition) + eDLLMode, //!< DLL mode (static file in C) + eSparseDLLMode //!< Sparse DLL mode (static file in C with block decomposition plus a binary file) + }; +private: + typedef map, int> deriv_id_table_t; + //! Maps a pair (symbol_id, lag) to a deriv ID + deriv_id_table_t deriv_id_table; + //! Maps a deriv ID to a pair (symbol_id, lag) + vector > inv_deriv_id_table; + + //! Maps a deriv_id to the column index of the static Jacobian + /*! Contains only endogenous, exogenous and exogenous deterministic */ + map dyn_jacobian_cols_table; + + //! Maximum lag and lead over all types of variables (positive values) + /*! Set by computeDerivID() */ + int max_lag, max_lead; + //! Maximum lag and lead over endogenous variables (positive values) + /*! Set by computeDerivID() */ + int max_endo_lag, max_endo_lead; + //! Maximum lag and lead over exogenous variables (positive values) + /*! Set by computeDerivID() */ + int max_exo_lag, max_exo_lead; + //! Maximum lag and lead over deterministic exogenous variables (positive values) + /*! Set by computeDerivID() */ + int max_exo_det_lag, max_exo_det_lead; + + //! Number of columns of static jacobian + /*! Set by computeDerivID() and computeDynJacobianCols() */ + int dynJacobianColsNbr; + + //! Derivatives of the jacobian w.r. to parameters + /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter. + Only non-null derivatives are stored in the map. + Variable and parameter indices are those of the getDerivID() method. + */ + second_derivatives_type params_derivatives; + + //! Temporary terms for the file containing parameters dervicatives + temporary_terms_type params_derivs_temporary_terms; + + typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_type; + first_chain_rule_derivatives_type first_chain_rule_derivatives; + + + //! Writes static model file (Matlab version) + void writeStaticMFile(const string &static_basename) const; + //! Writes static model file (C version) + /*! \todo add third derivatives handling */ + void writeStaticCFile(const string &static_basename) const; + //! Writes static model file when SparseDLL option is on + void writeSparseStaticMFile(const string &static_basename, const string &basename, const int mode) const; + //! Writes the Block reordred structure of the model in M output + void writeModelEquationsOrdered_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, 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, bool dynamic); + void BlockLinear(Model_Block *ModelBlock); + map_idx_type map_idx; + //! Build The incidence matrix form the modeltree + void BuildIncidenceMatrix(); + + void computeTemporaryTermsOrdered(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, map_idx_type &map_idx) const; + //! Write chain rule derivative code of an equation w.r. to a variable + void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const; + + virtual int computeDerivID(int symb_id, int lag); + //! Get the type corresponding to a derivation ID + SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException); + //! Get the lag corresponding to a derivation ID + int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException); + //! Get the symbol ID corresponding to a derivation ID + int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException); + //! Compute the column indices of the static Jacobian + void computeStatJacobianCols(); + //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables + void computeChainRuleJacobian(Model_Block *ModelBlock); + //! Computes derivatives of the Jacobian w.r. to parameters + void computeParamsDerivatives(); + //! Computes temporary terms for the file containing parameters derivatives + void computeParamsDerivativesTemporaryTerms(); + //! Collect only the first derivatives + map >, NodeID> collect_first_order_derivatives_endogenous(); + + //! Helper for writing the Jacobian elements in MATLAB and C + /*! Writes either (i+1,j+1) or [i+j*no_eq] */ + void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const; + + //! Helper for writing the sparse Hessian elements in MATLAB and C + /*! Writes either (i+1,j+1) or [i+j*NNZDerivatives[1]] */ + void hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const; + + //! Write chain rule derivative of a recursive equation w.r. to a variable + void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; + + +public: + StaticDllModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); + //! Mode in which the ModelTree is supposed to work (Matlab, DLL or SparseDLL) + mode_t mode; + //! Adds a variable node + /*! This implementation allows for non-zero lag */ + virtual NodeID AddVariable(const string &name, int lag = 0); + //! 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 and simul_NG from simulate.cc) + double markowitz; + //! Compute the minimum feedback set in the static model: + /*! 0 : all endogenous variables are considered as feedback variables + 1 : the variables belonging to a non linear equation are considered as feedback variables + 2 : the variables belonging to a non normalizable non linear equation are considered as feedback variables + default value = 0 */ + int mfs; + //! the file containing the model and the derivatives code + ofstream code_file; + //! Execute computations (variable sorting + derivation) + /*! + \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed) + \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true) + \param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true) + \param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed (implies jacobianExo = true) + \param eval_context evaluation context for normalization + \param no_tmp_terms if true, no temporary terms will be computed in the static 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 string &basename) const; + //! Write statements to be added to the main M-file, after computational tasks + void writeOutputPostComputing(ostream &output, 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 &static_basename, const string &bin_basename, + const int &num, int &u_count_int, bool &file_open) const; + //! Writes static model file + void writeStaticFile(const string &basename) const; + //! Writes file containing parameters derivatives + void writeParamsDerivativesFile(const string &basename) const; + + //! Writes LaTeX file with the equations of the static model + void writeLatexFile(const string &basename) const; + + virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); + virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException); +}; + +#endif