Commit d19ccced authored by ferhat's avatar ferhat
Browse files

- 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
parent 4ba6a505
...@@ -167,8 +167,11 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog ...@@ -167,8 +167,11 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
return check; return check;
} }
t_vtype t_vtype
BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const
{ {
int nb_endo = symbol_table.endo_nbr(); int nb_endo = symbol_table.endo_nbr();
vector<int> variable_blck(nb_endo), equation_blck(nb_endo); vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
...@@ -190,9 +193,9 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int ...@@ -190,9 +193,9 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
variable_blck[Index_Var_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue); 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); 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); 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++) for (int k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{ {
bool *Cur_IM = incidencematrix.Get_IM(k, eEndogenous); bool *Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
...@@ -203,19 +206,22 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int ...@@ -203,19 +206,22 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
int i_1 = Index_Var_IM[i]; int i_1 = Index_Var_IM[i];
for (int j = 0; j < nb_endo; j++) 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) if (k > Variable_Type[i_1].second)
Variable_Type[i_1] = make_pair(Variable_Type[i_1].first, k); Variable_Type[i_1] = make_pair(Variable_Type[i_1].first, k);
if (k < -Variable_Type[i_1].first) if (k < -Variable_Type[i_1].first)
Variable_Type[i_1] = make_pair(-k, Variable_Type[i_1].second); 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); return (Variable_Type);
} }
...@@ -255,7 +261,9 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo ...@@ -255,7 +261,9 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
components_set[component[i]].first.insert(i); 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<int> tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM); vector<int> tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
int order = prologue; int order = prologue;
...@@ -265,7 +273,8 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo ...@@ -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 //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++) 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); add_edge(i, i, G2);
//For each block, the minimum set of feedback variable is computed //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 ...@@ -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]; Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue];
order++; order++;
} }
} }
free(AMp); free(AMp);
free(SIM); free(SIM);
} }
void 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<int> &Index_Var_IM, vector<int> &Index_Equ_IM)
{ {
int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li; 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; 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 ...@@ -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].Type = type;
ModelBlock->Block_List[count_Block].Nb_Recursives = recurs_Size; 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 = 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].Temporary_InUse->clear();
ModelBlock->Block_List[count_Block].Simulation_Type = SimType; 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 = (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 = (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].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].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)); 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 ...@@ -668,6 +679,7 @@ BlockTriangular::Free_Block(Model_Block *ModelBlock) const
delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i]; delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation); free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
delete (ModelBlock->Block_List[blk].Temporary_InUse); delete (ModelBlock->Block_List[blk].Temporary_InUse);
delete ModelBlock->Block_List[blk].Chaine_Rule_Derivatives;
} }
free(ModelBlock->Block_List); free(ModelBlock->Block_List);
free(ModelBlock); free(ModelBlock);
...@@ -735,7 +747,7 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations, ...@@ -735,7 +747,7 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations,
BlockTriangular::Recompute_Deriavtives_Respect_to_Feedback_Variables( BlockTriangular::Recompute_Deriavtives_Respect_to_Feedback_Variables(
*/ */
t_type t_type
BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type) BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM)
{ {
int i = 0; int i = 0;
int count_equ = 0, blck_count_simult = 0; int count_equ = 0, blck_count_simult = 0;
...@@ -842,6 +854,66 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue ...@@ -842,6 +854,66 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
return (Type); return (Type);
} }
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>
BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
{
map<pair<pair<int, int>, pair<pair<int, int>, 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(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)
{*/
bool OK = true;
map<pair<pair<int, int>, pair<pair<int, int>, 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 eq<ModelBlock->Block_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(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)*/
if(var<ModelBlock->Block_List[blck].Nb_Recursives)
{
int eqs = ModelBlock->Block_List[blck].Equation[var];
for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; vars<ModelBlock->Block_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 void
BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool *IM_0, jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives) BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool *IM_0, jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives)
{ {
...@@ -850,6 +922,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, ...@@ -850,6 +922,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
int count_Block, count_Equ; int count_Block, count_Equ;
bool *SIM0, *SIM00; bool *SIM0, *SIM00;
SIM0 = (bool *) malloc(n * n * sizeof(bool)); SIM0 = (bool *) malloc(n * n * sizeof(bool));
memcpy(SIM0, IM_0, 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); 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, ...@@ -873,7 +946,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
//double bi=1e-13; //double bi=1e-13;
int suppressed = 0; int suppressed = 0;
vector<int> Index_Equ_IM_save(Index_Equ_IM); vector<int> Index_Equ_IM_save(Index_Equ_IM);
while (!OK && bi > 1e-14) while (!OK && bi > 1e-19)
{ {
int suppress = 0; int suppress = 0;
Index_Equ_IM = Index_Equ_IM_save; Index_Equ_IM = Index_Equ_IM_save;
...@@ -907,7 +980,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, ...@@ -907,7 +980,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
//bi/=1.07; //bi/=1.07;
bi /= 2; bi /= 2;
counted++; counted++;
if (bi > 1e-14) if (bi > 1e-19)
free(SIM00); free(SIM00);
} }
if (!OK) if (!OK)
...@@ -921,12 +994,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, ...@@ -921,12 +994,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
if (prologue+epilogue < n) 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); 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; i = 0;
j = 0; j = 0;
Nb_SimulBlocks = 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++) 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) 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, ...@@ -935,7 +1008,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
if (it->second.first > j) if (it->second.first > j)
{ {
j = it->second.first; 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, ...@@ -945,7 +1018,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
cout << Nb_TotalBlocks << " block(s) found:\n"; cout << Nb_TotalBlocks << " block(s) found:\n";
cout << " " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n"; cout << " " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n";
cout << " the largest simultaneous block has " << j << " equation(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->Size = Nb_TotalBlocks;
ModelBlock->Periods = periods; ModelBlock->Periods = periods;
...@@ -953,6 +1026,8 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, ...@@ -953,6 +1026,8 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
count_Equ = count_Block = 0; count_Equ = count_Block = 0;
for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++) for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
{ {
if (count_Equ < prologue) if (count_Equ < prologue)
...@@ -961,10 +1036,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, ...@@ -961,10 +1036,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
if (it->second.first == 1) if (it->second.first == 1)
Btype = PROLOGUE; Btype = PROLOGUE;
else else
Btype = SIMULTANS; {
Btype = SIMULTANS;
}
else else
Btype = EPILOGUE; 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);
} }
} }
......
...@@ -44,6 +44,8 @@ typedef vector<pair< int, int> > t_vtype; ...@@ -44,6 +44,8 @@ typedef vector<pair< int, int> > t_vtype;
typedef set<int> temporary_terms_inuse_type; typedef set<int> temporary_terms_inuse_type;
typedef vector<pair< pair<int, int>, pair<int, pair<int, int> > > > chaine_rule_derivatives_type;
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
struct IM_compact struct IM_compact
{ {
...@@ -71,6 +73,7 @@ struct Block ...@@ -71,6 +73,7 @@ struct Block
//temporary_terms_type *Temporary_terms; //temporary_terms_type *Temporary_terms;
temporary_terms_inuse_type *Temporary_InUse; temporary_terms_inuse_type *Temporary_InUse;
IM_compact *IM_lead_lag; IM_compact *IM_lead_lag;
chaine_rule_derivatives_type *Chaine_Rule_Derivatives;
int Code_Start, Code_Length; int Code_Start, Code_Length;
}; };
...@@ -92,7 +95,7 @@ private: ...@@ -92,7 +95,7 @@ private:
//! Find equations and endogenous variables belonging to the prologue and epilogue of the model //! 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<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM0); void Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM0);
//! Allocates and fills the Model structure describing the content of each block //! 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<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Finds a matching between equations and endogenous variables //! 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<int> &Index_Var_IM) const; bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector<int> &Index_Var_IM) const;
//! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables //! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables
...@@ -100,9 +103,9 @@ private: ...@@ -100,9 +103,9 @@ private:
//! determines the type of each equation of the model (could be evaluated or need to be solved) //! determines the type of each equation of the model (could be evaluated or need to be solved)
t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM); t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Tries to merge the consecutive blocks in a single block and determine the type of each block: recursive, simultaneous, ... //! 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<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type); t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Compute for each variable its maximum lead and lag in its block //! Compute for each variable its maximum lead and lag in its block
t_vtype Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const; t_vtype Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const;
public: public:
SymbolTable &symbol_table; SymbolTable &symbol_table;
/*Blocks blocks; /*Blocks blocks;
...@@ -115,6 +118,7 @@ public: ...@@ -115,6 +118,7 @@ public:
//! Frees the Model structure describing the content of each block //! Frees the Model structure describing the content of each block
void Free_Block(Model_Block* ModelBlock) const; void Free_Block(Model_Block* ModelBlock) const;
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives); void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives);
......
This diff is collapsed.
...@@ -116,6 +116,8 @@ private: ...@@ -116,6 +116,8 @@ private:
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException); int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the dynamic Jacobian //! Compute the column indices of the dynamic Jacobian
void computeDynJacobianCols(bool jacobianExo); 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 //! Computes derivatives of the Jacobian w.r. to parameters
void computeParamsDerivatives(); void computeParamsDerivatives();
//! Computes temporary terms for the file containing parameters derivatives //! Computes temporary terms for the file containing parameters derivatives
......
...@@ -26,6 +26,8 @@ namespace MFS ...@@ -26,6 +26,8 @@ namespace MFS
void void
Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G) 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); clear_vertex(vertex_to_eliminate, G);
remove_vertex(vertex_to_eliminate, G); remove_vertex(vertex_to_eliminate, G);
} }
...@@ -41,6 +43,7 @@ namespace MFS ...@@ -41,6 +43,7 @@ namespace MFS
void void
Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G) 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) if (in_degree (vertex_to_eliminate, G) > 0 && out_degree (vertex_to_eliminate, G) > 0)
{ {
AdjacencyList_type::in_edge_iterator it_in, in_end; AdjacencyList_type::in_edge_iterator it_in, in_end;
...@@ -66,11 +69,14 @@ namespace MFS ...@@ -66,11 +69,14 @@ namespace MFS
color[u] = gray_color; color[u] = gray_color;
graph_traits<AdjacencyList_type>::out_edge_iterator vi, vi_end; graph_traits<AdjacencyList_type>::out_edge_iterator vi, vi_end;
for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi) 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 if (has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
circuit_stack.push_back(v_index[target(*vi, g)]); {
return true; // cycle detected, return immediately
circuit_stack.push_back(v_index[target(*vi, g)]);
return true;
}
} }
else if (color[target(*vi, g)] == gray_color) else if (color[target(*vi, g)] == gray_color)
{ {
...@@ -209,6 +215,8 @@ namespace MFS ...@@ -209,6 +215,8 @@ namespace MFS
vector_vertex_descriptor vector_vertex_descriptor
Collect_Doublet(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G) 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::in_edge_iterator it_in, in_end;
AdjacencyList_type::out_edge_iterator it_out, out_end; AdjacencyList_type::out_edge_iterator it_out, out_end;
vector<AdjacencyList_type::vertex_descriptor> Doublet; vector<AdjacencyList_type::vertex_descriptor> Doublet;
...@@ -223,6 +231,7 @@ namespace MFS ...@@ -223,6 +231,7 @@ namespace MFS
bool bool
Vertex_Belong_to_a_Clique(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G) 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<AdjacencyList_type::vertex_descriptor> liste; vector<AdjacencyList_type::vertex_descriptor> liste;
bool agree = true; bool agree = true;
AdjacencyList_type::in_edge_iterator it_in, in_end; AdjacencyList_type::in_edge_iterator it_in, in_end;
...@@ -262,6 +271,7 @@ namespace MFS ...@@ -262,6 +271,7 @@ namespace MFS
bool bool
Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G) 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 something_has_been_done = false;
bool not_a_loop; bool not_a_loop;
int i; int i;
......
...@@ -65,6 +65,20 @@ ModelTree::computeJacobian(const set<int> &vars) ...@@ -65,6 +65,20 @@ ModelTree::computeJacobian(const set<int> &vars)
} }
} }
void
ModelTree::writeChaineRuleDerivative(ostream &output, int eq, int var, int lag,
ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
map<pair<int, pair<int, int> >, 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 void
ModelTree::computeHessian(const set<int> &vars) ModelTree::computeHessian(const set<int> &vars)
{ {
......
...@@ -47,6 +47,9 @@ protected: ...@@ -47,6 +47,9 @@ protected:
*/ */
first_derivatives_type first_derivatives; first_derivatives_type first_derivatives;