From 1a058fea4f1a08c060a9571fd33b2f71dc365766 Mon Sep 17 00:00:00 2001 From: ferhat Date: Fri, 15 May 2009 22:41:51 +0000 Subject: [PATCH] - Bugs correction in the new block decomposition (incorporating the feedback variables) - First draft of DynamicModel.cc files with feedback variables. TODO : - reduction of the Jacobian matrix - symbolic normalization of equations - application to the binary code evaluation (simulate.dll). git-svn-id: https://www.dynare.org/svn/dynare/trunk@2678 ac1d8469-bf42-47a9-8791-bf33cf982152 --- preprocessor/BlockTriangular.cc | 209 +++++++++++++++++++++----------- preprocessor/BlockTriangular.hh | 31 +++-- preprocessor/CodeInterpreter.hh | 11 ++ preprocessor/DynamicModel.cc | 75 +++++++++--- preprocessor/DynamicModel.hh | 3 +- preprocessor/DynareMain.cc | 4 +- preprocessor/ExprNode.hh | 5 +- 7 files changed, 239 insertions(+), 99 deletions(-) diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc index 16067918c..e52c7a6a9 100644 --- a/preprocessor/BlockTriangular.cc +++ b/preprocessor/BlockTriangular.cc @@ -126,7 +126,6 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog Vertices n to 2*n-1 are for equations (using equation no.) */ BipartiteGraph g(2 * n); - // Fill in the graph for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) @@ -155,24 +154,23 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog } return false; } + vector Index_Equ_IM_tmp(Index_Equ_IM); + bool *SIM; + SIM = (bool*)malloc(equation_number*equation_number*sizeof(bool)); + memcpy(SIM, IM, equation_number*equation_number*sizeof(bool)); for (int i = 0; i < n; i++) { - int j = Index_Equ_IM[mate_map[i]-n+prologue]; - Index_Equ_IM[mate_map[i]-n+prologue] = Index_Equ_IM[i+prologue]; - Index_Equ_IM[i+prologue] = j; + Index_Equ_IM[i + prologue] = Index_Equ_IM_tmp[mate_map[i] - n + prologue]; for (int k=0; k &Index_Equ_IM, vector &Index_Var_IM, vector > &blocks, bool verbose_) const +BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector &Index_Equ_IM, vector &Index_Var_IM, vector > &blocks, t_etype &Equation_Type, bool verbose_) const { int n = nb_var - prologue - epilogue; bool *AMp; @@ -180,20 +178,19 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo //transforms the incidence matrix of the complet model into an adjancent matrix of the non-recursive part of the model for (int i=prologue; i component(num_vertices(G2)), discover_time(num_vertices(G2)); - //write_graphviz(cout, G2); - int num = strong_components(G2, &component[0]); - blocks = vector >(num, make_pair(0,0)); + blocks = vector >(num , make_pair(0,0)); //This vector contains for each block: // - first set = equations belonging to the block, @@ -207,24 +204,33 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo components_set[component[i]].first.insert(i); } - //For each block, the minimum set of feedback variable is computed - // and the non-feedback variables are reordered to get - // a sub-recursive block without feedback variables + vector tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM); int order = prologue; bool *SIM; SIM = (bool*)malloc(nb_var*nb_var*sizeof(bool)); memcpy(SIM, IM, nb_var*nb_var*sizeof(bool)); + + //Add a loop on vertices which could not be normalized => force those vvertices to belong to the feedback set + for(int i=0; i feed_back_vertices; + //Print(G); AdjacencyList_type G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G); property_map::type v_index = get(vertex_index, G); components_set[i].second.first = feed_back_vertices; blocks[i].second = feed_back_vertices.size(); vector Reordered_Vertice; Reordered_Vertice = Reorder_the_recursive_variables(G, feed_back_vertices); + //First we have the recursive equations conditional on feedback variables for (vector::iterator its=Reordered_Vertice.begin();its != Reordered_Vertice.end(); its++) { Index_Equ_IM[order] = tmp_Index_Equ_IM[*its+prologue]; @@ -232,7 +238,7 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo order++; } components_set[i].second.second = Reordered_Vertice; - + //Second we have the equations related to the feedback variables for (set::iterator its=feed_back_vertices.begin();its != feed_back_vertices.end(); its++) { Index_Equ_IM[order] = tmp_Index_Equ_IM[v_index[vertex(*its, G)]+prologue]; @@ -240,10 +246,12 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo order++; } } + free(AMp); + free(SIM); } void -BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock) +BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock, t_etype &Equation_Type, int recurs_Size) { int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li; int *tmp_size, *tmp_size_other_endo, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_other_endo, *tmp_exo, tmp_nb_other_endo, tmp_nb_exo, nb_lead_lag_endo; @@ -252,14 +260,24 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block bool *IM, OK; int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo; + + cout << "Allocate Block=" << count_Block << " recurs_Size=" << recurs_Size << "\n"; ModelBlock->Periods = periods; ModelBlock->Block_List[count_Block].is_linear=true; ModelBlock->Block_List[count_Block].Size = size; ModelBlock->Block_List[count_Block].Type = type; + ModelBlock->Block_List[count_Block].Nb_Recursives = recurs_Size; ModelBlock->Block_List[count_Block].Temporary_InUse=new temporary_terms_inuse_type (); ModelBlock->Block_List[count_Block].Temporary_InUse->clear(); ModelBlock->Block_List[count_Block].Simulation_Type = SimType; ModelBlock->Block_List[count_Block].Equation = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); + ModelBlock->Block_List[count_Block].Equation_Type = (EquationType*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(EquationType)); + ModelBlock->Block_List[count_Block].Equation_Type_Var = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); + /*cout << "ModelBlock->Block_List[" << count_Block << "].Size = " << ModelBlock->Block_List[count_Block].Size << " E_UNKNOWN=" << E_UNKNOWN << "\n"; + t_etype eType(size); + cout << "eType[0].first=" << eType[0].first << " eType[0].second=" << eType[0].second << " eType.size()=" << eType.size() << "\n"; + ModelBlock->Block_List[count_Block].Equation_Type(eType); + cout << "after that\n";*/ ModelBlock->Block_List[count_Block].Variable = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation = (temporary_terms_type**)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(temporary_terms_type)); ModelBlock->Block_List[count_Block].Own_Derivative = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); @@ -289,6 +307,8 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]->clear(); ModelBlock->Block_List[count_Block].Equation[i] = Index_Equ_IM[*count_Equ]; ModelBlock->Block_List[count_Block].Variable[i] = Index_Var_IM[*count_Equ]; + ModelBlock->Block_List[count_Block].Equation_Type[i] = Equation_Type[Index_Equ_IM[*count_Equ]].first; + ModelBlock->Block_List[count_Block].Equation_Type_Var[i] = Equation_Type[Index_Equ_IM[*count_Equ]].second; i_1 = Index_Var_IM[*count_Equ]; for (k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++) { @@ -571,6 +591,8 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const free(ModelBlock->Block_List[blk].Exogenous); free(ModelBlock->Block_List[blk].Own_Derivative); free(ModelBlock->Block_List[blk].Other_Endogenous); + free(ModelBlock->Block_List[blk].Equation_Type); + free(ModelBlock->Block_List[blk].Equation_Type_Var); for (i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) { if (incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size*/) @@ -605,34 +627,96 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const free(ModelBlock); } - -t_type -BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector equations ) +t_etype +BlockTriangular::Equation_Type_determination(vector &equations, map, NodeID> &first_cur_endo_derivatives, vector &Index_Var_IM, vector &Index_Equ_IM) { - int i=0; NodeID lhs, rhs; ostringstream tmp_output; BinaryOpNode *eq_node; ostringstream tmp_s; temporary_terms_type temporary_terms; + EquationType Equation_Simulation_Type; + t_etype V_Equation_Simulation_Type(equations.size()); + for(unsigned int i = 0; i < equations.size(); i++) + { + temporary_terms.clear(); + int eq = Index_Equ_IM[i]; + int var = Index_Var_IM[i]; + eq_node = equations[eq]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + Equation_Simulation_Type = E_SOLVE; + int Var_To_Derivate = -1; + tmp_s.str(""); + tmp_output.str(""); + lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); + tmp_s << "y(it_, " << Index_Var_IM[i]+1 << ")"; + map, NodeID>::iterator derivative = first_cur_endo_derivatives.find(make_pair(eq, var)); + set > result; + //result.clear(); + derivative->second->collectEndogenous(result); + set >::const_iterator d_endo_variable = result.find(make_pair(var, 0)); + //Determine whether the equation could be evaluated rather than to be solved + if (tmp_output.str() == tmp_s.str() and d_endo_variable == result.end()) + Equation_Simulation_Type = E_EVALUATE; + else + { //the equation could be normalized by a permutation of the rhs and the lhs + tmp_output.str(""); + rhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); + if (tmp_output.str() == tmp_s.str() and d_endo_variable == result.end()) + { + Equation_Simulation_Type = E_EVALUATE_R; + //cout << "Equation " << eq << " is reversed\n"; + } + else + { //the equation could be normalized using the derivative independant of the endogenous variable + if (d_endo_variable == result.end()) + { + Equation_Simulation_Type = E_EVALUATE_S; + Var_To_Derivate = var; + } + } + } + /*cout << "-----------------------------------------------------------\n"; + lhs->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms); + cout << " = "; + rhs->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms); + cout << "% " << c_Equation_Type(Equation_Simulation_Type) << " " << var+1 << "\n";*/ + V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, Var_To_Derivate); + } + return(V_Equation_Simulation_Type); +} + +t_type +BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector &equations, t_etype &Equation_Type) +{ + int i=0; int count_equ = 0, blck_count_simult =0; - int Blck_Size; + int Blck_Size, Recurs_Size; int Lead, Lag; t_type Type; bool *Cur_IM; BlockSimulationType Simulation_Type , prev_Type=UNKNOWN; + int eq = 0; for ( i=0; iget_arg1(); - rhs = eq_node->get_arg2(); - tmp_s.str(""); - tmp_output.str(""); - lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); - tmp_s << "y(it_, " << Index_Var_IM[first_count_equ]+1 << ")"; - //Determine whether the equation could be evaluated rather than to be solved - if (tmp_output.str()==tmp_s.str()) + if(Equation_Type[Index_Equ_IM[eq]].first==E_EVALUATE or Equation_Type[Index_Equ_IM[eq]].first==E_EVALUATE_R) { - if (Simulation_Type==SOLVE_BACKWARD_SIMPLE) - Simulation_Type=EVALUATE_BACKWARD; - else if (Simulation_Type==SOLVE_FORWARD_SIMPLE) - Simulation_Type=EVALUATE_FORWARD; - } - else - { - tmp_output.str(""); - rhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); - if (tmp_output.str()==tmp_s.str()) - { - if (Simulation_Type==SOLVE_BACKWARD_SIMPLE) - Simulation_Type=EVALUATE_BACKWARD_R; - else if (Simulation_Type==SOLVE_FORWARD_SIMPLE) - Simulation_Type=EVALUATE_FORWARD_R; - } + if (Simulation_Type == SOLVE_BACKWARD_SIMPLE) + Simulation_Type = EVALUATE_BACKWARD; + else if (Simulation_Type == SOLVE_FORWARD_SIMPLE) + Simulation_Type = EVALUATE_FORWARD; } if (i > 0) { - if ( ((prev_Type == EVALUATE_FORWARD_R || prev_Type == EVALUATE_FORWARD) && (Simulation_Type == EVALUATE_FORWARD_R || Simulation_Type == EVALUATE_FORWARD)) - || ((prev_Type == EVALUATE_BACKWARD_R || prev_Type == EVALUATE_BACKWARD) && (Simulation_Type == EVALUATE_BACKWARD_R || Simulation_Type == EVALUATE_BACKWARD)) - ) + if ((prev_Type == EVALUATE_FORWARD and Simulation_Type == EVALUATE_FORWARD) + or (prev_Type == EVALUATE_BACKWARD and Simulation_Type == EVALUATE_BACKWARD)) { BlockSimulationType c_Type = (Type[Type.size()-1]).first; - int c_Size = (Type[Type.size()-1]).second; - Type[Type.size()-1]=make_pair(c_Type, ++c_Size); + int c_Size = (Type[Type.size()-1]).second.first; + Type[Type.size()-1]=make_pair(c_Type, make_pair(++c_Size, Type[Type.size()-1].second.second)); } else - Type.push_back(make_pair(Simulation_Type, Blck_Size)); + Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size))); } else - Type.push_back(make_pair(Simulation_Type, Blck_Size)); + Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size))); } else { - Type.push_back(make_pair(Simulation_Type, Blck_Size)); + Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size))); } prev_Type = Simulation_Type; + eq += Blck_Size; } return(Type); } void -BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0, jacob_map j_m, vector equations ) +BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0, jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map, NodeID> &first_cur_endo_derivatives) { int i, j, Nb_TotalBlocks, Nb_RecursBlocks, Nb_SimulBlocks; BlockType Btype; int count_Block, count_Equ; - //block_result_t* res; - //Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set)); bool* SIM0, *SIM00; + SIM0 = (bool*)malloc(n * n * sizeof(bool)); memcpy(SIM0,IM_0,n*n*sizeof(bool)); Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0); @@ -763,9 +825,11 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, bool OK=false; double bi=0.99999999; int suppressed=0; + vector Index_Equ_IM_save(Index_Equ_IM); while (!OK && bi>1e-14) { int suppress=0; + Index_Equ_IM = Index_Equ_IM_save; SIM0 = (bool*)malloc(n * n * sizeof(bool)); memset(SIM0,0,n*n*sizeof(bool)); SIM00 = (bool*)malloc(n * n * sizeof(bool)); @@ -804,12 +868,14 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, } + V_Equation_Type = Equation_Type_determination(equations, first_cur_endo_derivatives, Index_Var_IM, Index_Equ_IM); + cout << "Finding the optimal block decomposition of the model ...\n"; vector > blocks; if (prologue+epiloguefirst==SOLVE_FORWARD_COMPLETE || it->first==SOLVE_BACKWARD_COMPLETE || it->first==SOLVE_TWO_BOUNDARIES_COMPLETE) { Nb_SimulBlocks++; - if (it->second>j) + if (it->second.first>j) { - j=it->second; + j=it->second.first; Nb_fv = blocks[Nb_SimulBlocks-1].second; } } @@ -845,13 +911,13 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, if (count_Equsecond==1) + if (it->second.first==1) Btype = PROLOGUE; else Btype = SIMULTANS; else Btype = EPILOGUE; - Allocate_Block(it->second, &count_Equ, count_Block++, Btype, it->first, ModelBlock); + Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second); } } @@ -860,7 +926,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, // normalize each equation of the dynamic model // and find the optimal block triangular decomposition of the static model void -BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m, vector equations) +BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector &equations, t_etype &equation_simulation_type, map, NodeID> &first_cur_endo_derivatives) { bool* SIM, *SIM_0; bool* Cur_IM; @@ -885,7 +951,6 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ cout << "incidence matrix for the static model (unsorted) \n"; incidencematrix.Print_SIM(SIM, eEndogenous); } - //Index_Equ_IM = (simple*)malloc(symbol_table.endo_nbr() * sizeof(*Index_Equ_IM)); Index_Equ_IM = vector(symbol_table.endo_nbr()); for (i = 0;i < symbol_table.endo_nbr();i++) { @@ -903,7 +968,7 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ SIM_0 = (bool*)malloc(symbol_table.endo_nbr() * symbol_table.endo_nbr() * sizeof(*SIM_0)); for (i = 0;i < symbol_table.endo_nbr()*symbol_table.endo_nbr();i++) SIM_0[i] = Cur_IM[i]; - Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr(), prologue, epilogue, Index_Var_IM, Index_Equ_IM, SIM_0, j_m, equations); + Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr(), prologue, epilogue, Index_Var_IM, Index_Equ_IM, SIM_0, j_m, equations, equation_simulation_type, first_cur_endo_derivatives); free(SIM_0); free(SIM); } diff --git a/preprocessor/BlockTriangular.hh b/preprocessor/BlockTriangular.hh index 9b3e26963..d350efafd 100644 --- a/preprocessor/BlockTriangular.hh +++ b/preprocessor/BlockTriangular.hh @@ -27,8 +27,7 @@ #include "ModelNormalization.hh" #include "ModelBlocks.hh" #include "IncidenceMatrix.hh" - - +#include "Modeltree.hh" @@ -37,7 +36,9 @@ //! Sparse matrix of double to store the values of the Jacobian typedef map,double> jacob_map; -typedef vector > t_type; +typedef vector > > t_type; + +typedef vector > t_etype; //! Creates the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition class BlockTriangular @@ -46,13 +47,15 @@ private: //! Find equations and endogenous variables belonging to the prologue and epilogue of the model void Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM0); //! Allocates and fills the Model structure describing the content of each block - void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock); + void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock, t_etype &Equation_Type, int recurs_Size); //! Finds a matching between equations and endogenous variables bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector &Index_Var_IM) const; //! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables - void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector &Index_Equ_IM, vector &Index_Var_IM, vector > &blocks, bool verbose_) const; + void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector &Index_Equ_IM, vector &Index_Var_IM, vector > &blocks, t_etype &Equation_Type, bool verbose_) const; + //! determine the type of each equation of the model (couble evaluated or need to be solved) + t_etype Equation_Type_determination(vector &equations, map, NodeID> &first_cur_endo_derivatives, vector &Index_Var_IM, vector &Index_Equ_IM); //! Tries to merge the consecutive blocks in a single block and determine the type of each block: recursive, simultaneous, ... - t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector equations ); + t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector &equations, t_etype &Equation_Type); public: const SymbolTable &symbol_table; BlockTriangular(const SymbolTable &symbol_table_arg); @@ -63,8 +66,8 @@ public: Blocks blocks; Normalization normalization; IncidenceMatrix incidencematrix; - void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m, vector equations); - void Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0 , jacob_map j_m, vector equations); + void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map, NodeID> &first_cur_endo_derivatives); + void Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0 , jacob_map &j_m, vector &equations, t_etype &equation_simulation_type, map, NodeID> &first_cur_endo_derivatives); vector Index_Equ_IM; vector Index_Var_IM; int prologue, epilogue; @@ -128,5 +131,17 @@ public: break; } }; + inline static std::string c_Equation_Type(int type) + { + char c_Equation_Type[5][13]= + { + "E_UNKNOWN ", + "E_EVALUATE ", + "E_EVALUATE_R", + "E_EVALUATE_S", + "E_SOLVE " + }; + return(c_Equation_Type[type]); + }; }; #endif diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh index aa043499a..1ecc837c7 100644 --- a/preprocessor/CodeInterpreter.hh +++ b/preprocessor/CodeInterpreter.hh @@ -49,6 +49,17 @@ enum BlockType SIMULTAN = 3 //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 << " 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 @@ -357,15 +357,30 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin { case EVALUATE_BACKWARD: case EVALUATE_FORWARD: + evaluation: output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel - << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; + << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl; output << " "; - output << tmp_output.str(); - output << " = "; - rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); + if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE) + { + output << tmp_output.str(); + output << " = "; + rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); + } + else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_R) + { + rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); + output << " = "; + output << tmp_output.str(); + output << "; %reversed " << ModelBlock->Block_List[j].Equation_Type[i] << " \n"; + } + else + { + cerr << "Type missmatch for equation " << ModelBlock->Block_List[j].Equation[i]+1 << "\n"; + exit(-1); + } output << ";\n"; - break; - case EVALUATE_BACKWARD_R: + 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; @@ -374,21 +389,25 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin output << " = "; lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); output << ";\n"; - break; + 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; output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel - << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; - output << " " << "residual(" << i+1 << ") = ("; + << " (" << 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; case SOLVE_TWO_BOUNDARIES_COMPLETE: case SOLVE_TWO_BOUNDARIES_SIMPLE: + if(iBlock_List[j].Nb_Recursives) + goto evaluation; 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_) = ("; + << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl; + Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_) = -residual(" << i+1 << ", it_)"; + output << " residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ", it_) = ("; goto end; default: end: @@ -2164,6 +2183,31 @@ DynamicModel::BlockLinear(Model_Block *ModelBlock) } +map, NodeID> +DynamicModel::collect_first_order_derivatives_current_endogenous() +{ + map, NodeID> curr_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); + /*cout << "eq=" << eq << " var=" << symbol_table.getName(var) << " (" << var << ") lag=" << lag; + ExprNodeOutputType output_type=oMatlabDynamicModelSparse; + const temporary_terms_type temporary_terms; + cout << " derivative = "; + (it2->second)->writeOutput(cout, output_type, temporary_terms); + cout << "\n";*/ + if(lag==0) + curr_endo_derivatives[make_pair(eq, var)] = it2->second; + } + } + return curr_endo_derivatives; +} + void @@ -2225,9 +2269,10 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative cout << "The gross incidence matrix \n"; block_triangular.incidencematrix.Print_IM(eEndogenous); } - block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations); - + t_etype equation_simulation_type; + map, NodeID> first_cur_endo_derivatives = collect_first_order_derivatives_current_endogenous(); + block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations, equation_simulation_type, first_cur_endo_derivatives); BlockLinear(block_triangular.ModelBlock); if (!no_tmp_terms) computeTemporaryTermsOrdered(block_triangular.ModelBlock); diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh index 107db04da..2576c9269 100644 --- a/preprocessor/DynamicModel.hh +++ b/preprocessor/DynamicModel.hh @@ -109,7 +109,8 @@ private: void computeParamsDerivatives(); //! Computes temporary terms for the file containing parameters derivatives void computeParamsDerivativesTemporaryTerms(); - + //! Collect only the first derivatives + map, NodeID> collect_first_order_derivatives_current_endogenous(); public: DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); diff --git a/preprocessor/DynareMain.cc b/preprocessor/DynareMain.cc index e00ee49a0..14f31d5a6 100644 --- a/preprocessor/DynareMain.cc +++ b/preprocessor/DynareMain.cc @@ -25,7 +25,9 @@ using namespace std; #include #include - +#ifndef PACKAGE_VERSION + #define PACKAGE_VERSION 4. +#endif #include "macro/MacroDriver.hh" /* Prototype for second part of main function diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh index b21642bc2..7b096d2cf 100644 --- a/preprocessor/ExprNode.hh +++ b/preprocessor/ExprNode.hh @@ -383,7 +383,7 @@ struct IM_compact //! One block of the model struct Block { - int Size, Sized, nb_exo, nb_exo_det, nb_other_endo; + int Size, Sized, nb_exo, nb_exo_det, nb_other_endo, Nb_Recursives; BlockType Type; BlockSimulationType Simulation_Type; int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo; @@ -392,7 +392,8 @@ struct Block int Max_Lag_Exo, Max_Lead_Exo; bool is_linear; int *Equation, *Own_Derivative; - int *Variable, *Other_Endogenous, *Exogenous; + EquationType *Equation_Type; + int *Variable, *Other_Endogenous, *Exogenous, *Equation_Type_Var; temporary_terms_type **Temporary_Terms_in_Equation; //temporary_terms_type *Temporary_terms; temporary_terms_inuse_type *Temporary_InUse;