From 8421712f80289663d7d88ed7cba7c26bf9723d05 Mon Sep 17 00:00:00 2001 From: ferhat Date: Fri, 10 Jul 2009 15:05:09 +0000 Subject: [PATCH] - Complete implementation of feedback variables in dynamic model with sparse option - Normalization of linear in endogenous variable equations git-svn-id: https://www.dynare.org/svn/dynare/trunk@2834 ac1d8469-bf42-47a9-8791-bf33cf982152 --- BlockTriangular.cc | 113 ++++++-- BlockTriangular.hh | 10 +- DynamicModel.cc | 603 ++++++++++++++++++++++-------------------- DynamicModel.hh | 2 + MinimumFeedbackSet.cc | 18 +- ModelTree.cc | 14 + ModelTree.hh | 5 + 7 files changed, 453 insertions(+), 312 deletions(-) diff --git a/BlockTriangular.cc b/BlockTriangular.cc index 6d8c6b76..557f0b42 100644 --- a/BlockTriangular.cc +++ b/BlockTriangular.cc @@ -167,8 +167,11 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog return check; } + + + t_vtype -BlockTriangular::Get_Variable_LeadLag_By_Block(vector &components_set, int nb_blck_sim, int prologue, int epilogue) const +BlockTriangular::Get_Variable_LeadLag_By_Block(vector &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const { int nb_endo = symbol_table.endo_nbr(); vector variable_blck(nb_endo), equation_blck(nb_endo); @@ -190,9 +193,9 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector &components_set, int variable_blck[Index_Var_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue); equation_blck[Index_Equ_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue); } - //cout << "equation_blck[" << Index_Equ_IM[i] << "]=" << equation_blck[Index_Equ_IM[i]] << " variable_blck[" << Index_Var_IM[i] << "] = " << variable_blck[Index_Var_IM[i]] << "\n"; Variable_Type[i] = make_pair(0, 0); } + equation_lead_lag = Variable_Type; for (int k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++) { bool *Cur_IM = incidencematrix.Get_IM(k, eEndogenous); @@ -203,19 +206,22 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector &components_set, int int i_1 = Index_Var_IM[i]; for (int j = 0; j < nb_endo; j++) { - if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[Index_Equ_IM[ j]]) + int j_l = Index_Equ_IM[ j]; + if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[j_l]) { if (k > Variable_Type[i_1].second) Variable_Type[i_1] = make_pair(Variable_Type[i_1].first, k); if (k < -Variable_Type[i_1].first) Variable_Type[i_1] = make_pair(-k, Variable_Type[i_1].second); + if (k > equation_lead_lag[j_l].second) + equation_lead_lag[j_l] = make_pair(equation_lead_lag[j_l].first, k); + if (k < -equation_lead_lag[j_l].first) + equation_lead_lag[j_l] = make_pair(-k, equation_lead_lag[j_l].second); } } } } } - /*for(int i = 0; i < nb_endo; i++) - cout << "Variable_Type[" << Index_Equ_IM[i] << "].first = " << Variable_Type[Index_Equ_IM[i]].first << " Variable_Type[" << Index_Equ_IM[i] << "].second=" << Variable_Type[Index_Equ_IM[i]].second << "\n";*/ return (Variable_Type); } @@ -255,7 +261,9 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo components_set[component[i]].first.insert(i); } - V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue); + + t_vtype equation_lead_lag; + V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue, equation_lead_lag); vector tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM); int order = prologue; @@ -265,7 +273,8 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo //Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set for (int i = 0; i < n; i++) - if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second >= 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0) + if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second > 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0 + or equation_lead_lag[Index_Equ_IM[i+prologue]].second > 0 or equation_lead_lag[Index_Equ_IM[i+prologue]].first > 0) add_edge(i, i, G2); //For each block, the minimum set of feedback variable is computed @@ -297,13 +306,14 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue]; order++; } + } free(AMp); free(SIM); } void -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) +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, vector &Index_Var_IM, vector &Index_Equ_IM) { 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; @@ -318,11 +328,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block 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].Chaine_Rule_Derivatives = new chaine_rule_derivatives_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_Normalized = (NodeID *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID)); + ModelBlock->Block_List[count_Block].Equation_Normalized = (NodeID*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID)); 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)); @@ -668,6 +679,7 @@ BlockTriangular::Free_Block(Model_Block *ModelBlock) const delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i]; free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation); delete (ModelBlock->Block_List[blk].Temporary_InUse); + delete ModelBlock->Block_List[blk].Chaine_Rule_Derivatives; } free(ModelBlock->Block_List); free(ModelBlock); @@ -735,7 +747,7 @@ BlockTriangular::Equation_Type_determination(vector &equations, BlockTriangular::Recompute_Deriavtives_Respect_to_Feedback_Variables( */ t_type -BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector &equations, t_etype &Equation_Type) +BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector &equations, t_etype &Equation_Type, vector &Index_Var_IM, vector &Index_Equ_IM) { int i = 0; int count_equ = 0, blck_count_simult = 0; @@ -842,6 +854,66 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue return (Type); } +map, pair, int> >, int> +BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck) +{ + map, pair, int> >, int> Derivatives; + Derivatives.clear(); + int nb_endo = symbol_table.endo_nbr(); + /*ModelBlock.Block_List[Blck].first_order_determinstic_simulation_derivatives = new*/ + for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++) + { + bool *IM=incidencematrix.Get_IM(lag, eEndogenous); + if(IM) + { + for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++) + { + int eqr = ModelBlock->Block_List[blck].Equation[eq]; + for(int var = 0; var < ModelBlock->Block_List[blck].Size; var++) + { + int varr = ModelBlock->Block_List[blck].Variable[var]; + /*cout << "IM=" << IM << "\n"; + cout << "varr=" << varr << " eqr=" << eqr << " lag=" << lag << "\n";*/ + if(IM[varr+eqr*nb_endo]) + { + + /*if(eqBlock_List[blck].Nb_Recursives and varBlock_List[blck].Nb_Recursives) + {*/ + bool OK = true; + map, pair, int> >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag))); + if(its!=Derivatives.end()) + { + if(its->second == 2) + OK=false; + } + + if(OK) + { + if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eqBlock_List[blck].Nb_Recursives) + Derivatives[make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag))] = 1; + else + Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varr, var), lag))] = 0; + } + /*} + else if(eqBlock_List[blck].Nb_Recursives and varBlock_List[blck].Nb_Recursives)*/ + if(varBlock_List[blck].Nb_Recursives) + { + int eqs = ModelBlock->Block_List[blck].Equation[var]; + for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; varsBlock_List[blck].Size; vars++) + { + int varrs = ModelBlock->Block_List[blck].Variable[vars]; + if(Derivatives.find(make_pair(make_pair(eqs, var), make_pair(make_pair(varrs, vars), lag)))!=Derivatives.end()) + Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varrs, vars), lag))] = 2; + } + } + } + } + } + } + } + return(Derivatives); +} + 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, t_etype &V_Equation_Type, map >, NodeID> &first_order_endo_derivatives) { @@ -850,6 +922,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int count_Block, count_Equ; 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); @@ -873,7 +946,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, //double bi=1e-13; int suppressed = 0; vector Index_Equ_IM_save(Index_Equ_IM); - while (!OK && bi > 1e-14) + while (!OK && bi > 1e-19) { int suppress = 0; Index_Equ_IM = Index_Equ_IM_save; @@ -907,7 +980,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, //bi/=1.07; bi /= 2; counted++; - if (bi > 1e-14) + if (bi > 1e-19) free(SIM00); } if (!OK) @@ -921,12 +994,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, if (prologue+epilogue < n) Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, V_Equation_Type, false); - t_type Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type); + t_type Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type, Index_Var_IM, Index_Equ_IM); i = 0; j = 0; Nb_SimulBlocks = 0; - int Nb_fv = 0; + int Nb_feedback_variable = 0; for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++) { if (it->first == SOLVE_FORWARD_COMPLETE || it->first == SOLVE_BACKWARD_COMPLETE || it->first == SOLVE_TWO_BOUNDARIES_COMPLETE) @@ -935,7 +1008,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, if (it->second.first > j) { j = it->second.first; - Nb_fv = blocks[Nb_SimulBlocks-1].second; + Nb_feedback_variable = blocks[Nb_SimulBlocks-1].second; } } } @@ -945,7 +1018,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, cout << Nb_TotalBlocks << " block(s) found:\n"; cout << " " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n"; cout << " the largest simultaneous block has " << j << " equation(s)\n" - <<" and " << Nb_fv << " feedback variable(s).\n"; + <<" and " << Nb_feedback_variable << " feedback variable(s).\n"; ModelBlock->Size = Nb_TotalBlocks; ModelBlock->Periods = periods; @@ -953,6 +1026,8 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, count_Equ = count_Block = 0; + + for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++) { if (count_Equ < prologue) @@ -961,10 +1036,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, if (it->second.first == 1) Btype = PROLOGUE; else - Btype = SIMULTANS; + { + Btype = SIMULTANS; + } else Btype = EPILOGUE; - Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second); + Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second, Index_Var_IM, Index_Equ_IM); } } diff --git a/BlockTriangular.hh b/BlockTriangular.hh index 51df4321..bd30c3c7 100644 --- a/BlockTriangular.hh +++ b/BlockTriangular.hh @@ -44,6 +44,8 @@ typedef vector > t_vtype; typedef set temporary_terms_inuse_type; +typedef vector, pair > > > chaine_rule_derivatives_type; + //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model struct IM_compact { @@ -71,6 +73,7 @@ struct Block //temporary_terms_type *Temporary_terms; temporary_terms_inuse_type *Temporary_InUse; IM_compact *IM_lead_lag; + chaine_rule_derivatives_type *Chaine_Rule_Derivatives; int Code_Start, Code_Length; }; @@ -92,7 +95,7 @@ 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, t_etype &Equation_Type, int recurs_Size); + 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, vector &Index_Var_IM, vector &Index_Equ_IM); //! 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 @@ -100,9 +103,9 @@ private: //! determines the type of each equation of the model (could be evaluated or need to be solved) t_etype Equation_Type_determination(vector &equations, map >, NodeID> &first_order_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_etype &Equation_Type); + t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector &equations, t_etype &Equation_Type, vector &Index_Var_IM, vector &Index_Equ_IM); //! Compute for each variable its maximum lead and lag in its block - t_vtype Get_Variable_LeadLag_By_Block(vector &components_set, int nb_blck_sim, int prologue, int epilogue) const; + t_vtype Get_Variable_LeadLag_By_Block(vector &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const; public: SymbolTable &symbol_table; /*Blocks blocks; @@ -115,6 +118,7 @@ public: //! Frees the Model structure describing the content of each block void Free_Block(Model_Block* ModelBlock) const; + map, pair, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck); void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map >, NodeID> &first_order_endo_derivatives); diff --git a/DynamicModel.cc b/DynamicModel.cc index 370a4d54..75b7ce9e 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -105,6 +105,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock) ostringstream tmp_output; BinaryOpNode *eq_node; first_derivatives_type::const_iterator it; + first_chaine_rule_derivatives_type::const_iterator it_chr; ostringstream tmp_s; temporary_terms.clear(); @@ -132,6 +133,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock) //if(it!=first_derivatives.end()) it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx); } + } + for(i=0; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + { + pair< pair, pair > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i); + eq=it.first.second; + var=it.second.second.first; + lag=it.second.second.second; + it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag))); + //it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + it_chr->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx); } /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) { @@ -187,6 +198,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock) //if(it!=first_derivatives.end()) it->second->collectTemporary_terms(temporary_terms, ModelBlock, j); } + } + for(i=0; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + { + pair< pair, pair > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i); + eq=it.first.second; + var=it.second.second.first; + lag=it.second.second.second; + it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag))); + //it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + it_chr->second->collectTemporary_terms(temporary_terms, ModelBlock, j); } /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) { @@ -240,13 +261,13 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin ofstream output; //temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); int nze, nze_exo, nze_other_endo; - map recursive_variables; + //map recursive_variables; vector feedback_variables; //---------------------------------------------------------------------- //For each block for (j = 0;j < ModelBlock->Size;j++) { - recursive_variables.clear(); + //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; @@ -388,14 +409,17 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin lhs = eq_node->get_arg1(); rhs = eq_node->get_arg2(); tmp_output.str(""); - lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); + 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, oMatlabDynamicModelSparse, temporary_terms); + else + lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); switch (ModelBlock->Block_List[j].Simulation_Type) { 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 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl; +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) { @@ -418,9 +442,19 @@ evaluation: { rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); output << "\n "; - temporary_terms_type tt2; - tt2.clear(); - ModelBlock->Block_List[j].Equation_Normalized[i]->writeOutput(output , oMatlabDynamicModelSparse, temporary_terms/*tt2*/); + /*temporary_terms_type tt2; + tt2.clear();*/ + tmp_output.str(""); + eq_node = (BinaryOpNode *)ModelBlock->Block_List[j].Equation_Normalized[i]; + lhs = eq_node->get_arg1(); + rhs = eq_node->get_arg2(); + if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD) + lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); + else + lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); + output << tmp_output.str(); + output << " = "; + rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); } } else @@ -436,10 +470,10 @@ evaluation: case SOLVE_FORWARD_COMPLETE: if (iBlock_List[j].Nb_Recursives) { - if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S) + /*if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S) recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i]; else - recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]]; + recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];*/ goto evaluation; } feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]); @@ -451,10 +485,10 @@ evaluation: case SOLVE_TWO_BOUNDARIES_SIMPLE: if (iBlock_List[j].Nb_Recursives) { - if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S) + /*if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S) recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i]; else - recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]]; + recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];*/ goto evaluation; } feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]); @@ -618,79 +652,44 @@ end: m=ModelBlock->Block_List[j].Max_Lag; //cout << "\nDerivatives in Block " << j << "\n"; - for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + for(i=0; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + { + //Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag)))); + pair< pair, pair > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i); + int eqr=it.first.first; + int eq=it.first.second; + int varr=it.second.first; + int var=it.second.second.first; + k=it.second.second.second; + /*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]; int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - bool derivative_exist; ostringstream tmp_output; if (eqrBlock_List[j].Nb_Recursives) { if (varr>=ModelBlock->Block_List[j].Nb_Recursives) - { - /*tmp_output << " g1_tmp_r(" << eqr+1 << ", " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = "; - NodeID tmp_n; - if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE) - tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]]; - else - tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr]; - int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0); - NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables); - ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); - output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) - << "(" << k - << ") " << var+1 - << ", equation=" << eq+1 << endl;*/ - } - } - else - { - if (varr>=ModelBlock->Block_List[j].Nb_Recursives) - { + {*/ output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", " << varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = "; - /*writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);*/ - /*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE) - derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables); - else - derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables); - //if (derivative_exist) - output << tmp_output.str() << ";";*/ - NodeID tmp_n; - //if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE) - tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]]; - /*else - tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/ - //cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"; - //cout << "derivaive eq=" << eq << " var=" << var << " k0=" << k << "\n"; - int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0); - NodeID ChaineRule_Derivative = tmp_n->getChainRuleDerivative(deriv_id, recursive_variables); - ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); + writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); output << ";"; output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; } + /*} } - /*if (eqrBlock_List[j].Nb_Recursives or varrBlock_List[j].Nb_Recursives) - output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = "; - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getName(var) - << "(0) " << var+1 - << ", equation=" << eq+1 << endl;*/ - } + }*/ output << " end;\n"; - //output << " ya = y;\n"; break; case SOLVE_TWO_BOUNDARIES_SIMPLE: case SOLVE_TWO_BOUNDARIES_COMPLETE: output << " if ~jacobian_eval" << endl; - for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) { k=m-ModelBlock->Block_List[j].Max_Lag; for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++) @@ -698,163 +697,86 @@ end: int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - bool derivative_exist; - ostringstream tmp_output; - //cout << "ModelBlock->Block_List[" << j << "].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n"; - if (eqrBlock_List[j].Nb_Recursives) - { - /*if (varrBlock_List[j].Nb_Recursives) - { - if (k==0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1 - << ", " << varr+1 - << "+" << ModelBlock->Block_List[j].Size*ModelBlock->Block_List[j].Max_Lag << ")*y(it_, " << var+1 << ")"; - else if (k>0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1 - << ", " << varr+1 - << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_+" << k << ", " << var+1 << ")"; - else if (k<0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1 - << ", " << varr+1 - << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_" << k << ", " << var+1 << ")"; - if (k==0) - tmp_output << " g1_tmp_r(" << eqr+1 << ", " - << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = "; - else if (k>0) - tmp_output << " g1_tmp_r(" << eqr+1 << ", " - << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = "; - else if (k<0) - tmp_output << " g1_tmp_r(" << eqr+1 << ", " - << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = "; - if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE) - derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables); - else - { - BinaryOpNode* tt = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[eqr]; - derivative_exist = tt->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables); - } - if (derivative_exist) - output << tmp_output.str() << ";"; - else - { - //output << "1" << ";"; - if (ModelBlock->Block_List[j].Equation_Type[eqr] != E_EVALUATE) - { + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];*/ + for(i=0; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + { + //Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag)))); + pair< pair, pair > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i); + int eqr=it.first.first; + int eq=it.first.second; + int varr=it.second.first; + int var=it.second.second.first; + k=it.second.second.second; - } - } - //writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); - output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << " derivative_exist=" << derivative_exist << " varr+1=" << varr+1 << endl; - }*/ - } - else + + + //bool derivative_exist; + ostringstream tmp_output; + /*if (eqr>=ModelBlock->Block_List[j].Nb_Recursives) { if (varr>=ModelBlock->Block_List[j].Nb_Recursives) + {*/ + /*for(int equation = ModelBlock->Block_List[j].Nb_Recursives; equationBlock_List[j].Size; equation++) + { + int eq = ModelBlock->Block_List[j].Equation[equation]; + int eqr = equation - ModelBlock->Block_List[j].Nb_Recursives; + for(int variable = ModelBlock->Block_List[j].Nb_Recursives; variableBlock_List[j].Size; variable++) { - if (k==0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+Per_K_)*y(it_, " << var+1 << ")"; - else if (k==1) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+Per_y_)*y(it_+1, " << var+1 << ")"; - else if (k>0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")"; - else if (k<0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives - << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; - if (k==0) - tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = "; - else if (k==1) - tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = "; - else if (k>0) - tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = "; - else if (k<0) - tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " - << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = "; - /*NodeID tmp_n; - if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE) - tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]]; - else - tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/ - /*int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),k); - //cout << "-------------------------------------------------------------------------------------\n"; - //cout << "derivaive eq=" << eq << " var=" << var << " k=" << k << "\n"; - //output << " " << tmp_output.str(); - map recursive_variables_save(recursive_variables); - NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables, var, k); - recursive_variables = recursive_variables_save; - //ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); - */ - /*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE) - derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables); - else - derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables); - if (derivative_exist) - output << tmp_output.str() << ";";*/ + int var = ModelBlock->Block_List[j].Variable[variable]; + int varr = variable - ModelBlock->Block_List[j].Nb_Recursives;*/ + //cout << "eqr=" << eqr << " varr=" << varr; + //cout << "k=" << k << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << " ModelBlock->Block_List[j].Equation[eq]=" << ModelBlock->Block_List[j].Equation[eq] << "\n"; + if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives) + { + + if (k==0) + Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives + << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives + << "+Per_K_)*y(it_, " << varr+1 << ")"; + else if (k==1) + Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives + << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives + << "+Per_y_)*y(it_+1, " << varr+1 << ")"; + else if (k>0) + Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives + << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives + << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << varr+1 << ")"; + else if (k<0) + Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives + << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives + << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << varr+1 << ")"; + if (k==0) + tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " + << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = "; + else if (k==1) + tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " + << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = "; + else if (k>0) + tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " + << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = "; + else if (k<0) + tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, " + << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = "; + + + output << " " << tmp_output.str(); + + writeChaineRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms); - /*output << tmp_output.str();*/ - //output << ";"; - //output << "\n%"; - output << tmp_output.str(); - writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms); output << ";"; - - output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - /*else - { - if (k==0) - tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", " - << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = "; - else if (k>0) - tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", " - << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = "; - else if (k<0) - tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", " - << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = "; - if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE) - derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables); - else - derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables); - if (derivative_exist) - output << tmp_output.str() << ";"; - - if (k==0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives - << ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag - << ")*y(it_, " << var+1 << ")"; - else if (k>0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives - << ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag - << ")*y(it_+" << k << ", " << var+1 << ")"; - else if (k<0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives - << ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag - << ")*y(it_" << k << ", " << var+1 << ")"; - output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; - }*/ - } + output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr)) + << "(" << k << ") " << varr+1 + << ", equation=" << eqr+1 << endl; + } + //cout << " done\n"; + /* } + }*/ #ifdef CONDITION output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++) { @@ -975,12 +897,10 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model bool lhs_rhs_done; Uff Uf[symbol_table.endo_nbr()]; map reference_count; - map ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number; + //map ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number; + vector feedback_variables; int prev_Simulation_Type=-1; - //SymbolicGaussElimination SGE; bool file_open=false; - //temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); - //---------------------------------------------------------------------- string main_name=file_name; main_name+=".cod"; code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate ); @@ -993,52 +913,26 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model code_file.write(&FDIMT, sizeof(FDIMT)); k=temporary_terms.size(); code_file.write(reinterpret_cast(&k),sizeof(k)); - //search for successive and identical blocks - i=k=k0=0; - ModelBlock_Aggregated_Count=-1; - for (j = 0;j < ModelBlock->Size;j++) - { - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) - && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD - ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD - /*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R - ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R */)) - { - } - else - { - k=k0=0; - ModelBlock_Aggregated_Count++; - } - k0+=ModelBlock->Block_List[j].Size; - ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0; - ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k; - prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; - } - ModelBlock_Aggregated_Count++; + +ModelBlock_Aggregated_Count = ModelBlock->Size; //For each block - j=0; - for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++) + + for (j = 0; j < ModelBlock->Size ;j++) { - k1=j; - if (k0>0) + feedback_variables.clear(); + if (j>0) code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK)); - v=ModelBlock_Aggregated_Number[k0]; + v=ModelBlock->Block_List[j].Size; code_file.write(reinterpret_cast(&v),sizeof(v)); v=ModelBlock->Block_List[j].Simulation_Type; code_file.write(reinterpret_cast(&v),sizeof(v)); - for (k=0; kBlock_List[j].Size;i++) { - for (i=0; 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])); - } - j++; + 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])); } - j=k1; if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE) { @@ -1051,8 +945,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model 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)); - //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) - //{ int u_count_int=0; Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open, ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE); @@ -1061,8 +953,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model file_open=true; //} } - for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++) - { //For a block composed of a single equation determines whether we have to evaluate or to solve the equation if (ModelBlock->Block_List[j].Size==1) { @@ -1076,7 +966,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model // The equations for (i = 0;i < ModelBlock->Block_List[j].Size;i++) { - //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); //The Temporary terms temporary_terms_type tt2; #ifdef DEBUGC @@ -1116,6 +1005,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model } switch (ModelBlock->Block_List[j].Simulation_Type) { +evaluation: case EVALUATE_BACKWARD: case EVALUATE_FORWARD: if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE) @@ -1133,19 +1023,13 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model lhs->compile(code_file, true, temporary_terms, map_idx); } break; - /*case EVALUATE_BACKWARD_R: - case EVALUATE_FORWARD_R: - lhs->compile(code_file, false, temporary_terms, map_idx); - rhs->compile(code_file, true, temporary_terms, map_idx); - break;*/ case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: - v=ModelBlock->Block_List[j].Equation[i]; - Uf[v].eqr=i; - Uf[v].Ufl=NULL; - goto end; case SOLVE_TWO_BOUNDARIES_COMPLETE: case SOLVE_TWO_BOUNDARIES_SIMPLE: + 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].eqr=i; Uf[v].Ufl=NULL; @@ -1167,6 +1051,7 @@ end: } code_file.write(&FENDEQU, sizeof(FENDEQU)); // The Jacobian if we have to solve the block + bool feedback_variable = (feedback_variables.size()>0); if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD /*&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R @@ -1183,30 +1068,73 @@ end: break; case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: - m=ModelBlock->Block_List[j].Max_Lag; - for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + if(feedback_variable) { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - int v=ModelBlock->Block_List[j].Equation[eqr]; - if (!Uf[v].Ufl) + int u = feedback_variables.size(); + for(i=0; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) { - Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); - Uf[v].Ufl_First=Uf[v].Ufl; + //Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag)))); + pair< pair, pair > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i); + int eqr=it.first.first; + int eq=it.first.second; + int varr=it.second.first; + int var=it.second.second.first; + int v=ModelBlock->Block_List[j].Equation[eqr]; + k=it.second.second.second; + if (!Uf[v].Ufl) + { + Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl_First=Uf[v].Ufl; + } + else + { + Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl=Uf[v].Ufl->pNext; + } + Uf[v].Ufl->pNext=NULL; + Uf[v].Ufl->u=u; + Uf[v].Ufl->var=var; + compileDerivative(code_file, eq, var, 0, map_idx); + code_file.write(&FSTPU, sizeof(FSTPU)); + code_file.write(reinterpret_cast(&u), sizeof(u)); + u++; + /*output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", " + << varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = "; + writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + output << ";"; + output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) + << "(" << k + << ") " << var+1 + << ", equation=" << eq+1 << endl;*/ } - else + } + else + { + m=ModelBlock->Block_List[j].Max_Lag; + for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { - Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); - Uf[v].Ufl=Uf[v].Ufl->pNext; + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; + int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + int v=ModelBlock->Block_List[j].Equation[eqr]; + if (!Uf[v].Ufl) + { + Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl_First=Uf[v].Ufl; + } + else + { + Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l)); + Uf[v].Ufl=Uf[v].Ufl->pNext; + } + Uf[v].Ufl->pNext=NULL; + Uf[v].Ufl->u=u; + Uf[v].Ufl->var=var; + compileDerivative(code_file, eq, var, 0, map_idx); + code_file.write(&FSTPU, sizeof(FSTPU)); + code_file.write(reinterpret_cast(&u), sizeof(u)); } - Uf[v].Ufl->pNext=NULL; - Uf[v].Ufl->u=u; - Uf[v].Ufl->var=var; - compileDerivative(code_file, eq, var, 0, map_idx); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast(&u), sizeof(u)); } for (i = 0;i < ModelBlock->Block_List[j].Size;i++) { @@ -1340,8 +1268,6 @@ end: prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; } - j++; - } } code_file.write(&FENDBLOCK, sizeof(FENDBLOCK)); code_file.write(&FEND, sizeof(FEND)); @@ -1993,7 +1919,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const hessian_output << "v2"; hessianHelper(hessian_output, k, 2, output_type); - hessian_output << "=v2"; + hessian_output << "=v2"; hessianHelper(hessian_output, k-1, 2, output_type); hessian_output << ";" << endl; @@ -2041,7 +1967,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const k += k2; } - + if (mode == eStandardMode) { DynamicOutput << "%" << endl @@ -2597,6 +2523,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative BlockLinear(block_triangular.ModelBlock); if (!no_tmp_terms) computeTemporaryTermsOrdered(block_triangular.ModelBlock); + + computeChaineRuleJacobian(block_triangular.ModelBlock); } else if (!no_tmp_terms) @@ -2817,6 +2745,107 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti return it->second; } + +void +DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock) +{ + //cout << "computeChaineRuleJacobian\n"; + //clock_t t1 = clock(); + map recursive_variables; + first_chaine_rule_derivatives.clear(); + for(int blck = 0; blckSize; blck++) + { + //cout << "blck=" << blck << "\n"; + 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) + { + //cout << "SOLVE_TWO_BOUNDARIES_COMPLETE \n"; + ModelBlock->Block_List[blck].Chaine_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]]; + } + //cout << "After recursive_alloc\n"; + map, pair,int> > , int > Derivatives = block_triangular.get_Derivatives(ModelBlock, blck); + for(map, pair,int> > , int >::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++) + { + int eqr = it->first.first.first; + int eq = it->first.first.second; + int varr = it->first.second.first.first; + int var = it->first.second.first.second; + int lag = it->first.second.second; + int Deriv_type = it->second; + if(Deriv_type == 0) + first_chaine_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_chaine_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 and eqBlock_List[blck].Nb_Recursives) + first_chaine_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_chaine_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].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag)))); + } + + + + /*for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++) + { + + for(int eq = 0; 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]->getChaineRuleDerivative(recursive_variables, varr, lag); + if (d1 == Zero) + continue; + first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = d1; + ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag)))); + } + } + }*/ + + + } + 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) + { + //cout << "SOLVE_FORWARD_SIMPLE \n"; + ModelBlock->Block_List[blck].Chaine_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_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1; + ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, 0)))); + } + } + } + } + //cout << "elapsed time in milliseconds = " << 1000.0*(double(clock()) - double(t1))/double(CLOCKS_PER_SEC) << "\n"; +} + + + void DynamicModel::computeParamsDerivatives() { diff --git a/DynamicModel.hh b/DynamicModel.hh index c0fb0602..7ecefae5 100644 --- a/DynamicModel.hh +++ b/DynamicModel.hh @@ -116,6 +116,8 @@ private: int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException); //! Compute the column indices of the dynamic Jacobian void computeDynJacobianCols(bool jacobianExo); + //! Computes chaine rule derivatives of the Jacobian w.r. to endogenous variables + void computeChaineRuleJacobian(Model_Block *ModelBlock); //! Computes derivatives of the Jacobian w.r. to parameters void computeParamsDerivatives(); //! Computes temporary terms for the file containing parameters derivatives diff --git a/MinimumFeedbackSet.cc b/MinimumFeedbackSet.cc index 33dcd3d2..d39f89da 100644 --- a/MinimumFeedbackSet.cc +++ b/MinimumFeedbackSet.cc @@ -26,6 +26,8 @@ namespace MFS void Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G) { + /*clear all in and out edges of vertex_to_eliminate + and remove vertex_to_eliminate from the graph*/ clear_vertex(vertex_to_eliminate, G); remove_vertex(vertex_to_eliminate, G); } @@ -41,6 +43,7 @@ namespace MFS void Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G) { + /*before the vertex i suppression replace all edges e_k_i and e_i_j by e_k_j*/ if (in_degree (vertex_to_eliminate, G) > 0 && out_degree (vertex_to_eliminate, G) > 0) { AdjacencyList_type::in_edge_iterator it_in, in_end; @@ -66,11 +69,14 @@ namespace MFS color[u] = gray_color; graph_traits::out_edge_iterator vi, vi_end; for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi) - if (color[target(*vi, g)] == white_color && has_cycle_dfs(g, target(*vi, g), color, circuit_stack)) + if (color[target(*vi, g)] == white_color) { - // cycle detected, return immediately - circuit_stack.push_back(v_index[target(*vi, g)]); - return true; + if (has_cycle_dfs(g, target(*vi, g), color, circuit_stack)) + { + // cycle detected, return immediately + circuit_stack.push_back(v_index[target(*vi, g)]); + return true; + } } else if (color[target(*vi, g)] == gray_color) { @@ -209,6 +215,8 @@ namespace MFS vector_vertex_descriptor Collect_Doublet(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G) { + /*collect all doublet (for each edge e_i_k there is an edge e_k_i with k!=i) in the graph + and return the vector of doublet*/ AdjacencyList_type::in_edge_iterator it_in, in_end; AdjacencyList_type::out_edge_iterator it_out, out_end; vector Doublet; @@ -223,6 +231,7 @@ namespace MFS bool Vertex_Belong_to_a_Clique(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G) { + /*Detect all the clique (all vertex in a clique are related to each other) in the graph*/ vector liste; bool agree = true; AdjacencyList_type::in_edge_iterator it_in, in_end; @@ -262,6 +271,7 @@ namespace MFS bool Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G) { + /*Graph reduction: eliminating purely intermediate variables or variables outside of any circuit*/ bool something_has_been_done = false; bool not_a_loop; int i; diff --git a/ModelTree.cc b/ModelTree.cc index 7ad95f20..c590e68e 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -65,6 +65,20 @@ ModelTree::computeJacobian(const set &vars) } } +void +ModelTree::writeChaineRuleDerivative(ostream &output, int eq, int var, int lag, + ExprNodeOutputType output_type, + const temporary_terms_type &temporary_terms) const +{ + map >, NodeID>::const_iterator it = first_chaine_rule_derivatives.find(make_pair(eq, make_pair(var, lag))); + if (it != first_chaine_rule_derivatives.end()) + (it->second)->writeOutput(output, output_type, temporary_terms); + else + output << 0; +} + + + void ModelTree::computeHessian(const set &vars) { diff --git a/ModelTree.hh b/ModelTree.hh index c4c759fd..7037778c 100644 --- a/ModelTree.hh +++ b/ModelTree.hh @@ -47,6 +47,9 @@ protected: */ first_derivatives_type first_derivatives; + typedef map< pair< int, pair< int, int> >, NodeID> first_chaine_rule_derivatives_type; + first_chaine_rule_derivatives_type first_chaine_rule_derivatives; + typedef map >, NodeID> second_derivatives_type; //! Second order derivatives /*! First index is equation number, second and third are variables w.r. to which is computed the derivative. @@ -80,6 +83,8 @@ protected: //! Write derivative of an equation w.r. to a variable void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; + //! Write chaine rule derivative of a recursive equation w.r. to a variable + void writeChaineRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; //! Computes temporary terms (for all equations and derivatives) void computeTemporaryTerms(bool is_matlab); //! Writes temporary terms