From a67a20bf668c3a912c8baba7e525068aa70ef10d Mon Sep 17 00:00:00 2001 From: ferhat Date: Fri, 14 Nov 2008 16:07:47 +0000 Subject: [PATCH] - new Incidence_Matrix class - lead and lag on exogenous variables - corrections in dr1_sparse and dr11_sparse - minor corrections in simulate git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2255 ac1d8469-bf42-47a9-8791-bf33cf982152 --- BlockTriangular.cc | 928 +++++++++++++++++++++++++------------ ComputingTasks.cc | 6 +- ExprNode.cc | 49 +- ModFile.cc | 2 +- ModelTree.cc | 764 +++++++++++++++++------------- ParsingDriver.cc | 18 +- include/BlockTriangular.hh | 56 ++- include/ExprNode.hh | 20 +- include/ModelTree.hh | 5 +- 9 files changed, 1177 insertions(+), 671 deletions(-) diff --git a/BlockTriangular.cc b/BlockTriangular.cc index b21c48c1..205122d7 100644 --- a/BlockTriangular.cc +++ b/BlockTriangular.cc @@ -31,158 +31,346 @@ using namespace std; #include "BlockTriangular.hh" //------------------------------------------------------------------------------ +/*BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) : + symbol_table(symbol_table_arg), + normalization(symbol_table_arg)*/ BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) : symbol_table(symbol_table_arg), - normalization(symbol_table_arg) + normalization(symbol_table_arg), + incidencematrix(symbol_table_arg) { bt_verbose = 0; ModelBlock = NULL; - Model_Max_Lead = 0; - Model_Max_Lag = 0; periods = 0; } + +IncidenceMatrix::IncidenceMatrix(const SymbolTable &symbol_table_arg) : + symbol_table(symbol_table_arg) +{ + Model_Max_Lead = Model_Max_Lead_Endo = Model_Max_Lead_Exo = 0; + Model_Max_Lag = Model_Max_Lag_Endo = Model_Max_Lag_Exo = 0; + +} //------------------------------------------------------------------------------ //For a lead or a lag build the Incidence Matrix structures List_IM* -BlockTriangular::Build_IM(int lead_lag) +IncidenceMatrix::Build_IM(int lead_lag, SymbolType type) { - List_IM* pIM = new List_IM; int i; - Last_IM->pNext = pIM; - pIM->IM = (bool*)malloc(endo_nbr * endo_nbr * sizeof(pIM->IM[0])); - for(i = 0;i < endo_nbr*endo_nbr;i++) - pIM->IM[i] = 0; - pIM->lead_lag = lead_lag; - if(lead_lag > 0) + List_IM* pIM = new List_IM; + if(type==eEndogenous) { - if(lead_lag > Model_Max_Lead) - Model_Max_Lead = lead_lag; + Last_IM->pNext = pIM; + pIM->IM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(pIM->IM[0])); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) + pIM->IM[i] = 0; + pIM->lead_lag = lead_lag; + if(lead_lag > 0) + { + if(lead_lag > Model_Max_Lead_Endo) + { + Model_Max_Lead_Endo = lead_lag; + if(lead_lag > Model_Max_Lead) + Model_Max_Lead = lead_lag; + } + } + else + { + if( -lead_lag > Model_Max_Lag_Endo) + { + Model_Max_Lag_Endo = -lead_lag; + if(-lead_lag > Model_Max_Lag) + Model_Max_Lag = -lead_lag; + } + } + pIM->pNext = NULL; + Last_IM = pIM; } else - { - if( -lead_lag > Model_Max_Lag) - Model_Max_Lag = -lead_lag; + { //eExogenous + Last_IM_X->pNext = pIM; + pIM->IM = (bool*)malloc(symbol_table.exo_nbr * symbol_table.endo_nbr * sizeof(pIM->IM[0])); + for(i = 0;i < symbol_table.exo_nbr*symbol_table.endo_nbr;i++) + pIM->IM[i] = 0; + pIM->lead_lag = lead_lag; + if(lead_lag > 0) + { + if(lead_lag > Model_Max_Lead_Exo) + { + Model_Max_Lead_Exo = lead_lag; + if(lead_lag > Model_Max_Lead) + Model_Max_Lead = lead_lag; + } + } + else + { + if( -lead_lag > Model_Max_Lag_Exo) + { + Model_Max_Lag_Exo = -lead_lag; + if(-lead_lag > Model_Max_Lag) + Model_Max_Lag = -lead_lag; + } + } + pIM->pNext = NULL; + Last_IM_X = pIM; } - pIM->pNext = NULL; - Last_IM = pIM; return (pIM); } //------------------------------------------------------------------------------ // initialize all the incidence matrix structures void -BlockTriangular::init_incidence_matrix(int nb_endo) +IncidenceMatrix::init_incidence_matrix() { int i; - endo_nbr = nb_endo; First_IM = new List_IM; - First_IM->IM = (bool*)malloc(nb_endo * nb_endo * sizeof(First_IM->IM[0])); - for(i = 0;i < nb_endo*nb_endo;i++) + First_IM->IM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(First_IM->IM[0])); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) First_IM->IM[i] = 0; First_IM->lead_lag = 0; First_IM->pNext = NULL; Last_IM = First_IM; - //cout << "init_incidence_matrix done \n"; + + First_IM_X = new List_IM; + First_IM_X->IM = (bool*)malloc(symbol_table.exo_nbr * symbol_table.endo_nbr * sizeof(First_IM_X->IM[0])); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.exo_nbr;i++) + First_IM_X->IM[i] = 0; + First_IM_X->lead_lag = 0; + First_IM_X->pNext = NULL; + Last_IM_X = First_IM_X; + } void -BlockTriangular::Free_IM(List_IM* First_IM) const +IncidenceMatrix::Free_IM() const { -#ifdef DEBUG - cout << "Free_IM\n"; -#endif List_IM *Cur_IM, *SFirst_IM; Cur_IM = SFirst_IM = First_IM; while(Cur_IM) { - First_IM = Cur_IM->pNext; + SFirst_IM = Cur_IM->pNext; free(Cur_IM->IM); delete Cur_IM; - Cur_IM = First_IM; + Cur_IM = SFirst_IM; + } + Cur_IM = SFirst_IM = First_IM_X; + while(Cur_IM) + { + SFirst_IM = Cur_IM->pNext; + free(Cur_IM->IM); + delete Cur_IM; + Cur_IM = SFirst_IM; } - //free(SFirst_IM); - //delete SFirst_IM; } //------------------------------------------------------------------------------ -// Return the inceidence matrix related to a lead or a lag +// Return the incidence matrix related to a lead or a lag List_IM* -BlockTriangular::Get_IM(int lead_lag) +IncidenceMatrix::Get_IM(int lead_lag, SymbolType type) const { List_IM* Cur_IM; - Cur_IM = First_IM; + if(type==eEndogenous) + Cur_IM = First_IM; + else + Cur_IM = First_IM_X; while ((Cur_IM != NULL) && (Cur_IM->lead_lag != lead_lag)) Cur_IM = Cur_IM->pNext; return (Cur_IM); } bool* -BlockTriangular::bGet_IM(int lead_lag) const +IncidenceMatrix::bGet_IM(int lead_lag, SymbolType type) const { List_IM* Cur_IM; - Cur_IM = First_IM; + if(type==eEndogenous) + Cur_IM = First_IM; + else + Cur_IM = First_IM_X; while ((Cur_IM != NULL) && (Cur_IM->lead_lag != lead_lag)) { Cur_IM = Cur_IM->pNext; } - if((Cur_IM->lead_lag != lead_lag) || (Cur_IM==NULL)) - { - cout << "the incidence matrix with lag " << lead_lag << " does not exist !!"; - exit(EXIT_FAILURE); - } - return (Cur_IM->IM); + if(Cur_IM) + return (Cur_IM->IM); + else + return(0); } //------------------------------------------------------------------------------ // Fill the incidence matrix related to a lead or a lag void -BlockTriangular::fill_IM(int equation, int variable_endo, int lead_lag) +IncidenceMatrix::fill_IM(int equation, int variable, int lead_lag, SymbolType type) { List_IM* Cur_IM; - //cout << "equation=" << equation << " variable_endo=" << variable_endo << " lead_lag=" << lead_lag << "\n"; - Cur_IM = Get_IM(lead_lag); - if(equation >= endo_nbr) + Cur_IM = Get_IM(lead_lag, type); + if(equation >= symbol_table.endo_nbr) { - cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n"; - system("PAUSE"); + cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << symbol_table.endo_nbr << ")\n"; exit(EXIT_FAILURE); } if (!Cur_IM) - Cur_IM = Build_IM(lead_lag); - Cur_IM->IM[equation*endo_nbr + variable_endo] = 1; + Cur_IM = Build_IM(lead_lag, type); + if(type==eEndogenous) + Cur_IM->IM[equation*symbol_table.endo_nbr + variable] = 1; + else + Cur_IM->IM[equation*symbol_table.exo_nbr + variable] = 1; } //------------------------------------------------------------------------------ // unFill the incidence matrix related to a lead or a lag void -BlockTriangular::unfill_IM(int equation, int variable_endo, int lead_lag) +IncidenceMatrix::unfill_IM(int equation, int variable, int lead_lag, SymbolType type) { List_IM* Cur_IM; - //cout << "lead_lag=" << lead_lag << "\n"; - Cur_IM = Get_IM(lead_lag); - /*if(equation >= endo_nbr) - { - cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n"; - system("PAUSE"); - exit(EXIT_FAILURE); - }*/ + Cur_IM = Get_IM(lead_lag, type); if (!Cur_IM) - Cur_IM = Build_IM(lead_lag); - Cur_IM->IM[equation*endo_nbr + variable_endo] = 0; - /*system("pause");*/ + Cur_IM = Build_IM(lead_lag, type); + if(type==eEndogenous) + Cur_IM->IM[equation*symbol_table.endo_nbr + variable] = 0; + else + Cur_IM->IM[equation*symbol_table.exo_nbr + variable] = 0; +} + +List_IM* +IncidenceMatrix::Get_First(SymbolType type) const +{ + if(type==eEndogenous) + return(First_IM); + else + return(First_IM_X); } +//------------------------------------------------------------------------------ +//For a lead or a lag build the Incidence Matrix structures +/*List_IM* +IncidenceMatrix::Build_IM_X(int lead_lag) +{ + List_IM* pIM_X = new List_IM; + int i; + Last_IM_X->pNext = pIM_X; + pIM_X->IM = (bool*)malloc(exo_nbr * endo_nbr * sizeof(pIM_X->IM[0])); + for(i = 0;i < exo_nbr*endo_nbr;i++) + pIM_X->IM[i] = 0; + pIM_X->lead_lag = lead_lag; + if(lead_lag > 0) + { + if(lead_lag > Model_Max_Lead_Exo) + { + Model_Max_Lead_Exo = lead_lag; + if(lead_lag > Model_Max_Lead) + Model_Max_Lead = lead_lag; + } + } + else + { + if( -lead_lag > Model_Max_Lag_Exo) + { + Model_Max_Lag_Exo = -lead_lag; + if(-lead_lag > Model_Max_Lag) + Model_Max_Lag = -lead_lag; + } + } + pIM_X->pNext = NULL; + Last_IM_X = pIM_X; + return (pIM_X); +} + +//------------------------------------------------------------------------------ +// initialize all the incidence matrix structures +void +IncidenceMatrix::init_incidence_matrix_X(int nb_exo) +{ + int i; + //cout << "init_incidence_matrix_X nb_exo = " << nb_exo << " endo_nbr=" << endo_nbr << "\n"; + exo_nbr = nb_exo; + First_IM_X = new List_IM; + First_IM_X->IM = (bool*)malloc(nb_exo * endo_nbr * sizeof(First_IM_X->IM[0])); + for(i = 0;i < endo_nbr*nb_exo;i++) + First_IM_X->IM[i] = 0; + First_IM_X->lead_lag = 0; + First_IM_X->pNext = NULL; + Last_IM_X = First_IM_X; +} + + +void +IncidenceMatrix::Free_IM_X(List_IM* First_IM_X) const +{ + List_IM *Cur_IM_X, *SFirst_IM_X; + Cur_IM_X = SFirst_IM_X = First_IM_X; + while(Cur_IM_X) + { + First_IM_X = Cur_IM_X->pNext; + free(Cur_IM_X->IM); + delete Cur_IM_X; + Cur_IM_X = First_IM_X; + } +} + + +//------------------------------------------------------------------------------ +// Return the inceidence matrix related to a lead or a lag +List_IM* +IncidenceMatrix::Get_IM_X(int lead_lag) +{ + List_IM* Cur_IM_X; + Cur_IM_X = First_IM_X; + while ((Cur_IM_X != NULL) && (Cur_IM_X->lead_lag != lead_lag)) + Cur_IM_X = Cur_IM_X->pNext; + return (Cur_IM_X); +} + +bool* +IncidenceMatrix::bGet_IM_X(int lead_lag) const +{ + List_IM* Cur_IM_X; + Cur_IM_X = First_IM_X; + while ((Cur_IM_X != NULL) && (Cur_IM_X->lead_lag != lead_lag)) + { + Cur_IM_X = Cur_IM_X->pNext; + } + if(Cur_IM_X) + return (Cur_IM_X->IM); + else + return(0); +} + + +//------------------------------------------------------------------------------ +// Fill the incidence matrix related to a lead or a lag +void +IncidenceMatrix::fill_IM_X(int equation, int variable_exo, int lead_lag) +{ + List_IM* Cur_IM_X; + Cur_IM_X = Get_IM_X(lead_lag); + if(equation >= endo_nbr) + { + cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n"; + exit(EXIT_FAILURE); + } + if (!Cur_IM_X) + { + Cur_IM_X = Build_IM_X(lead_lag); + } + Cur_IM_X->IM[equation*exo_nbr + variable_exo] = true; +} +*/ + //------------------------------------------------------------------------------ //Print azn incidence matrix void -BlockTriangular::Print_SIM(bool* IM, int n) const +IncidenceMatrix::Print_SIM(bool* IM, SymbolType type) const { - int i, j; - for(i = 0;i < n;i++) + int i, j, n; + if(type == eEndogenous) + n = symbol_table.endo_nbr; + else + n = symbol_table.exo_nbr; + for(i = 0;i < symbol_table.endo_nbr;i++) { cout << " "; for(j = 0;j < n;j++) @@ -194,15 +382,18 @@ BlockTriangular::Print_SIM(bool* IM, int n) const //------------------------------------------------------------------------------ //Print all incidence matrix void -BlockTriangular::Print_IM(int n) const +IncidenceMatrix::Print_IM(SymbolType type) const { List_IM* Cur_IM; - Cur_IM = First_IM; + if(type == eEndogenous) + Cur_IM = First_IM; + else + Cur_IM = First_IM_X; cout << "-------------------------------------------------------------------\n"; while(Cur_IM) { cout << "Incidence matrix for lead_lag = " << Cur_IM->lead_lag << "\n"; - Print_SIM(Cur_IM->IM, n); + Print_SIM(Cur_IM->IM, type); Cur_IM = Cur_IM->pNext; } } @@ -211,7 +402,7 @@ BlockTriangular::Print_IM(int n) const //------------------------------------------------------------------------------ // Swap rows and columns of the incidence matrix void -BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) +IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const { int tmp_i, j; bool tmp_b; @@ -269,15 +460,13 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n if ((k == 1) && IM0[Index_Equ_IM[i].index*n + Index_Var_IM[l].index]) { modifie = 1; - swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n); + incidencematrix.swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n); (*prologue)++; } } } *epilogue = 0; modifie = 1; - /*print_SIM(IM,n); - print_SIM(IM*/ while(modifie) { modifie = 0; @@ -295,7 +484,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n if ((k == 1) && IM0[Index_Equ_IM[l].index*n + Index_Var_IM[i].index]) { modifie = 1; - swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n); + incidencematrix.swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n); (*epilogue)++; } } @@ -306,11 +495,12 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n void BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock) { - int i, j, k, l, ls, m, i_1, Lead, Lag, size_list_lead_var, first_count_equ, i1; - int *list_lead_var, *tmp_size, *tmp_var, *tmp_endo, nb_lead_lag_endo; + int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1; + int *tmp_size, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_exo, tmp_nb_exo, nb_lead_lag_endo; List_IM *Cur_IM; bool *IM, OK; ModelBlock->Periods = periods; + int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo; if ((type == PROLOGUE) || (type == EPILOGUE)) { for(i = 0;i < size;i++) @@ -321,128 +511,221 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN; ModelBlock->Block_List[*count_Block].Temporary_terms=new temporary_terms_type (); ModelBlock->Block_List[*count_Block].Temporary_terms->clear(); - list_lead_var = (int*)malloc(Model_Max_Lead * endo_nbr * sizeof(int)); - size_list_lead_var = 0; - tmp_endo = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); - tmp_size = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); + tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); tmp_var = (int*)malloc(sizeof(int)); - memset(tmp_size, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); - memset(tmp_endo, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); + tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); nb_lead_lag_endo = Lead = Lag = 0; - Cur_IM = First_IM; + Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0; + Cur_IM = incidencematrix.Get_First(eEndogenous); while(Cur_IM) { k = Cur_IM->lead_lag; - i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr; + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr; if(k > 0) { - if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index]) + if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index]) { nb_lead_lag_endo++; - tmp_size[Model_Max_Lag + k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if(k > Lead) - { - Lead = k; - list_lead_var[size_list_lead_var] = Index_Var_IM[*count_Equ].index + size * (k - 1); - size_list_lead_var++; - } + Lead = k; } } else { k = -k; - if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index]) + if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index]) { - tmp_size[Model_Max_Lag - k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; nb_lead_lag_endo++; if(k > Lag) - { - Lag = k; - } + Lag = k; } } Cur_IM = Cur_IM->pNext; } + + Lag_Endo = Lag; + Lead_Endo = Lead; + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + tmp_nb_exo = 0; + + + Cur_IM = incidencematrix.Get_First(eExogenous); + k = Cur_IM->lead_lag; + while(Cur_IM) + { + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(j=0;jIM[i_1 + j]) + { + if(!tmp_exo[j]) + { + tmp_exo[j] = 1; + tmp_nb_exo++; + } + if(k>0 && k>Lead_Exo) + Lead_Exo = k; + else if(k<0 && (-k)>Lag_Exo) + Lag_Exo = -k; + if(k>0 && k>Lead) + Lead = k; + else if(k<0 && (-k)>Lag) + Lag = -k; + tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++; + } + Cur_IM = Cur_IM->pNext; + if(Cur_IM) + k = Cur_IM->lead_lag; + } + + + ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo; + ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); + k = 0; + for(j=0;jBlock_List[*count_Block].Exogenous[k] = j; + k++; + } + + ModelBlock->Block_List[*count_Block].nb_exo_det = 0; + ModelBlock->Block_List[*count_Block].Max_Lag = Lag; ModelBlock->Block_List[*count_Block].Max_Lead = Lead; - free(list_lead_var); + ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo; + ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo; + ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo; + ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo; ModelBlock->Block_List[*count_Block].Equation = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Variable = (int*)malloc(sizeof(int)); - ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Equation[0] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].Variable[0] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].Variable_Sorted[0] = -1; - ModelBlock->in_Block_Equ[Index_Equ_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = 0; if ((Lead > 0) && (Lag > 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE; else if((Lead > 0) && (Lag == 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_BACKWARD_SIMPLE; else ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE; + + Cur_IM = incidencematrix.Get_First(eExogenous); + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + tmp_nb_exo = 0; + while(Cur_IM) + { + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(j=0;jIM[i_1 + j] && (!tmp_exo[j])) + { + tmp_exo[j] = 1; + tmp_nb_exo++; + } + Cur_IM = Cur_IM->pNext; + } + ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo; + ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; - if(nb_lead_lag_endo) - { - ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - } + + + k = 0; + for(j=0;jBlock_List[*count_Block].Exogenous[k] = j; + k++; + } ls = l = 1; i1 = 0; for(int li = 0;li < Lead + Lag + 1;li++) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[Model_Max_Lag - Lag + li]; - if(tmp_size[Model_Max_Lag - Lag + li]) + if(incidencematrix.Model_Max_Lag_Endo - Lag + li>=0) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[Model_Max_Lag - Lag + li]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].size=" << tmp_size[Model_Max_Lag_Endo - Lag + li] << "\n"; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l; - IM = bGet_IM(li - Lag); - if(IM == NULL) + IM = incidencematrix.bGet_IM(li - Lag, eEndogenous); + if(IM) { - cout << "Error IM(" << li - Lag << ") doesn't exist\n"; - exit(EXIT_FAILURE); - } - if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*endo_nbr] && nb_lead_lag_endo) - { - ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = li - Lag; - tmp_var[0] = i1; - i1++; - } - m = 0; - i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr; - if(IM[Index_Var_IM[*count_Equ].index + i_1]) - { - if(li == Lag) + if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*symbol_table.endo_nbr] && nb_lead_lag_endo) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls; - ls++; + tmp_var[0] = i1; + i1++; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[0]]; - l++; - m++; + m = 0; + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr; + if(IM[Index_Var_IM[*count_Equ].index + i_1]) + { + if(li == Lag) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls; + ls++; + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index; + l++; + m++; + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1; - } + } + else + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = 0; + if(incidencematrix.Model_Max_Lag_Exo - Lag + li>=0) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li]; + //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].Exogenous= malloc(" << tmp_size_exo[Model_Max_Lag_Exo - Lag + li] << ")\n"; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + IM = incidencematrix.bGet_IM(li - Lag, eExogenous); + if(IM) + { + m = 0; + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(k = 0; kBlock_List[*count_Block].Exogenous[k]+i_1]) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous[m] = k; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index[m] = ModelBlock->Block_List[*count_Block].Exogenous[k]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X_Index[m] = Index_Equ_IM[*count_Equ].index; + m++; + } + } + } + } + else + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = 0; } (*count_Equ)++; (*count_Block)++; free(tmp_size); + free(tmp_size_exo); free(tmp_endo); + free(tmp_exo); free(tmp_var); } } @@ -456,25 +739,25 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN; ModelBlock->Block_List[*count_Block].Equation = (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].Variable_Sorted = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); + //ModelBlock->Block_List[*count_Block].Variable_Sorted = (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)); Lead = Lag = 0; first_count_equ = *count_Equ; tmp_var = (int*)malloc(size * sizeof(int)); - tmp_endo = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); - tmp_size = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); - memset(tmp_size, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); - memset(tmp_endo, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); + tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); nb_lead_lag_endo = 0; + Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0; + for(i = 0;i < size;i++) { ModelBlock->Block_List[*count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].Variable[i] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].Variable_Sorted[i] = -1; - ModelBlock->in_Block_Equ[Index_Equ_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = i; - Cur_IM = First_IM; + Cur_IM = incidencematrix.Get_First(eEndogenous); i_1 = Index_Var_IM[*count_Equ].index; while(Cur_IM) { @@ -484,12 +767,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { for(j = 0;j < size;j++) { - if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*endo_nbr]) + if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr]) { - tmp_size[Model_Max_Lag + k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if (!OK) { - tmp_endo[Model_Max_Lag + k]++; + tmp_endo[incidencematrix.Model_Max_Lag + k]++; nb_lead_lag_endo++; OK = true; } @@ -503,12 +786,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc k = -k; for(j = 0;j < size;j++) { - if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*endo_nbr]) + if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr]) { - tmp_size[Model_Max_Lag - k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; if (!OK) { - tmp_endo[Model_Max_Lag - k]++; + tmp_endo[incidencematrix.Model_Max_Lag - k]++; nb_lead_lag_endo++; OK = true; } @@ -537,87 +820,158 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc else ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE; } + + + Lag_Endo = Lag; + Lead_Endo = Lead; + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + tmp_nb_exo = 0; + for(i = 0;i < size;i++) + { + Cur_IM = incidencematrix.Get_First(eExogenous); + k = Cur_IM->lead_lag; + while(Cur_IM) + { + i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.exo_nbr; + for(j=0;jIM[i_1 + j]) + { + if(!tmp_exo[j]) + { + tmp_exo[j] = 1; + tmp_nb_exo++; + } + if(k>0 && k>Lead_Exo) + Lead_Exo = k; + else if(k<0 && (-k)>Lag_Exo) + Lag_Exo = -k; + if(k>0 && k>Lead) + Lead = k; + else if(k<0 && (-k)>Lag) + Lag = -k; + tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++; + } + Cur_IM = Cur_IM->pNext; + if(Cur_IM) + k = Cur_IM->lead_lag; + } + } + + + ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo; + ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); + k = 0; + for(j=0;jBlock_List[*count_Block].Exogenous[k] = j; + k++; + } + + ModelBlock->Block_List[*count_Block].nb_exo_det = 0; + ModelBlock->Block_List[*count_Block].Max_Lag = Lag; ModelBlock->Block_List[*count_Block].Max_Lead = Lead; + ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo; + ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo; + ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo; + ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo; ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); ls = l = size; i1 = 0; ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; - ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int)); for(i = 0;i < Lead + Lag + 1;i++) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[Model_Max_Lag - Lag + i]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_endo[Model_Max_Lag - Lag + i]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l; - IM = bGet_IM(i - Lag); - if(IM != NULL) + if(incidencematrix.Model_Max_Lag_Endo-Lag+i>=0) { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); } else + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = 0; + if(incidencematrix.Model_Max_Lag_Exo-Lag+i>=0) { - cout << "Error IM(" << i - Lag << ") doesn't exist\n"; - exit(EXIT_FAILURE); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); } - for(j = first_count_equ;j < size + first_count_equ;j++) + else + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l; + IM = incidencematrix.bGet_IM(i - Lag, eEndogenous); + if(IM) { - i_1 = Index_Var_IM[j].index; - m = 0; - for(k = first_count_equ;k < size + first_count_equ;k++) - if(IM[i_1 + Index_Equ_IM[k].index*endo_nbr]) - m++; - if(m > 0) + for(j = first_count_equ;j < size + first_count_equ;j++) { - ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = i_1; - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = i - Lag; - tmp_var[j - first_count_equ] = i1; - i1++; + i_1 = Index_Var_IM[j].index; + m = 0; + for(k = first_count_equ;k < size + first_count_equ;k++) + if(IM[i_1 + Index_Equ_IM[k].index*symbol_table.endo_nbr]) + m++; + if(m > 0) + { + tmp_var[j - first_count_equ] = i1; + i1++; + } + } + m = 0; + for(j = first_count_equ;j < size + first_count_equ;j++) + { + i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr; + for(k = first_count_equ;k < size + first_count_equ;k++) + if(IM[Index_Var_IM[k].index + i_1]) + { + if(i == Lag) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls; + ls++; + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index; + l++; + m++; + } + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1; + } + IM = incidencematrix.bGet_IM(i - Lag, eExogenous); + if(IM) + { + m = 0; + for(j = first_count_equ;j < size + first_count_equ;j++) + { + i_1 = Index_Equ_IM[j].index * symbol_table.exo_nbr; + for(k = 0; kBlock_List[*count_Block].Exogenous[k]+i_1]) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous[m] = k; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous_Index[m] = ModelBlock->Block_List[*count_Block].Exogenous[k]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X[m] = j - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X_Index[m] = Index_Equ_IM[j].index; + m++; + } + } } } - m = 0; - for(j = first_count_equ;j < size + first_count_equ;j++) - { - i_1 = Index_Equ_IM[j].index * endo_nbr; - for(k = first_count_equ;k < size + first_count_equ;k++) - if(IM[Index_Var_IM[k].index + i_1]) - { - if(i == Lag) - { - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls; - ls++; -#ifdef DEBUG - printf("j=%d ModelBlock->Block_List[%d].Variable[%d]=%d Index_Var_IM[%d].index=%d", j, *count_Block, j - first_count_equ, ModelBlock->Block_List[*count_Block].Variable[j - first_count_equ], k, Index_Var_IM[k].index); - if(ModelBlock->Block_List[*count_Block].Variable[j - first_count_equ] == Index_Var_IM[k].index) - { - ModelBlock->Block_List[*count_Block].Own_Derivative[j - first_count_equ]=l; - printf(" l=%d OK\n",l); - } - else - printf("\n"); -#endif - } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[k - first_count_equ]]; - l++; - m++; - } - } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1; } (*count_Block)++; free(tmp_size); + free(tmp_size_exo); free(tmp_endo); + free(tmp_exo); free(tmp_var); } } @@ -627,11 +981,6 @@ void BlockTriangular::Free_Block(Model_Block* ModelBlock) const { int blk, i; -#ifdef DEBUG - - cout << "Free_Block\n"; -#endif - for(blk = 0;blk < ModelBlock->Size;blk++) { @@ -639,22 +988,26 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const { free(ModelBlock->Block_List[blk].Equation); free(ModelBlock->Block_List[blk].Variable); - free(ModelBlock->Block_List[blk].Variable_Sorted); + free(ModelBlock->Block_List[blk].Exogenous); free(ModelBlock->Block_List[blk].Own_Derivative); - if(ModelBlock->Block_List[blk].Nb_Lead_Lag_Endo) - { - free(ModelBlock->Block_List[blk].variable_dyn_index); - free(ModelBlock->Block_List[blk].variable_dyn_leadlag); - } for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) { - free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index); + if(ModelBlock->Block_List[blk].IM_lead_lag[i].size) + { + free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); + } + if(ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo) + { + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X); + } } free(ModelBlock->Block_List[blk].IM_lead_lag); delete(ModelBlock->Block_List[blk].Temporary_terms); @@ -663,28 +1016,31 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const { free(ModelBlock->Block_List[blk].Equation); free(ModelBlock->Block_List[blk].Variable); - free(ModelBlock->Block_List[blk].Variable_Sorted); + free(ModelBlock->Block_List[blk].Exogenous); free(ModelBlock->Block_List[blk].Own_Derivative); - free(ModelBlock->Block_List[blk].variable_dyn_index); - free(ModelBlock->Block_List[blk].variable_dyn_leadlag); for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) { - free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index); + if(incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 && ModelBlock->Block_List[blk].IM_lead_lag[i].size) + { + free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); + } + if(incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 && ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo) + { + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X); + } } free(ModelBlock->Block_List[blk].IM_lead_lag); delete(ModelBlock->Block_List[blk].Temporary_terms); } } - free(ModelBlock->in_Block_Equ); - free(ModelBlock->in_Block_Var); - free(ModelBlock->in_Equ_of_Block); - free(ModelBlock->in_Var_of_Block); free(ModelBlock->Block_List); free(ModelBlock); free(Index_Equ_IM); @@ -700,40 +1056,19 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int i, j, Nb_TotalBlocks, Nb_RecursBlocks; int count_Block, count_Equ; block_result_t* res; - //List_IM * p_First_IM, *p_Cur_IM, *Cur_IM; Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set)); bool* SIM0, *SIM00; - /*p_First_IM = (List_IM*)malloc(sizeof(*p_First_IM)); - p_Cur_IM = p_First_IM; - Cur_IM = First_IM; - i = endo_nbr * endo_nbr; - while(Cur_IM) - { - p_Cur_IM->lead_lag = Cur_IM->lead_lag; - p_Cur_IM->IM = (bool*)malloc(i * sizeof(bool)); - memcpy ( p_Cur_IM->IM, Cur_IM->IM, i); - Cur_IM = Cur_IM->pNext; - if(Cur_IM) - { - p_Cur_IM->pNext = (List_IM*)malloc(sizeof(*p_Cur_IM)); - p_Cur_IM = p_Cur_IM->pNext; - } - else - p_Cur_IM->pNext = NULL; - }*/ SIM0 = (bool*)malloc(n * n * sizeof(bool)); - //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n"; memcpy(SIM0,IM_0,n*n*sizeof(bool)); Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0); - //cout << "free SIM0=" << SIM0 << "\n"; free(SIM0); if(bt_verbose) { cout << "prologue : " << *prologue << " epilogue : " << *epilogue << "\n"; cout << "IM_0\n"; - Print_SIM(IM_0, n); + incidencematrix.Print_SIM(IM_0, eEndogenous); cout << "IM\n"; - Print_SIM(IM, n); + incidencematrix.Print_SIM(IM, eEndogenous); for(i = 0;i < n;i++) cout << "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index << " Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i].index << "\n"; } @@ -745,11 +1080,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, if(mixing) { double* max_val=(double*)malloc(n*sizeof(double)); - //cout << "n=" << n << "\n"; memset(max_val,0,n*sizeof(double)); for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) { - //cout << "iter->first.first=" << iter->first.first << "\n"; if(fabs(iter->second)>max_val[iter->first.first]) max_val[iter->first.first]=fabs(iter->second); } @@ -763,17 +1096,13 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, { int suppress=0; SIM0 = (bool*)malloc(n * n * sizeof(bool)); - //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n"; memset(SIM0,0,n*n*sizeof(bool)); SIM00 = (bool*)malloc(n * n * sizeof(bool)); - //cout << "Allocate SIM00=" << SIM00 << " size=" << n * n * sizeof(bool) << "\n"; memset(SIM00,0,n*n*sizeof(bool)); - //cout << "n*n=" << n*n << "\n"; for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) { if(fabs(iter->second)>bi) { - //cout << "iter->first.first*n+iter->first.second=" << iter->first.first*n+iter->first.second << "\n"; SIM0[iter->first.first*n+iter->first.second]=1; if(!IM_0[iter->first.first*n+iter->first.second]) { @@ -783,14 +1112,11 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, else suppress++; } - //cout << "n*n=" << n*n << "\n"; for(i = 0;i < n;i++) for(j = 0;j < n;j++) { - //cout << "Index_Equ_IM[i].index * n + Index_Var_IM[j].index=" << Index_Equ_IM[i].index * n + Index_Var_IM[j].index << "\n"; SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index]; } - //cout << "free SIM0=" << SIM0 << "\n"; free(SIM0); if(suppress!=suppressed) { @@ -850,15 +1176,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, cout << " the largest simultaneous block has " << j << " equation(s). \n"; ModelBlock->Size = Nb_TotalBlocks; ModelBlock->Periods = periods; - ModelBlock->in_Block_Equ = (int*)malloc(n * sizeof(int)); - ModelBlock->in_Block_Var = (int*)malloc(n * sizeof(int)); - ModelBlock->in_Equ_of_Block = (int*)malloc(n * sizeof(int)); - ModelBlock->in_Var_of_Block = (int*)malloc(n * sizeof(int)); ModelBlock->Block_List = (Block*)malloc(sizeof(ModelBlock->Block_List[0]) * Nb_TotalBlocks); - blocks.block_result_to_IM(res, IM, *prologue, endo_nbr, Index_Equ_IM, Index_Var_IM); - //Free_IM(p_First_IM); + blocks.block_result_to_IM(res, IM, *prologue, symbol_table.endo_nbr, Index_Equ_IM, Index_Var_IM); count_Equ = count_Block = 0; - //Print_IM(endo_nbr); if (*prologue) Allocate_Block(*prologue, &count_Equ, &count_Block, PROLOGUE, ModelBlock); for(j = 0;j < res->n_sets;j++) @@ -887,11 +1207,11 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ List_IM* Cur_IM; int i; //First create a static model incidence matrix - SIM = (bool*)malloc(endo_nbr * endo_nbr * sizeof(*SIM)); - for(i = 0;i < endo_nbr*endo_nbr;i++) + SIM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(*SIM)); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) { SIM[i] = 0; - Cur_IM = First_IM; + Cur_IM = incidencematrix.Get_First(eEndogenous); while(Cur_IM) { SIM[i] = (SIM[i]) || (Cur_IM->IM[i]); @@ -901,26 +1221,26 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ if(bt_verbose) { cout << "incidence matrix for the static model (unsorted) \n"; - Print_SIM(SIM, endo_nbr); + incidencematrix.Print_SIM(SIM, eEndogenous); } - Index_Equ_IM = (simple*)malloc(endo_nbr * sizeof(*Index_Equ_IM)); - for(i = 0;i < endo_nbr;i++) + Index_Equ_IM = (simple*)malloc(symbol_table.endo_nbr * sizeof(*Index_Equ_IM)); + for(i = 0;i < symbol_table.endo_nbr;i++) { Index_Equ_IM[i].index = i; } - Index_Var_IM = (simple*)malloc(endo_nbr * sizeof(*Index_Var_IM)); - for(i = 0;i < endo_nbr;i++) + Index_Var_IM = (simple*)malloc(symbol_table.endo_nbr * sizeof(*Index_Var_IM)); + for(i = 0;i < symbol_table.endo_nbr;i++) { Index_Var_IM[i].index = i; } if(ModelBlock != NULL) Free_Block(ModelBlock); ModelBlock = (Model_Block*)malloc(sizeof(*ModelBlock)); - Cur_IM = Get_IM(0); - SIM_0 = (bool*)malloc(endo_nbr * endo_nbr * sizeof(*SIM_0)); - for(i = 0;i < endo_nbr*endo_nbr;i++) + Cur_IM = incidencematrix.Get_IM(0, eEndogenous); + 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->IM[i]; - Normalize_and_BlockDecompose(SIM, ModelBlock, endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m); + Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m); free(SIM_0); free(SIM); } diff --git a/ComputingTasks.cc b/ComputingTasks.cc index c886f83b..a56f953e 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -169,7 +169,11 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename) const { options_list.writeOutput(output); symbol_list.writeOutput("var_list_", output); - output << "info = stoch_simul(var_list_);\n"; + output << "if(options_.model_mode)\n"; + output << " info = stoch_simul_sparse(var_list_);\n"; + output << "else\n"; + output << " info = stoch_simul(var_list_);\n"; + output << "end\n"; } ForecastStatement::ForecastStatement(const SymbolList &symbol_list_arg, diff --git a/ExprNode.cc b/ExprNode.cc index 592b5d21..8b4bf616 100644 --- a/ExprNode.cc +++ b/ExprNode.cc @@ -164,6 +164,12 @@ NumConstNode::collectEndogenous(set > &result) const { } +void +NumConstNode::collectExogenous(set > &result) const +{ +} + + VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg) : ExprNode(datatree_arg), symb_id(symb_id_arg), @@ -473,6 +479,14 @@ VariableNode::collectEndogenous(set > &result) const result.insert(make_pair(symb_id, lag)); } +void +VariableNode::collectExogenous(set > &result) const +{ + if (type == eExogenous) + result.insert(make_pair(symb_id, lag)); +} + + UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) : ExprNode(datatree_arg), arg(arg_arg), @@ -864,6 +878,13 @@ UnaryOpNode::collectEndogenous(set > &result) const arg->collectEndogenous(result); } +void +UnaryOpNode::collectExogenous(set > &result) const +{ + arg->collectExogenous(result); +} + + BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, BinaryOpcode op_code_arg, const NodeID arg2_arg) : ExprNode(datatree_arg), @@ -1322,6 +1343,14 @@ BinaryOpNode::collectEndogenous(set > &result) const arg2->collectEndogenous(result); } + +void +BinaryOpNode::collectExogenous(set > &result) const +{ + arg1->collectExogenous(result); + arg2->collectExogenous(result); +} + TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) : ExprNode(datatree_arg), @@ -1334,7 +1363,7 @@ TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, // Non-null derivatives are the union of those of the arguments // Compute set union of arg{1,2,3}->non_null_derivatives - set non_null_derivatives_tmp; + set non_null_derivatives_tmp; set_union(arg1->non_null_derivatives.begin(), arg1->non_null_derivatives.end(), arg2->non_null_derivatives.begin(), @@ -1590,6 +1619,15 @@ TrinaryOpNode::collectEndogenous(set > &result) const arg3->collectEndogenous(result); } +void +TrinaryOpNode::collectExogenous(set > &result) const +{ + arg1->collectExogenous(result); + arg2->collectExogenous(result); + arg3->collectExogenous(result); +} + + UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg, int symb_id_arg, const vector &arguments_arg) : @@ -1650,6 +1688,15 @@ UnknownFunctionNode::collectEndogenous(set > &result) const (*it)->collectEndogenous(result); } +void +UnknownFunctionNode::collectExogenous(set > &result) const +{ + for(vector::const_iterator it = arguments.begin(); + it != arguments.end(); it++) + (*it)->collectExogenous(result); +} + + double UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (EvalException) { diff --git a/ModFile.cc b/ModFile.cc index 495b030a..f0d06425 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -69,7 +69,7 @@ ModFile::checkPass() exit(EXIT_FAILURE); } - if (mod_file_struct.simul_present && stochastic_statement_present) + if (mod_file_struct.simul_present && stochastic_statement_present && model_tree.mode==0) { cerr << "ERROR: A .mod file cannot contain both a simul command and one of {stoch_simul, estimation, forecast, osr, ramsey_policy}" << endl; exit(EXIT_FAILURE); diff --git a/ModelTree.cc b/ModelTree.cc index 59fda032..116c61f0 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -53,9 +53,10 @@ ModelTree::equation_number() const void ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, - const temporary_terms_type &temporary_terms) const + const temporary_terms_type &temporary_terms, + SymbolType type) const { - first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag))); + first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(type, symb_id, lag))); if (it != first_derivatives.end()) (it->second)->writeOutput(output, output_type, temporary_terms); else @@ -67,14 +68,9 @@ ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, { first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag))); if (it != first_derivatives.end()) - { - /*NodeID Id = it->second;*/ - (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx); - } + (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx); else - { - code_file.write(&FLDZ, sizeof(FLDZ)); - } + code_file.write(&FLDZ, sizeof(FLDZ)); } @@ -219,6 +215,35 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t } } + +void +ModelTree::BuildIncidenceMatrix() +{ + set > endogenous, exogenous; + for(int eq = 0; eq < (int) equations.size(); eq++) + { + BinaryOpNode *eq_node = equations[eq]; + endogenous.clear(); + NodeID Id = eq_node->arg1; + Id->collectEndogenous(endogenous); + Id = eq_node->arg2; + Id->collectEndogenous(endogenous); + for(set >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++) + { + block_triangular.incidencematrix.fill_IM(eq, it_endogenous->first, it_endogenous->second, eEndogenous); + } + exogenous.clear(); + Id = eq_node->arg1; + Id->collectExogenous(exogenous); + Id = eq_node->arg2; + Id->collectExogenous(exogenous); + for(set >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) + { + block_triangular.incidencematrix.fill_IM(eq, it_exogenous->first, it_exogenous->second, eExogenous); + } + } +} + void ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const { @@ -244,7 +269,7 @@ void ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) { map reference_count, first_occurence; - int i, j, m, eq, var, lag/*, prev_size=0*/; + int i, j, m, eq, var, lag; temporary_terms_type vect; ostringstream tmp_output; BinaryOpNode *eq_node; @@ -265,6 +290,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) tmp_output.str(""); lhs->writeOutput(tmp_output, oCDynamicModelSparseDLL, temporary_terms); tmp_s << "y[Per_y_+" << ModelBlock->Block_List[j].Variable[0] << "]"; + //Determine whether the equation could be evaluated rather than to be solved if (tmp_output.str()==tmp_s.str()) { if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE) @@ -285,6 +311,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) } } } + // Compute the temporary terms reordered for(i = 0;i < ModelBlock->Block_List[j].Size;i++) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; @@ -335,12 +362,11 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) for(second_derivatives_type::iterator it = second_derivatives.begin(); it != second_derivatives.end(); it++) it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, 0, ModelBlock, map_idx); - /*New*/ + // Add a mapping form node ID to temporary terms order j=0; for(temporary_terms_type::const_iterator it = temporary_terms.begin(); it != temporary_terms.end(); it++) map_idx[(*it)->idx]=j++; - /*EndNew*/ } void @@ -355,6 +381,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock ostringstream Uf[symbol_table.endo_nbr]; map reference_count; int prev_Simulation_Type=-1, count_derivates=0; + int jacobian_max_endo_col; temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); //---------------------------------------------------------------------- //Temporary variables declaration @@ -416,9 +443,9 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n"; else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) - output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, g1, g2, g3, y_index, jacobian_eval)\n"; + output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n"; else - output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods, g1, g2, g3)\n"; + output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, periods, jacobian_eval, g1, g2, g3, y_kmin, y_size)\n"; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl @@ -434,13 +461,13 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock temporary_terms_type tt2; - if(ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE) + if(ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) { int nze; for(nze=0,m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) nze+=ModelBlock->Block_List[j].IM_lead_lag[m].size; - //output << " Jacobian_Size=" << ModelBlock->Block_List[j].Size << "*(y_kmin+" << ModelBlock->Block_List[j].Max_Lead << " +periods);\n"; - //output << " g1=spalloc( y_size*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; + output << " g2=0;g3=0;\n"; + output << " b = [];\n"; output << " for it_ = y_kmin+1:(periods+y_kmin)\n"; output << " Per_y_=it_*y_size;\n"; output << " Per_J_=(it_-y_kmin-1)*y_size;\n"; @@ -467,7 +494,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock // 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); string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ; output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; @@ -502,7 +528,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock goto end; case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: - Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = residual(" << i+1 << ")"; + Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = residual(" << i+1 << ")"; output << sps << "residual(" << i+1 << ") = ("; goto end; case SOLVE_TWO_BOUNDARIES_COMPLETE: @@ -543,23 +569,42 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { if(ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]==ModelBlock->Block_List[j].Variable[0]) { - //output << " g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; - //output << " g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; - output << " g1(" << ModelBlock->Block_List[j].Equation[0]+1 << ", " << ModelBlock->Block_List[j].Variable[0]+1 + (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms); + output << " g1(" << ModelBlock->Block_List[j].Equation[0]+1 << ", " + << ModelBlock->Block_List[j].Variable[0]+1 + + (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr + << ")="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) - << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) << ") " << ModelBlock->Block_List[j].Variable[0]+1 << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; } } } + jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr; + for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1(" << eq+1 << ", " + << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr/*ModelBlock->Block_List[j].nb_exo*/ << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); + output << "; % variable=" << symbol_table.getNameByID(eExogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) { output << " else\n"; output << " g1="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms); + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) << ") " << ModelBlock->Block_List[j].Variable[0]+1 @@ -576,6 +621,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: count_derivates++; + output << " b = [];\n"; for(m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++) { k=m-ModelBlock->Block_List[j].Max_Lag; @@ -583,17 +629,35 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { 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 varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - output << " g1(" << eqr+1 << ", " << varr+1+m*ModelBlock->Block_List[j].Size << ")="; - writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; + //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + output << " g1(" << eq+1 << ", " << var+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) << "(" << /*variable_table.getLag(variable_table.getSymbolID(var))*/k << ") " << var+1 << ", equation=" << eq+1 << endl; } } - output << " else\n"; + jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr; + for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1(" << eq+1 << ", " << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); + output << "; % variable=" << symbol_table.getNameByID(eExogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + output << " else" << endl; + m=ModelBlock->Block_List[j].Max_Lag; for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { @@ -601,22 +665,20 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock 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]; - //Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")"; Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << ", " << varr+1 << ")*y(it_, " << var+1 << ")"; - //output << " u(" << u+1 << ") = "; output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; - writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms); + writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 << ", equation=" << eq+1 << endl; } - output << " end;\n"; for(i = 0;i < ModelBlock->Block_List[j].Size;i++) output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + output << " end;\n"; break; case SOLVE_TWO_BOUNDARIES_SIMPLE: case SOLVE_TWO_BOUNDARIES_COMPLETE: - output << " g2=0;g3=0;\n"; + output << " if ~jacobian_eval" << endl; for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) { k=m-ModelBlock->Block_List[j].Max_Lag; @@ -624,7 +686,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - //int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; if (k==0) @@ -637,14 +698,14 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; //output << " u(" << u+1 << "+Per_u_) = "; if(k==0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; else if(k==1) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; else if(k>0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; else if(k<0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; - writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; @@ -656,7 +717,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock } for(i = 0;i < ModelBlock->Block_List[j].Size;i++) { - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; #ifdef CONDITION output << " if (fabs(condition(" << i+1 << "))Block_List[j].Size;i++) output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; #endif + + output << " else" << endl; + for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;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]; + output << " g1(" << eqr+1 << ", " << varr+1+(m-ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lag_Endo)*ModelBlock->Block_List[j].Size << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size; + for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1(" << eqr+1 << ", " + << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); + output << "; % variable (exogenous)=" << symbol_table.getNameByID(eExogenous, var) + << "(" << k << ") " << var+1 << " " << varr+1 + << ", equation=" << eq+1 << endl; + } + } + output << " end;\n"; output << " end;\n"; break; + default: + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; } @@ -689,7 +789,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock void ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *ModelBlock, const string &static_basename) const { - int i,j,k,m, var, eq; + int i,j,k,m, var, eq, g1_index = 1; string tmp_s, sps; ostringstream tmp_output, global_output; NodeID lhs=NULL, rhs=NULL; @@ -734,9 +834,15 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode ||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 )) - skip_the_head=true; + { + skip_the_head=true; + g1_index++; + } else - skip_the_head=false; + { + skip_the_head=false; + g1_index = 1; + } if (!skip_the_head) { if (j>0) @@ -749,9 +855,9 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R ) - output << "function [y] = " << static_basename << "_" << j+1 << "(y, x)\n"; + output << "function [y, g1] = " << static_basename << "_" << j+1 << "(y, x, jacobian_eval)\n"; else - output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n"; + output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x, jacobian_eval)\n"; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl @@ -773,14 +879,17 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode memset(IM, 0, n*n*sizeof(bool)); for(m=-ModelBlock->Block_List[j].Max_Lag;m<=ModelBlock->Block_List[j].Max_Lead;m++) { - IMl=block_triangular.bGet_IM(m); - for(i=0;iBlock_List[j].Equation[i]; - for(k=0;kBlock_List[j].Variable[k]; - IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var]; + eq=ModelBlock->Block_List[j].Equation[i]; + for(k=0;kBlock_List[j].Variable[k]; + IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var]; + } } } } @@ -814,7 +923,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode // 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); + //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ; output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; @@ -862,108 +971,114 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode } } // The Jacobian if we have to solve the block - if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD_R) + output << " " << sps << "% Jacobian " << endl; + switch(ModelBlock->Block_List[j].Simulation_Type) { - output << " " << sps << "% Jacobian " << endl; - switch(ModelBlock->Block_List[j].Simulation_Type) + case EVALUATE_BACKWARD: + case EVALUATE_FORWARD: + case EVALUATE_BACKWARD_R: + case EVALUATE_FORWARD_R: + output << " if(jacobian_eval)\n"; + output << " g1( " << g1_index << ", " << g1_index << ")="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) + << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; + output << " end\n"; + break; + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: + output << " g1(1)="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) + << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; + break; + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + output << " g2=0;g3=0;\n"; + m=ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FORWARD_SIMPLE: - output << " g1(1)="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) - << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) - << ") " << ModelBlock->Block_List[j].Variable[0]+1 - << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; - break; - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_COMPLETE: - output << " g2=0;g3=0;\n"; - m=ModelBlock->Block_List[j].Max_Lag; + 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]; + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")"; + output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; + writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) + << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + break; + case SOLVE_TWO_BOUNDARIES_COMPLETE: + case SOLVE_TWO_BOUNDARIES_SIMPLE: + output << " g2=0;g3=0;\n"; + for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; for(i=0;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[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]; - //Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")"; - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")"; - //output << " u(" << u+1 << ") = "; - output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; - writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms); + if(!IM[eqr*ModelBlock->Block_List[j].Size+varr]) + { + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 + << ", " << varr+1 << ")*y( " << var+1 << ")"; + IM[eqr*ModelBlock->Block_List[j].Size+varr]=1; + } + output << " g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + "; + writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; - break; - case SOLVE_TWO_BOUNDARIES_COMPLETE: - output << " g2=0;g3=0;\n"; - for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for(i=0;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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - //output << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n"; - if(!IM[eqr*ModelBlock->Block_List[j].Size+varr]) - { - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 - << ", " << varr+1 << ")*y( " << var+1 << ")"; - IM[eqr*ModelBlock->Block_List[j].Size+varr]=1; - } - output << " g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + "; - writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; #ifdef CONDITION - output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++) - { - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; -#ifdef CONDITION - output << " if (fabs(condition(" << i+1 << "))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++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n"; - } - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; -#endif - break; } + output << " if(~jacobian_eval)\n"; + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; +#ifdef CONDITION + output << " if (fabs(condition(" << i+1 << "))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++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; + int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n"; + } + } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; +#endif + break; + default: + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; free(IM); } output << "return;\n\n\n"; - //free(IM); } @@ -1031,7 +1146,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; } ModelBlock_Aggregated_Count++; - //cout << "ModelBlock_Aggregated_Count=" << ModelBlock_Aggregated_Count << "\n"; //For each block j=0; for(k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++) @@ -1044,7 +1158,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl code_file.write(reinterpret_cast(&v),sizeof(v)); v=ModelBlock->Block_List[j].Simulation_Type; code_file.write(reinterpret_cast(&v),sizeof(v)); - //cout << "FBEGINBLOCK j=" << j << " size=" << ModelBlock_Aggregated_Number[k0] << " type=" << v << "\n"; for(k=0; kBlock_List[j].Size;i++) @@ -1079,7 +1192,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl } for(k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++) { - //For a block composed of a single equation determines wether we have to evaluate or to solve the equation + //For a block composed of a single equation determines whether we have to evaluate or to solve the equation if (ModelBlock->Block_List[j].Size==1) { lhs_rhs_done=true; @@ -1089,10 +1202,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl } else lhs_rhs_done=false; - /*if (ModelBlock->Block_List[j].Size==1) - lhs_rhs_done=true; - else - lhs_rhs_done=false;*/ //The Temporary terms temporary_terms_type tt2; i=0; @@ -1124,7 +1233,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl // 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); + //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); if (!lhs_rhs_done) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; @@ -1343,6 +1452,8 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n"; #endif break; + default: + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; @@ -1734,13 +1845,10 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size; int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j]; int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j]; - /*cout << " ! IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr+k1*block_triangular.ModelBlock->Block_List[num].Size << "), " << k1 << ")] = " << u << ";\n"; - cout << " ? IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr << "), " << k1 << ")] = " << u << ";\n";*/ SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); SaveCode.write(reinterpret_cast(&u), sizeof(u)); - //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n"; u_count_int++; } } @@ -1748,13 +1856,12 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b { int eqr1=j; int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods - +block_triangular.Model_Max_Lead); + +block_triangular.incidencematrix.Model_Max_Lead_Endo); int k1=0; SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); - //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n"; u_count_int++; } for(j=0;jBlock_List[num].Size;j++) @@ -1775,7 +1882,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b { string filename; ofstream mStaticModelFile; - int i, k, prev_Simulation_Type; + int i, k, prev_Simulation_Type, ga_index = 1; bool /*printed = false,*/ skip_head, open_par=false; filename = static_basename + ".m"; mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); @@ -1812,6 +1919,12 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; } mStaticModelFile << " ];\n"; + mStaticModelFile << " y_index_eq=["; + for(int ik=0;ikBlock_List[i].Size;ik++) + { + mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; + } + mStaticModelFile << " ];\n"; k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) @@ -1825,21 +1938,27 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b case EVALUATE_FORWARD_R: case EVALUATE_BACKWARD_R: if(!skip_head) - mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n"; + { + ga_index = 1; + mStaticModelFile << " [y, ga]=" << static_basename << "_" << i + 1 << "(y, x, 1);\n"; + } + else + ga_index++; mStaticModelFile << " residual(y_index)=ys(y_index)-y(y_index);\n"; + mStaticModelFile << " g1(y_index_eq, y_index) = ga(" << ga_index << ", " << ga_index << ");\n"; break; case SOLVE_FORWARD_COMPLETE: case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_SIMPLE: case SOLVE_BACKWARD_SIMPLE: case SOLVE_TWO_BOUNDARIES_COMPLETE: - mStaticModelFile << " [r, g1]=" << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " [r, g1(y_index_eq, y_index)]=" << static_basename << "_" << i + 1 << "(y, x, 1);\n"; mStaticModelFile << " residual(y_index)=r;\n"; break; } prev_Simulation_Type=k; } - mStaticModelFile << " varargout{1}=residual;\n"; + mStaticModelFile << " varargout{1}=residual';\n"; mStaticModelFile << " varargout{2}=g1;\n"; mStaticModelFile << " return;\n"; mStaticModelFile << " end;\n"; @@ -1867,7 +1986,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b { mStaticModelFile << " end\n"; } - mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x, 0);\n"; } open_par=false; } @@ -1880,32 +1999,16 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b open_par=false; mStaticModelFile << " g1=0;\n"; mStaticModelFile << " r=0;\n"; - /*mStaticModelFile << " for it_=y_kmin+1:periods+y_kmin\n"; - mStaticModelFile << " cvg=0;\n"; - mStaticModelFile << " iter=0;\n"; - mStaticModelFile << " Per_y_=it_*y_size;\n"; - mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; - mStaticModelFile << " y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r/g1;\n"; - mStaticModelFile << " cvg=((r[0]*r[0])maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n"; mStaticModelFile << " y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mStaticModelFile << " cvg=((r*r)maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; - mStaticModelFile << " [L, U] = LU(g1);\n"; - mStaticModelFile << " y(it_, y_index) = U\\(L\\b);\n"; - mStaticModelFile << " cvg=((r'*r)maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n"; mStaticModelFile << " max_res=max(abs(r));\n"; mStaticModelFile << " cvg=(max_resBlock_List[i].Size;ik++) { tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; } - //mDynamicModelFile << " ];\n"; if(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R) { if(i==block_triangular.ModelBlock->Size-1) @@ -2127,8 +2198,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string case SOLVE_FORWARD_SIMPLE: case SOLVE_BACKWARD_SIMPLE: mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n"; - mDynamicModelFile << " y_index = " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n"; - mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, y_index, 1);\n"; + mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n"; mDynamicModelFile << " residual(y_index_eq)=r;\n"; tmp_eq.str(""); tmp.str(""); @@ -2136,41 +2206,38 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string case SOLVE_FORWARD_COMPLETE: case SOLVE_BACKWARD_COMPLETE: mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; - tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1; - mDynamicModelFile << " y_index_c = y_index;\n"; - mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n"; - mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n"; - mDynamicModelFile << " end;\n"; - //tmp_i=variable_table.max_lag+variable_table.max_lead+1; - mDynamicModelFile << " ga = [];\n"; - mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n"; - //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, ga, g2, g3);\n"; - mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, ga, g2, g3);\n"; - mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n"; + mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n"; mDynamicModelFile << " residual(y_index_eq)=r;\n"; break; case SOLVE_TWO_BOUNDARIES_COMPLETE: case SOLVE_TWO_BOUNDARIES_SIMPLE: - //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, 1);\n"; + int j; mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; + mDynamicModelFile << " y_index = ["; + for(j=0;jBlock_List[i].Size;ik++) + { + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr; + } + int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1; + for(j=0;jBlock_List[i].nb_exo;ik++) + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr+symbol_table.endo_nbr*tmp_i; + mDynamicModelFile << " ];\n"; tmp.str(""); tmp_eq.str(""); - tmp_i=variable_table.max_lag+variable_table.max_lead+1; mDynamicModelFile << " ga = [];\n"; - mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n"; - mDynamicModelFile << " y_index_c = y_index;\n"; - tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1; - mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n"; - mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n"; - mDynamicModelFile << " end;\n"; - //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n"; - mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n"; + j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1) + + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1); + mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " << + block_triangular.ModelBlock->Block_List[i].Size*j << ");\n"; + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; + mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-" << variable_table.max_lag << ", 1, ga, g2, g3, " << variable_table.max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n"; if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead) - mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n"; + mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n"; else - mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n"; + mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n"; mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n"; break; } @@ -2256,7 +2323,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " iter=0;\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n"; + mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n"; mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mDynamicModelFile << " cvg=((r*r)maxit_),\n"; - mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n"; + mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n"; mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mDynamicModelFile << " cvg=((r*r)Block_List[i].Size)) { - /*if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - if (!printed) - { - printed = true; - } - SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, symbol_table.endo_nbr); - Nb_SGE++; - mDynamicModelFile << " Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";*/ if (open_par) mDynamicModelFile << " end\n"; open_par=false; @@ -2348,17 +2405,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " return;\n"; mDynamicModelFile << " end\n"; mDynamicModelFile << " end\n"; + /*cerr << "Not implemented block SOLVE_FORWARD_COMPLETE" << endl; exit(EXIT_FAILURE);*/ } else if ((k == SOLVE_BACKWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size)) { - /*if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, symbol_table.endo_nbr); - Nb_SGE++; - mDynamicModelFile << " Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";*/ if (open_par) mDynamicModelFile << " end\n"; open_par=false; @@ -2398,7 +2450,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string printed = true; } Nb_SGE++; - //cout << "new_SGE=" << new_SGE << "\n"; mDynamicModelFile << " cvg=0;\n"; mDynamicModelFile << " iter=0;\n"; @@ -2410,13 +2461,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string } mDynamicModelFile << " ];\n"; mDynamicModelFile << " Blck_size=" << block_triangular.ModelBlock->Block_List[i].Size << ";\n"; - /*mDynamicModelFile << " if(options_.simulation_method==2 | options_.simulation_method==3),\n"; - mDynamicModelFile << " [r, g1]= " << basename << "_static(y, x);\n"; - mDynamicModelFile << " [L1,U1] = lu(g1,1e-5);\n"; - mDynamicModelFile << " I = speye(periods);\n"; - mDynamicModelFile << " L1=kron(I,L1);\n"; - mDynamicModelFile << " U1=kron(I,U1);\n"; - mDynamicModelFile << " end;\n";*/ mDynamicModelFile << " y_kmin_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lag << ";\n"; mDynamicModelFile << " y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n"; mDynamicModelFile << " lambda=options_.slowc;\n"; @@ -2428,12 +2472,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; mDynamicModelFile << " Jacobian_Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*(y_kmin+" << block_triangular.ModelBlock->Block_List[i].Max_Lead << " +periods);\n"; mDynamicModelFile << " g1=spalloc( length(y_index)*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; - /*mDynamicModelFile << " cpath=path;\n"; - mDynamicModelFile << " addpath(fullfile(matlabroot,'toolbox','matlab','sparfun'));\n";*/ - mDynamicModelFile << " bicgstabh=@bicgstab;\n"; - //mDynamicModelFile << " path(cpath);\n"; mDynamicModelFile << sp << " reduced = 0;\n"; - //mDynamicModelFile << " functions(bicgstabh)\n"; if (!block_triangular.ModelBlock->Block_List[i].is_linear) { sp=" "; @@ -2443,7 +2482,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string { sp=""; } - mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\n"; + mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, periods, 0, g1, g2, g3, y_kmin, Blck_size);\n"; mDynamicModelFile << sp << " g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);\n"; mDynamicModelFile << sp << " b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';\n"; mDynamicModelFile << sp << " if(~isreal(r))\n"; @@ -2525,7 +2564,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << sp << " flag1=1;\n"; mDynamicModelFile << sp << " while(flag1>0)\n"; mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n"; - mDynamicModelFile << sp << " [za,flag1] = bicgstabh(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; + mDynamicModelFile << sp << " [za,flag1] = bicgstab(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; mDynamicModelFile << sp << " if (flag1>0 | reduced)\n"; mDynamicModelFile << sp << " if(flag1==1)\n"; mDynamicModelFile << sp << " disp(['No convergence inside BICGSTAB after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n"; @@ -2791,11 +2830,12 @@ ModelTree::writeOutput(ostream &output) const */ output << "M_.lead_lag_incidence = ["; // Loop on endogenous variables + int lag = 0; for(int endoID = 0; endoID < symbol_table.endo_nbr; endoID++) { output << "\n\t"; // Loop on periods - for(int lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++) + for(lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++) { // Print variableID if exists with current period, otherwise print 0 try @@ -2819,7 +2859,7 @@ ModelTree::writeOutput(ostream &output) const bool skip_the_head; int k=0; int count_lead_lag_incidence = 0; - int max_lead, max_lag; + int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo; for(int j = 0;j < block_triangular.ModelBlock->Size;j++) { //For a block composed of a single equation determines wether we have to evaluate or to solve the equation @@ -2835,9 +2875,17 @@ ModelTree::writeOutput(ostream &output) const k++; count_lead_lag_incidence = 0; int Block_size=block_triangular.ModelBlock->Block_List[j].Size; + //int Block_exo_nbr=block_triangular.ModelBlock->Block_List[j].exo_nbr; max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ; max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead; + max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ; + max_lead_endo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Endo; + max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ; + max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo; bool evaluate=false; + vector exogenous;//, tmp_exogenous(block_triangular.exo_nbr,-1); + vector::iterator it_exogenous; + exogenous.clear(); ostringstream tmp_s, tmp_s_eq; tmp_s.str(""); tmp_s_eq.str(""); @@ -2846,6 +2894,13 @@ ModelTree::writeOutput(ostream &output) const tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1; tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1; } + for(int i=0;iBlock_List[j].nb_exo;i++) + { + int ii=block_triangular.ModelBlock->Block_List[j].Exogenous[i]; + for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) /*cout << "*it_exogenous=" << *it_exogenous << "\n"*/; + if(it_exogenous==exogenous.end() || exogenous.begin()==exogenous.end()) + exogenous.push_back(ii); + } if (block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R @@ -2864,78 +2919,74 @@ ModelTree::writeOutput(ostream &output) const max_lag =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag ; if(max_leadBlock_List[j+Block_size].Max_Lead) max_lead=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead; - //cout << "block_triangular.ModelBlock->Block_List[" << j+Block_size << "].Size=" << block_triangular.ModelBlock->Block_List[j+Block_size].Size << "\n"; + if(max_lag_endo Block_List[j+Block_size].Max_Lag_Endo ) + max_lag_endo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Endo ; + if(max_lead_endoBlock_List[j+Block_size].Max_Lead_Endo) + max_lead_endo=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Endo; + if(max_lag_exo Block_List[j+Block_size].Max_Lag_Exo ) + max_lag_exo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Exo ; + if(max_lead_exoBlock_List[j+Block_size].Max_Lead_Exo) + max_lead_exo=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Exo; for(int i=0;iBlock_List[j+Block_size].Size;i++) { tmp_s << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Variable[i]+1; tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Equation[i]+1; } + for(int i=0;iBlock_List[j+Block_size].nb_exo;i++) + { + int ii=block_triangular.ModelBlock->Block_List[j+Block_size].Exogenous[i]; + //for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) cout << "*it_exogenous=" << *it_exogenous << "\n"; + if(it_exogenous==exogenous.end()) + exogenous.push_back(ii); + } Block_size+=block_triangular.ModelBlock->Block_List[j+Block_size].Size; } - //cout << "i=" << i << " max_lag=" << max_lag << " max_lead=" << max_lead << "\n"; } } output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n"; - //output << "M_.block_structure.block(" << k << ").size = " << block_triangular.ModelBlock->Block_List[j].Size << ";\n"; output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_lead = " << max_lead << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag_endo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead_endo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_exo_lag = " << max_lag_exo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_exo_lead = " << max_lead_exo << ";\n"; output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n"; output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n"; output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n"; + output << "M_.block_structure.block(" << k << ").exogenous = ["; + int i=0; + for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) + if(*it_exogenous>=0) + { + output << " " << *it_exogenous+1; + i++; + } + output << "];\n"; + output << "M_.block_structure.block(" << k << ").exo_nbr = " << i << ";\n"; + + output << "M_.block_structure.block(" << k << ").exo_det_nbr = " << block_triangular.ModelBlock->Block_List[j].nb_exo_det << ";\n"; + tmp_s.str(""); - cout << "begining of lead_lag_incidence\n"; bool done_IM=false; if(!evaluate) { output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n"; - for(int l=-max_lag;lBlock_List[j].Size;l_var++) + tmp_IM=block_triangular.incidencematrix.bGet_IM(l, eEndogenous); + if(tmp_IM) { - for(int l_equ=0;l_equBlock_List[j].Size;l_equ++) - if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) - { - count_lead_lag_incidence++; - if(tmp_s.str().length()) - tmp_s << " "; - tmp_s << count_lead_lag_incidence; - done_IM=true; - break; - } - if(!done_IM) - tmp_s << " 0"; - done_IM=false; - } - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; - tmp_s.str(""); - } - } - else - { - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; - for(int l=-max_lag;lBlock_List[ii].Size;l_var++) + for(int l_var=0;l_varBlock_List[j].Size;l_var++) { - for(int l_equ=0;l_equBlock_List[ii].Size;l_equ++) - if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) + for(int l_equ=0;l_equBlock_List[j].Size;l_equ++) + if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) { - //if(not_increm && l==-max_lag) - count_lead_lag_incidence++; - not_increm=false; + count_lead_lag_incidence++; if(tmp_s.str().length()) tmp_s << " "; - //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; tmp_s << count_lead_lag_incidence; done_IM=true; break; @@ -2944,10 +2995,53 @@ ModelTree::writeOutput(ostream &output) const tmp_s << " 0"; done_IM=false; } - ii++; + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; + tmp_s.str(""); + } + } + } + else + { + bool done_some_where; + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; + for(int l=-max_lag_endo;lBlock_List[ii].Size;l_var++) + { + for(int l_equ=0;l_equBlock_List[ii].Size;l_equ++) + if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) + { + //if(not_increm && l==-max_lag) + count_lead_lag_incidence++; + not_increm=false; + if(tmp_s.str().length()) + tmp_s << " "; + //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; + tmp_s << count_lead_lag_incidence; + done_IM=true; + break; + } + if(!done_IM) + tmp_s << " 0"; + else + done_some_where = true; + done_IM=false; + } + ii++; + } + //if(done_some_where) + output << tmp_s.str() << "\n"; + tmp_s.str(""); } - output << tmp_s.str() << "\n"; - tmp_s.str(""); } output << "];\n"; } @@ -2955,13 +3049,13 @@ ModelTree::writeOutput(ostream &output) const prev_Simulation_Type=block_triangular.ModelBlock->Block_List[j].Simulation_Type; } - for(int j=-block_triangular.Model_Max_Lag;jfirst.second) == eEndogenous) { NodeID Id = it->second; - double val; + double val = 0; try { val = Id->eval(eval_context); @@ -3046,7 +3140,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_ int k1=variable_table.getLag(it->first.second); if (a_variable_lag!=k1) { - IM=block_triangular.bGet_IM(k1); + IM=block_triangular.incidencematrix.bGet_IM(k1, eEndogenous); a_variable_lag=k1; } if (k1==0) @@ -3058,7 +3152,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_ { if(block_triangular.bt_verbose) cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr << ")\n"; - block_triangular.unfill_IM(eq, var, k1); + block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous); i++; } } @@ -3159,6 +3253,8 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term if (mode == eSparseDLLMode || mode == eSparseMode) { + block_triangular.incidencematrix.init_incidence_matrix(); + BuildIncidenceMatrix(); jacob_map j_m; evaluateJacobian(eval_context, &j_m); @@ -3166,7 +3262,7 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term if (block_triangular.bt_verbose) { cout << "The gross incidence matrix \n"; - block_triangular.Print_IM( symbol_table.endo_nbr); + block_triangular.incidencematrix.Print_IM(eEndogenous); } block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m); BlockLinear(block_triangular.ModelBlock); @@ -3208,7 +3304,8 @@ ModelTree::writeDynamicFile(const string &basename) const case eSparseMode: writeSparseDynamicMFile(basename + "_dynamic", basename, mode); block_triangular.Free_Block(block_triangular.ModelBlock); - block_triangular.Free_IM(block_triangular.First_IM); + block_triangular.incidencematrix.Free_IM(); + //block_triangular.Free_IM_X(block_triangular.First_IM_X); break; case eDLLMode: writeDynamicCFile(basename + "_dynamic"); @@ -3216,7 +3313,8 @@ ModelTree::writeDynamicFile(const string &basename) const case eSparseDLLMode: writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL); block_triangular.Free_Block(block_triangular.ModelBlock); - block_triangular.Free_IM(block_triangular.First_IM); + block_triangular.incidencematrix.Free_IM(); + //block_triangular.Free_IM_X(block_triangular.First_IM_X); break; } } diff --git a/ParsingDriver.cc b/ParsingDriver.cc index 1516a026..c1109e1e 100644 --- a/ParsingDriver.cc +++ b/ParsingDriver.cc @@ -169,7 +169,7 @@ ParsingDriver::add_model_variable(string *name, string *olag) if (type == eUnknownFunction) error("Symbol " + *name + " is a function name unknown to Dynare. It cannot be used inside model."); - if (type == eExogenous && lag != 0) + if (type == eExogenous && lag != 0 && (model_tree->mode != eSparseDLLMode && model_tree->mode != eSparseMode)) warning("Exogenous variable " + *name + " has lead/lag " + *olag); if (type == eModelLocalVariable && lag != 0) @@ -179,9 +179,13 @@ ParsingDriver::add_model_variable(string *name, string *olag) NodeID id = model_tree->AddVariable(*name, lag); - if ((type == eEndogenous) && (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode)) - model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); - + /*if (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode) + { + if (type == eEndogenous) + model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); + if (type == eExogenous) + model_tree->block_triangular.fill_IM_X(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); + }*/ delete name; delete olag; return id; @@ -360,14 +364,16 @@ void ParsingDriver::sparse_dll() { model_tree->mode = eSparseDLLMode; - model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/ } void ParsingDriver::sparse() { model_tree->mode = eSparseMode; - model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/ } void diff --git a/include/BlockTriangular.hh b/include/BlockTriangular.hh index 894b8ccb..fd9e7898 100644 --- a/include/BlockTriangular.hh +++ b/include/BlockTriangular.hh @@ -36,41 +36,61 @@ struct List_IM bool* IM; }; + +//! create and manage the incidence matrix +class IncidenceMatrix //: SymbolTable +{ + //friend class BlockTriangular; +public: +const SymbolTable &symbol_table; + IncidenceMatrix(const SymbolTable &symbol_table_arg); + List_IM* Build_IM(int lead_lag, SymbolType type); + List_IM* Get_IM(int lead_lag, SymbolType type) const; + bool* bGet_IM(int lead_lag, SymbolType type) const; + void fill_IM(int equation, int variable_endo, int lead_lag, SymbolType type); + void unfill_IM(int equation, int variable_endo, int lead_lag, SymbolType type); + void init_incidence_matrix(); + void Free_IM() const; + List_IM* Get_First(SymbolType type) const; + void Print_IM(SymbolType type) const; + void Print_SIM(bool* IM, SymbolType type) const; + + void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const; +private: + List_IM *First_IM, *Last_IM, *First_IM_X, *Last_IM_X ; +public: + int Model_Max_Lead, Model_Max_Lag; + int Model_Max_Lead_Endo, Model_Max_Lag_Endo, Model_Max_Lead_Exo, Model_Max_Lag_Exo; +}; + + //! Matrix of doubles for representing jacobian typedef map,double> jacob_map; //! Create the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition class BlockTriangular { + //friend class IncidenceMatrix; public: - BlockTriangular(const SymbolTable &symbol_table_arg); const SymbolTable &symbol_table; + BlockTriangular(const SymbolTable &symbol_table_arg); + //BlockTriangular(const IncidenceMatrix &incidence_matrix_arg); + //const SymbolTable &symbol_table; Blocks blocks; Normalization normalization; - List_IM* Build_IM(int lead_lag); - List_IM* Get_IM(int lead_lag); - bool* bGet_IM(int lead_lag) const; - void fill_IM(int equation, int variable_endo, int lead_lag); - void unfill_IM(int equation, int variable_endo, int lead_lag); - void init_incidence_matrix(int nb_endo); - void Print_IM(int n) const; - void Free_IM(List_IM* First_IM) const; - void Print_SIM(bool* IM, int n) const; + IncidenceMatrix incidencematrix; void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m); bool Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0 , jacob_map j_m); void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0); - void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n); void Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock); void Free_Block(Model_Block* ModelBlock) const; - List_IM *First_IM ; - List_IM *Last_IM ; simple *Index_Equ_IM; simple *Index_Var_IM; int prologue, epilogue; - int Model_Max_Lead, Model_Max_Lag, periods; bool bt_verbose; - int endo_nbr; + //int endo_nbr, exo_nbr; Model_Block* ModelBlock; + int periods; inline static std::string BlockType0(int type) { switch (type) @@ -98,14 +118,14 @@ public: { case EVALUATE_FORWARD: case EVALUATE_FORWARD_R: - return ("EVALUATE FORWARD "); + return ("EVALUATE FORWARD "); break; case EVALUATE_BACKWARD: case EVALUATE_BACKWARD_R: return ("EVALUATE BACKWARD "); break; case SOLVE_FORWARD_SIMPLE: - return ("SOLVE FORWARD SIMPLE "); + return ("SOLVE FORWARD SIMPLE "); break; case SOLVE_BACKWARD_SIMPLE: return ("SOLVE BACKWARD SIMPLE "); @@ -114,7 +134,7 @@ public: return ("SOLVE TWO BOUNDARIES SIMPLE "); break; case SOLVE_FORWARD_COMPLETE: - return ("SOLVE FORWARD COMPLETE "); + return ("SOLVE FORWARD COMPLETE "); break; case SOLVE_BACKWARD_COMPLETE: return ("SOLVE BACKWARD COMPLETE "); diff --git a/include/ExprNode.hh b/include/ExprNode.hh index 55a2f9e1..d3d24e35 100644 --- a/include/ExprNode.hh +++ b/include/ExprNode.hh @@ -143,6 +143,7 @@ public: /*! Endogenous are stored as integer pairs of the form (symb_id, lag) They are added to the set given in argument */ virtual void collectEndogenous(set > &result) const = 0; + virtual void collectExogenous(set > &result) const = 0; virtual void computeTemporaryTerms(map &reference_count, temporary_terms_type &temporary_terms, map &first_occurence, @@ -179,6 +180,7 @@ public: NumConstNode(DataTree &datatree_arg, int id_arg); virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; }; @@ -198,6 +200,7 @@ public: VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg); virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms = temporary_terms_type()) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; }; @@ -223,6 +226,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -250,6 +254,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -282,6 +287,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -307,6 +313,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; }; @@ -314,21 +321,22 @@ public: //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model struct IM_compact { - int size, u_init, u_finish, nb_endo; - int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index; + int size, u_init, u_finish, nb_endo, size_exo; + int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index; }; //! One block of the model struct Block { - int Size, Sized; + int Size, Sized, nb_exo, nb_exo_det; BlockType Type; BlockSimulationType Simulation_Type; int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo; + int Max_Lag_Endo, Max_Lead_Endo; + int Max_Lag_Exo, Max_Lead_Exo; bool is_linear; int *Equation, *Own_Derivative; - int *Variable, *Variable_Sorted, *dVariable; - int *variable_dyn_index, *variable_dyn_leadlag; + int *Variable, *Exogenous; temporary_terms_type *Temporary_terms; IM_compact *IM_lead_lag; int Code_Start, Code_Length; @@ -339,7 +347,7 @@ struct Model_Block { int Size, Periods; Block* Block_List; - int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block; + //int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block; }; #endif diff --git a/include/ModelTree.hh b/include/ModelTree.hh index 34645599..058bb917 100644 --- a/include/ModelTree.hh +++ b/include/ModelTree.hh @@ -26,6 +26,7 @@ using namespace std; #include #include #include +#include #include "SymbolTable.hh" #include "NumericalConstants.hh" @@ -82,12 +83,14 @@ private: //! Computes derivatives of ModelTree void derive(int order); //! Write derivative of an equation w.r. to a variable - void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; + void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, SymbolType type) const; //! Write derivative code of an equation w.r. to a variable void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const; //! Computes temporary terms void computeTemporaryTerms(int order); void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock); + //! Build The incidence matrix form the modeltree + void BuildIncidenceMatrix(); //! Writes temporary terms void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const; //! Writes model local variables