diff --git a/BlockTriangular.cc b/BlockTriangular.cc
index bce845494b1a80b95131462eb0405da26b1b30fb..80631141ac679fed169715ecfc8676244a500b5d 100644
--- a/BlockTriangular.cc
+++ b/BlockTriangular.cc
@@ -88,17 +88,22 @@ BlockTriangular::init_incidence_matrix(int nb_endo)
void
-BlockTriangular::Free_IM(List_IM* First_IM)
+BlockTriangular::Free_IM(List_IM* First_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;
free(Cur_IM->IM);
+ delete Cur_IM;
Cur_IM = First_IM;
}
- free(SFirst_IM);
+ //free(SFirst_IM);
+ //delete SFirst_IM;
}
//------------------------------------------------------------------------------
@@ -223,7 +228,7 @@ BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Inde
SIM[pos2*n + j] = tmp_b;
}
}
- /* ...and variables (colomn)*/
+ /* ...and variables (column)*/
if(pos1 != pos3)
{
tmp_i = Index_Var_IM[pos1].index;
@@ -241,7 +246,7 @@ BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Inde
//------------------------------------------------------------------------------
// Find the prologue and the epilogue of the model
void
-BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM)
+BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0)
{
bool modifie = 1;
int i, j, k, l = 0;
@@ -261,7 +266,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
l = j;
}
}
- if ((k == 1) /* && (l==i)*/)
+ 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);
@@ -273,6 +278,8 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
}
*epilogue = 0;
modifie = 1;
+ /*print_SIM(IM,n);
+ print_SIM(IM*/
while(modifie)
{
modifie = 0;
@@ -287,7 +294,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
l = j;
}
}
- if ((k == 1) /* && (l==i)*/)
+ 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);
@@ -330,10 +337,10 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int
while(Cur_IM)
{
k = Cur_IM->lead_lag;
- i_1 = Index_Var_IM[*count_Equ].index * endo_nbr;
+ i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr;
if(k > 0)
{
- if(Cur_IM->IM[i_1 + Index_Equ_IM[ /*j*/*count_Equ].index])
+ if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index])
{
nb_lead_lag_endo++;
tmp_size[Model_Max_Lag + k]++;
@@ -348,7 +355,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int
else
{
k = -k;
- if(Cur_IM->IM[i_1 + Index_Equ_IM[ /*j*/*count_Equ].index])
+ if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index])
{
tmp_size[Model_Max_Lag - k]++;
nb_lead_lag_endo++;
@@ -374,39 +381,45 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int
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;
Index_Equ_IM[*count_Equ].block = *count_Block;
- cout << "Lead=" << Lead << " Lag=" << Lag << "\n";
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_FOREWARD_SIMPLE;
+ 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].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
- cout << "alloc ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag = (" << (Lead + Lag + 1) * sizeof(IM_compact) << ")\n";
- 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;
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));
- ls = l = 1;
- i1 = 0;
- for(i = 0;i < Lead + Lag + 1;i++)
+ }
+ 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])
{
- 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_size[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);
+ 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));
+ ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l;
+ IM = bGet_IM(li - Lag);
if(IM == NULL)
{
- cout << "Error IM(" << i - Lag << ") doesn't exist\n";
+ cout << "Error IM(" << li - Lag << ") doesn't exist\n";
exit( -1);
}
- if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*endo_nbr])
+ 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] = i - Lag;
+ ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = li - Lag;
tmp_var[0] = i1;
i1++;
}
@@ -414,31 +427,28 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int
i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr;
if(IM[Index_Var_IM[*count_Equ].index + i_1])
{
- if(i == Lag)
+ if(li == Lag)
{
- ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls;
+ ModelBlock->Block_List[*count_Block].IM_lead_lag[li].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] = 0;
- ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = 0;
- ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[*count_Equ].index;
- ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[*count_Equ].index;
- ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[0]];
+ 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++;
}
- ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1;
- }
+ ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1;
+ }
}
- 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_FOREWARD_SIMPLE;
(*count_Equ)++;
(*count_Block)++;
free(tmp_size);
free(tmp_endo);
+ free(tmp_var);
}
}
else
@@ -535,7 +545,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int
}
ModelBlock->Block_List[*count_Block].Max_Lag = Lag;
ModelBlock->Block_List[*count_Block].Max_Lead = Lead;
- cout << "alloc ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag = (" << (Lead + Lag + 1) * sizeof(IM_compact) << ")\n";
ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact));
ls = l = size;
i1 = 0;
@@ -621,7 +630,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int
void
-BlockTriangular::Free_Block(Model_Block* ModelBlock)
+BlockTriangular::Free_Block(Model_Block* ModelBlock) const
{
int blk, i;
#ifdef DEBUG
@@ -637,12 +646,33 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock)
free(ModelBlock->Block_List[blk].Equation);
free(ModelBlock->Block_List[blk].Variable);
free(ModelBlock->Block_List[blk].Variable_Sorted);
+ 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);
+ }
+ free(ModelBlock->Block_List[blk].IM_lead_lag);
+ delete(ModelBlock->Block_List[blk].Temporary_terms);
}
else
{
free(ModelBlock->Block_List[blk].Equation);
free(ModelBlock->Block_List[blk].Variable);
free(ModelBlock->Block_List[blk].Variable_Sorted);
+ 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);
@@ -651,8 +681,10 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock)
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);
}
free(ModelBlock->Block_List[blk].IM_lead_lag);
+ delete(ModelBlock->Block_List[blk].Temporary_terms);
}
}
free(ModelBlock->in_Block_Equ);
@@ -661,6 +693,8 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock)
free(ModelBlock->in_Var_of_Block);
free(ModelBlock->Block_List);
free(ModelBlock);
+ free(Index_Equ_IM);
+ free(Index_Var_IM);
}
//------------------------------------------------------------------------------
@@ -672,10 +706,10 @@ 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;
+ //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_First_IM = (List_IM*)malloc(sizeof(*p_First_IM));
p_Cur_IM = p_First_IM;
Cur_IM = First_IM;
i = endo_nbr * endo_nbr;
@@ -692,82 +726,110 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
}
else
p_Cur_IM->pNext = NULL;
- }
- Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM);
+ }*/
+ 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);
+ cout << "IM\n";
Print_SIM(IM, n);
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";
}
- if(Do_Normalization)
+ if(*prologue+*epilogue<n)
{
- cout << "Normalizing the model ...\n";
- if(mixing)
+ if(Do_Normalization)
{
- double* max_val=(double*)malloc(n*sizeof(double));
- memset(max_val,0,n*sizeof(double));
- for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
- {
- if(fabs(iter->second)>max_val[iter->first.first])
- max_val[iter->first.first]=fabs(iter->second);
- }
- for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
- iter->second/=max_val[iter->first.first];
- free(max_val);
- bool OK=false;
- double bi=0.99999999;
- int suppressed=0;
- while(!OK && bi>1e-14)
+ cout << "Normalizing the model ...\n";
+ if(mixing)
{
- int suppress=0;
- SIM0 = (bool*)malloc(n * n * sizeof(bool));
- memset(SIM0,0,n*n*sizeof(bool));
- SIM00 = (bool*)malloc(n * n * sizeof(bool));
- memset(SIM00,0,n*n*sizeof(bool));
+ 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++ )
- {
- if(fabs(iter->second)>bi)
- {
- SIM0[iter->first.first*n+iter->first.second]=1;
- if(!IM_0[iter->first.first*n+iter->first.second])
+ {
+ //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);
+ }
+ for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
+ iter->second/=max_val[iter->first.first];
+ free(max_val);
+ bool OK=false;
+ double bi=0.99999999;
+ int suppressed=0;
+ while(!OK && bi>1e-14)
+ {
+ 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])
+ {
+ cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n";
+ }
+ }
+ else
+ suppress++;
+ }
+ //cout << "n*n=" << n*n << "\n";
+ for(i = 0;i < n;i++)
+ for(j = 0;j < n;j++)
{
- cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n";
+ //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];
}
- }
- else
- suppress++;
- }
- for(i = 0;i < n;i++)
- for(j = 0;j < n;j++)
- SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index];
- free(SIM0);
- if(suppress!=suppressed)
- OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
- suppressed=suppress;
+ //cout << "free SIM0=" << SIM0 << "\n";
+ free(SIM0);
+ if(suppress!=suppressed)
+ {
+ OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
+ if(!OK)
+ normalization.Free_Equation(n-*prologue-*epilogue,Equation_gr);
+ }
+ suppressed=suppress;
+ if(!OK)
+ bi/=1.07;
+ if(bi>1e-14)
+ free(SIM00);
+ }
if(!OK)
- bi/=1.07;
- if(bi>1e-14)
- free(SIM00);
- }
- if(!OK)
- {
- normalization.Set_fp_verbose(true);
- OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
- cout << "Error\n";
- exit(-1);
+ {
+ normalization.Set_fp_verbose(true);
+ OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
+ cout << "Error\n";
+ exit(-1);
+ }
}
+ else
+ normalization.Normalize(n, *prologue, *epilogue, IM, Index_Equ_IM, Equation_gr, 0, 0);
}
else
- normalization.Normalize(n, *prologue, *epilogue, IM, Index_Equ_IM, Equation_gr, 0, 0);
+ normalization.Gr_to_IM_basic(n, *prologue, *epilogue, IM, Equation_gr, false);
}
- else
- normalization.Gr_to_IM_basic(n, *prologue, *epilogue, IM, Equation_gr, false);
cout << "Finding the optimal block decomposition of the model ...\n";
if(bt_verbose)
blocks.Print_Equation_gr(Equation_gr);
res = blocks.sc(Equation_gr);
+ normalization.Free_Equation(n-*prologue-*epilogue,Equation_gr);
+ free(Equation_gr);
if(bt_verbose)
blocks.block_result_print(res);
if ((*prologue) || (*epilogue))
@@ -775,8 +837,10 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
else
j = 0;
for(i = 0;i < res->n_sets;i++)
- if ((res->sets_f[i] - res->sets_s[i] + 1) > j)
- j = res->sets_f[i] - res->sets_s[i] + 1;
+ {
+ if ((res->sets_f[i] - res->sets_s[i] + 1) > j)
+ j = res->sets_f[i] - res->sets_s[i] + 1;
+ }
Nb_RecursBlocks = *prologue + *epilogue;
Nb_TotalBlocks = res->n_sets + Nb_RecursBlocks;
cout << Nb_TotalBlocks << " block(s) found:\n";
@@ -790,8 +854,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
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);
+ //Free_IM(p_First_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++)
@@ -803,6 +868,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
}
if (*epilogue)
Allocate_Block(*epilogue, &count_Equ, &count_Block, EPILOGUE, ModelBlock);
+ blocks.block_result_free(res);
return 0;
}
@@ -852,6 +918,8 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_
for(i = 0;i < endo_nbr*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);
+ free(SIM_0);
+ free(SIM);
if(bt_verbose)
for(i = 0;i < endo_nbr;i++)
cout << "Block=" << Index_Equ_IM[i].block << " Equ=" << Index_Equ_IM[i].index << " Var= " << Index_Var_IM[i].index << " " << symbol_table.getNameByID(eEndogenous, Index_Var_IM[i].index) << "\n";
diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index f139a019216f3fbb54a442e5a0ea8f27b7837ea1..26f8fff991b3c2380381bd7a5f3311f599652af5 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -68,6 +68,23 @@ CheckStatement::checkPass(ModFileStructure &mod_file_struct)
mod_file_struct.check_present = true;
}
+Model_InfoStatement::Model_InfoStatement(const OptionsList &options_list_arg) :
+ options_list(options_list_arg)
+{
+}
+
+void Model_InfoStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+ //mod_file_struct.model_info_present = true;
+}
+
+void Model_InfoStatement::writeOutput(ostream &output, const string &basename) const
+{
+ options_list.writeOutput(output);
+ output << "model_info();\n";
+}
+
+
SimulStatement::SimulStatement(const OptionsList &options_list_arg) :
options_list(options_list_arg)
{
diff --git a/DynareBison.yy b/DynareBison.yy
index 85c0c086508524c31609701a9178c10507e1c58d..2fee306ccab465a99371bae0c110eac33a342d4f 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -183,6 +183,7 @@ statement : declaration
| dynatype
| dynasave
| model_comparison
+ | model_info
| planner_objective
| ramsey_policy
| bvar_density
@@ -612,6 +613,10 @@ check_options_list : check_options_list COMMA check_options
check_options : o_solve_algo;
+model_info : MODEL_INFO ';'
+ { driver.model_info(); }
+ ;
+
simul : SIMUL ';'
{ driver.simulate(); }
| SIMUL '(' simul_options_list ')' ';'
diff --git a/DynareFlex.ll b/DynareFlex.ll
index e222fe5df7d118f486202c2f24935689226920a1..6a2ebd90feb7aeb6f5d8a61e5f639a1b82439d5f 100644
--- a/DynareFlex.ll
+++ b/DynareFlex.ll
@@ -104,6 +104,7 @@ int sigma_e = 0;
<INITIAL>periods {BEGIN DYNARE_STATEMENT; return token::PERIODS;}
<INITIAL>cutoff {BEGIN DYNARE_STATEMENT; return token::CUTOFF;}
<INITIAL>markowitz {BEGIN DYNARE_STATEMENT; return token::MARKOWITZ;}
+<INITIAL>model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;}
<INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
<INITIAL>prior_analysis {BEGIN DYNARE_STATEMENT; return token::PRIOR_ANALYSIS;}
<INITIAL>posterior_analysis {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_ANALYSIS;}
diff --git a/DynareMain.cc b/DynareMain.cc
index 61dc9bbb808177954a180073c46f63d2b56afcc5..7a32a3a0afca6f08deef02c884608c703c396932 100644
--- a/DynareMain.cc
+++ b/DynareMain.cc
@@ -86,6 +86,7 @@ main(int argc, char** argv)
}
// Do the rest
+
main2(macro_output, basename, debug, clear_all);
return 0;
diff --git a/DynareMain2.cc b/DynareMain2.cc
index 53a5584d2219384b0cda3f24777f41c5291fb027..7061e71e266f4a66d9c502f7ff937b8d8e697d94 100644
--- a/DynareMain2.cc
+++ b/DynareMain2.cc
@@ -28,10 +28,10 @@ void
main2(stringstream &in, string &basename, bool debug, bool clear_all)
{
ParsingDriver p;
-
+ //cout << "OK\n";
// Do parsing and construct internal representation of mod file
ModFile *mod_file = p.parse(in, debug);
-
+ //cout << "OK1\n";
// Run checking pass
mod_file->checkPass();
diff --git a/ModFile.cc b/ModFile.cc
index b09b789a27f421b754aa355eea9c57053941bb8f..5a9699e01fb6a8cca6ec15edef0379c2452f7344 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -150,6 +150,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
<< "warning warning_old_state" << endl;
mOutputFile << "logname_ = '" << basename << ".log';" << endl;
mOutputFile << "diary " << basename << ".log" << endl;
+ mOutputFile << "options_.model_mode = " << model_tree.mode << ";\n";
if (model_tree.equation_number() > 0)
diff --git a/ModelBlocks.cc b/ModelBlocks.cc
index f37117a2ae1ff9aea999db6f46862ceee01fb2bf..a7a6ea12bb47343bcac3e042b90b52d2a6e6d8d4 100644
--- a/ModelBlocks.cc
+++ b/ModelBlocks.cc
@@ -207,6 +207,9 @@ Blocks::sc(Equation_set *g)
void
Blocks::block_result_free(block_result_t *r)
{
+#ifdef DEBUG
+ cout << "block_result_free\n";
+#endif
free(r->vertices);
free(r->sets_s);
free(r->sets_f);
diff --git a/ModelNormalization.cc b/ModelNormalization.cc
index fbd8cdb05ba8594c1db8a30c5cded222d67f6922..75dfe62b57eaf47b53430c524c168af38b91e605 100644
--- a/ModelNormalization.cc
+++ b/ModelNormalization.cc
@@ -283,9 +283,23 @@ void
Normalization::Gr_to_IM_basic(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation, bool transpose)
{
int i, j, edges, n;
- Edge *e1;
+ Edge *e1, *e2;
n = n0 - prologue - epilogue;
Equation->size = n;
+ if(Equation->Number)
+ {
+ for(i = 0;i < n;i++)
+ {
+ e1 = Equation->Number[i].First_Edge;
+ while(e1 != NULL)
+ {
+ e2 = e1->next;
+ free(e1);
+ e1 = e2;
+ }
+ }
+ free(Equation->Number);
+ }
Equation->Number = (Equation_vertex*)malloc(n * sizeof(Equation_vertex));
edges = 0;
if(transpose)
@@ -382,8 +396,12 @@ Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* In
}
free(SIM);
free(Index_Equ_IM_tmp);
+ //cout << "mixing=" << mixing << "\n";
if(mixing)
- Gr_to_IM_basic(n0, prologue, epilogue, IM, Equation, true);
+ {
+ //Free_Equation(n,Equation);
+ Gr_to_IM_basic(n0, prologue, epilogue, IM, Equation, true);
+ }
else
{
// In this step we :
@@ -453,17 +471,21 @@ Normalization::Free_Equation(int n, Equation_set* Equation)
while(e1 != NULL)
{
e2 = e1->next;
+ free(e1);
e1 = e2;
}
}
free(Equation->Number);
- free(Equation);
+ //free(Equation);
}
void
Normalization::Free_Other(Variable_set* Variable)
{
//free unused space
+#ifdef DEBUG
+ cout << "Free_Other\n";
+#endif
free(Local_Heap);
free(Variable->Number);
free(Variable);
@@ -526,6 +548,7 @@ Normalization::Normalize(int n, int prologue, int epilogue, bool* IM, simple* In
}
}
Free_Other(Variable);
+ //Free_All(n,Equation,Variable);
#ifdef DEBUG
cout << "end of Normalize\n";
#endif
diff --git a/ModelTree.cc b/ModelTree.cc
index a50e0932d862beffe7920272fd5ed7e94156907f..efd3661e2745769a99ceb726b13a297b8fc1d3f6 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -335,7 +335,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
if (order == 2)
for(second_derivatives_type::iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
- it->second->computeTemporaryTerms(reference_count, temporary_terms, false);
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, 0, ModelBlock, map_idx);
/*New*/
j=0;
for(temporary_terms_type::const_iterator it = temporary_terms.begin();
@@ -391,7 +391,7 @@ ModelTree::writeModelEquationsOrdered_C(ostream &output, Model_Block *ModelBlock
}
else
lhs_rhs_done=false;
- if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
&& (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
@@ -593,7 +593,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
bool OK, lhs_rhs_done, skip_the_head;
ostringstream Uf[symbol_table.endo_nbr];
map<NodeID, int> reference_count;
- int prev_Simulation_Type=-1;
+ int prev_Simulation_Type=-1, count_derivates=0;
temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//----------------------------------------------------------------------
//Temporary variables declaration
@@ -610,12 +610,9 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
/*tmp_output << "[" << block_triangular.periods + variable_table.max_lag+variable_table.max_lead << "]";*/
}
-
- if (tmp_output.str().length()>0)
- {
- global_output << " global " << tmp_output.str() << " M_ ;\n";
- }
+ global_output << " global " << tmp_output.str() << " M_ ;\n";
//For each block
+ int gen_blocks=0;
for(j = 0;j < ModelBlock->Size;j++)
{
//For a block composed of a single equation determines wether we have to evaluate or to solve the equation
@@ -630,7 +627,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
}
else
lhs_rhs_done=false;
- if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
&& (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
@@ -640,6 +637,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
skip_the_head=false;
if (!skip_the_head)
{
+ count_derivates=0;
+ gen_blocks++;
if (j>0)
{
output << "return;\n\n\n";
@@ -650,12 +649,12 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R)
- output << "function [y] = " << dynamic_basename << "_" << j+1 << "(y, x, it_)\n";
+ output << "function [y, g1, g2, g3] = " << 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_FOREWARD_SIMPLE)
- output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_)\n";
+ output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, g1, g2, g3, y_index, jacobian_eval)\n";
else
- output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods)\n";
+ output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods, g1, g2, g3)\n";
output << " % ////////////////////////////////////////////////////////////////////////" << endl
<< " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type)
<< " //" << endl
@@ -676,8 +675,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
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 << " 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 << " 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";
@@ -706,8 +705,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
{
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] << " variable : " << sModel
- << " (" << ModelBlock->Block_List[j].Variable[i] << ")" << endl;
+ output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
+ << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
if (!lhs_rhs_done)
{
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
@@ -762,101 +761,142 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
- && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R)
- {
+ && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R
+ && ModelBlock->Block_List[j].Simulation_Type!=SOLVE_FOREWARD_SIMPLE
+ && ModelBlock->Block_List[j].Simulation_Type!=SOLVE_BACKWARD_SIMPLE)
output << " " << sps << "% Jacobian " << endl;
- switch(ModelBlock->Block_List[j].Simulation_Type)
+ else
+ {
+ output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl;
+ }
+ switch(ModelBlock->Block_List[j].Simulation_Type)
+ {
+ case SOLVE_BACKWARD_SIMPLE:
+ case SOLVE_FOREWARD_SIMPLE:
+ case EVALUATE_BACKWARD:
+ case EVALUATE_FOREWARD:
+ case EVALUATE_BACKWARD_R:
+ case EVALUATE_FOREWARD_R:
+ count_derivates++;
+ for(m=0;m<ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
{
- case SOLVE_BACKWARD_SIMPLE:
- case SOLVE_FOREWARD_SIMPLE:
- output << " g1(1)=";
+ k=m-ModelBlock->Block_List[j].Max_Lag;
+ for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
+ {
+ 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 << ")=";
+ writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, 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;
+ }
+ }
+ }
+ if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
+ || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FOREWARD_SIMPLE)
+ {
+ output << " else\n";
+ output << " g1=";
writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, 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]
- << ", equation=" << ModelBlock->Block_List[j].Equation[0] << endl;
- break;
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FOREWARD_COMPLETE:
- m=ModelBlock->Block_List[j].Max_Lag;
+ << ") " << ModelBlock->Block_List[j].Variable[0]+1
+ << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
+ }
+ if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
+ || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
+ || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
+ || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R
+ || ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
+ || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FOREWARD_SIMPLE)
+ output << " end;" << endl;
+ break;
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FOREWARD_COMPLETE:
+ m=ModelBlock->Block_List[j].Max_Lag;
+ for(i=0;i<ModelBlock->Block_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];
+ Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")";
+ output << " u(" << u+1 << ") = ";
+ writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms);
+ 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:
+ 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;i<ModelBlock->Block_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")";
- output << " u(" << u+1 << ") = ";
- writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms);
+ int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
+ if (k==0)
+ Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")";
+ else if (k==1)
+ Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_)*y(it_+1, " << var+1 << ")";
+ else if (k>0)
+ Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")";
+ else if (k<0)
+ 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_) = ";
+ else if(k==1)
+ 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 << ")) = ";
+ 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 << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
- << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var
- << ", equation=" << eq << 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;i<ModelBlock->Block_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];
- if (k==0)
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")";
- else if (k>0)
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k << "-1))*y(it_+" << k << ", " << var+1 << ")";
- else if (k<0)
- 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_) = ";
- else if(k>0)
- 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 << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
- << "(" << k << ") " << var
- << ", equation=" << eq << endl;
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
#ifdef CONDITION
- output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
- output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
+ output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
+ output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
#endif
- }
}
- for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\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 << "))<fabs(u(" << i << "+Per_u_)))\n";
- output << " condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
+ output << " if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
+ output << " condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
#endif
- }
+ }
#ifdef CONDITION
- for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+ for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+ {
+ k=m-ModelBlock->Block_List[j].Max_Lag;
+ for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for(i=0;i<ModelBlock->Block_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";
- }
+ 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
- output << " end;\n";
- break;
}
+ 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
+ output << " end;\n";
+ break;
}
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
}
@@ -890,10 +930,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
tmp_output << " ";
(*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms);
}
- if (tmp_output.str().length()>0)
- {
- global_output << " global " << tmp_output.str() << " M_ ;\n";
- }
+ global_output << " global " << tmp_output.str() << " M_ ;\n";
//For each block
for(j = 0;j < ModelBlock->Size;j++)
{
@@ -909,7 +946,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
}
else
lhs_rhs_done=false;
- if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
&& (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
@@ -925,7 +962,13 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
}
else
output << "\n\n";
- output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n";
+ if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
+ ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
+ ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
+ ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R )
+ output << "function [y] = " << static_basename << "_" << j+1 << "(y, x)\n";
+ else
+ output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n";
output << " % ////////////////////////////////////////////////////////////////////////" << endl
<< " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " "
<< BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl
@@ -939,43 +982,36 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
output << " end\n";
}
-
temporary_terms_type tt2;
int n=ModelBlock->Block_List[j].Size;
int n1=symbol_table.endo_nbr;
- //cout << "n1=" << n1 << "\n";
IM=(bool*)malloc(n*n*sizeof(bool));
memset(IM, 0, n*n*sizeof(bool));
- //cout << "ModelBlock->Block_List[j].Max_Lead" << ModelBlock->Block_List[j].Max_Lag << " ModelBlock->Block_List[j].Max_Lag=" << ModelBlock->Block_List[j].Max_Lead << "\n";
for(m=-ModelBlock->Block_List[j].Max_Lag;m<=ModelBlock->Block_List[j].Max_Lead;m++)
{
- //cout << "bGet_IM(" << m << ")\n";
IMl=block_triangular.bGet_IM(m);
- //cout <<"OK\n";
for(i=0;i<n;i++)
{
eq=ModelBlock->Block_List[j].Equation[i];
for(k=0;k<n;k++)
{
var=ModelBlock->Block_List[j].Variable[k];
- //cout << "eq=" << eq << " var=" << var << "\n";
IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var];
- /*if(IM[i*n+k])
- cout << " ->i=" << i << " j=" << j << "\n";*/
}
}
- //cout << "done\n";
}
for(nze=0, i=0;i<n*n;i++)
{
nze+=IM[i];
}
- cout << "nze=" << nze << "\n";
memset(IM, 0, n*n*sizeof(bool));
if( ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R)
- output << " g1=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
+ {
+ output << " g1=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
+ output << " residual=zeros(" << ModelBlock->Block_List[j].Size << ",1);\n";
+ }
sps="";
if (ModelBlock->Block_List[j].Temporary_terms->size())
output << " " << sps << "% //Temporary variables" << endl;
@@ -997,8 +1033,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
{
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] << " variable : "
- << sModel << " (" << ModelBlock->Block_List[j].Variable[i] << ")" << endl;
+ output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : "
+ << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
if (!lhs_rhs_done)
{
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
@@ -1057,8 +1093,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
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]
- << ", equation=" << ModelBlock->Block_List[j].Equation[0] << endl;
+ << ") " << ModelBlock->Block_List[j].Variable[0]+1
+ << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
break;
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FOREWARD_COMPLETE:
@@ -1073,8 +1109,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
output << " u(" << u+1 << ") = ";
writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
- << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var
- << ", equation=" << eq << endl;
+ << "(" << 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";
@@ -1091,8 +1127,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
//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";
- cout << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n";
+ //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
@@ -1102,8 +1137,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
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
- << ", equation=" << eq << endl;
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
#ifdef CONDITION
output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
@@ -1138,9 +1173,10 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
}
}
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
+ free(IM);
}
output << "return;\n\n\n";
- free(IM);
+ //free(IM);
}
@@ -1190,7 +1226,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
ModelBlock_Aggregated_Count=-1;
for(j = 0;j < ModelBlock->Size;j++)
{
- if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
&& (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
@@ -1266,10 +1302,10 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
}
else
lhs_rhs_done=false;
- if (ModelBlock->Block_List[j].Size==1)
+ /*if (ModelBlock->Block_List[j].Size==1)
lhs_rhs_done=true;
else
- lhs_rhs_done=false;
+ lhs_rhs_done=false;*/
//The Temporary terms
temporary_terms_type tt2;
i=0;
@@ -1994,7 +2030,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
- cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n";
+ //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n";
u_count_int++;
}
}
@@ -2008,7 +2044,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
- cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n";
+ //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n";
u_count_int++;
}
for(j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
@@ -2052,9 +2088,10 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
mStaticModelFile << " if(length(varargin)>0)\n";
mStaticModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
mStaticModelFile << " global it_;\n";
- mStaticModelFile << " y=varargin{1}(y_kmin,:);\n";
+ mStaticModelFile << " y=varargin{1}(:);\n";
mStaticModelFile << " ys=y;\n";
- mStaticModelFile << " x=varargin{2}(y_kmin,:);\n";
+ mStaticModelFile << " g1=[];\n";
+ mStaticModelFile << " x=varargin{2}(:);\n";
mStaticModelFile << " residual=zeros(1, " << symbol_table.endo_nbr << ");\n";
prev_Simulation_Type=-1;
for(i=0;i<block_triangular.ModelBlock->Size;i++)
@@ -2066,7 +2103,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
}
mStaticModelFile << " ];\n";
k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if (prev_Simulation_Type==k &&
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
skip_head=true;
else
@@ -2078,7 +2115,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
case EVALUATE_FOREWARD_R:
case EVALUATE_BACKWARD_R:
if(!skip_head)
- mStaticModelFile << " " << static_basename << "_" << i + 1 << "(y, x);\n";
+ mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n";
mStaticModelFile << " residual(y_index)=ys(y_index)-y(y_index);\n";
break;
case SOLVE_FOREWARD_COMPLETE:
@@ -2092,168 +2129,152 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
}
prev_Simulation_Type=k;
}
- mStaticModelFile << " varargout{1}=residual;\n";
- mStaticModelFile << " varargout{2}=g1;\n";
- mStaticModelFile << " return;\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " %it is the deterministic simulation of the block decomposed static model\n";
- mStaticModelFile << " periods=options_.periods;\n";
- mStaticModelFile << " maxit_=options_.maxit_;\n";
- mStaticModelFile << " solve_tolf=options_.solve_tolf;\n";
- mStaticModelFile << " y=oo_.steady_state;\n";
- mStaticModelFile << " x=oo_.exo_steady_state;\n";
- prev_Simulation_Type=-1;
- for(i = 0;i < block_triangular.ModelBlock->Size;i++)
- {
- k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if (prev_Simulation_Type==k &&
- (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
- skip_head=true;
- else
- skip_head=false;
- if ((k == EVALUATE_FOREWARD || k == EVALUATE_FOREWARD_R || k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (!skip_head)
- {
- if (open_par)
- {
- mStaticModelFile << " end\n";
- }
- mStaticModelFile << " " << static_basename << "_" << i + 1 << "(y, x);\n";
- }
- open_par=false;
- }
+ mStaticModelFile << " varargout{1}=residual;\n";
+ mStaticModelFile << " varargout{2}=g1;\n";
+ mStaticModelFile << " return;\n";
+ mStaticModelFile << " end;\n";
+ mStaticModelFile << " %it is the deterministic simulation of the block decomposed static model\n";
+ mStaticModelFile << " periods=options_.periods;\n";
+ mStaticModelFile << " maxit_=options_.maxit_;\n";
+ mStaticModelFile << " solve_tolf=options_.solve_tolf;\n";
+ mStaticModelFile << " y=oo_.steady_state;\n";
+ mStaticModelFile << " x=oo_.exo_steady_state;\n";
+ mStaticModelFile << " varargout{2}=1;\n";
+ prev_Simulation_Type=-1;
+ for(i = 0;i < block_triangular.ModelBlock->Size;i++)
+ {
+ k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
+ (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
+ skip_head=true;
+ else
+ skip_head=false;
+ if ((k == EVALUATE_FOREWARD || k == EVALUATE_FOREWARD_R || k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (!skip_head)
+ {
+ if (open_par)
+ {
+ mStaticModelFile << " end\n";
+ }
+ mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n";
+ }
+ open_par=false;
+ }
else if ((k == SOLVE_FOREWARD_SIMPLE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
+ {
+ if (open_par)
+ {
mStaticModelFile << " end\n";
- }
- 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])<solve_tolf);\n";
- mStaticModelFile << " iter=iter+1;\n";
- mStaticModelFile << " end\n";
- mStaticModelFile << " if cvg==0\n";
- mStaticModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
- mStaticModelFile << " return;\n";
- mStaticModelFile << " end\n";
- mStaticModelFile << " end\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);\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)<solve_tolf);\n";
- mStaticModelFile << " iter=iter+1;\n";
- mStaticModelFile << " end\n";
- mStaticModelFile << " if cvg==0\n";
- mStaticModelFile << " fprintf('Convergence not achieved in block " << i << ", after %d iterations\\n',iter);\n";
- mStaticModelFile << " return;\n";
- mStaticModelFile << " end\n";
- }
- else if ((k == SOLVE_FOREWARD_COMPLETE || k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_TWO_BOUNDARIES_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- mStaticModelFile << "end\n";
- }
- open_par=false;
- mStaticModelFile << " y_index=[";
- for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
- {
- mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
- }
- mStaticModelFile << " ];\n";
- mStaticModelFile << " g1=0;g2=0;g3=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, 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)<solve_tolf);\n";
- mStaticModelFile << " iter=iter+1;\n";
- mStaticModelFile << " end\n";
- mStaticModelFile << " if cvg==0\n";
- mStaticModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
- mStaticModelFile << " return;\n";
- mStaticModelFile << " end\n";
- mStaticModelFile << " end\n";*/
- mStaticModelFile << " cvg=0;\n";
- mStaticModelFile << " iter=0;\n";
- /*mStaticModelFile << " Per_y_=it_*y_size;\n";*/
- mStaticModelFile << " lambda=1;\n";
- mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n";
- mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n";
- mStaticModelFile << " max_res=max(abs(r));\n";
- mStaticModelFile << " if(iter>0)\n";
- mStaticModelFile << " if(~isreal(max_res) | max_resa<max_res)\n";
- mStaticModelFile << " if(lambda>1e-6)\n";
- mStaticModelFile << " lambda=lambda/2;\n";
- mStaticModelFile << " y(y_index)=y_save+lambda*dx;\n";
- mStaticModelFile << " continue;\n";
- mStaticModelFile << " else\n";
- mStaticModelFile << " disp(['No convergence after ' num2str(iter,'%d') ' iterations']);\n";
- mStaticModelFile << " return;\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " else\n";
- mStaticModelFile << " if(lambda<1)\n";
- mStaticModelFile << " lambda=max(lambda*2, 1);\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " max_resa=max_res;\n";
- mStaticModelFile << " cvg=(max_res<solve_tolf);\n";
- mStaticModelFile << " if (cvg==0),\n";
- mStaticModelFile << " spparms('autommd',0);\n";
- mStaticModelFile << " q = colamd(g1);\n";
- mStaticModelFile << " z = g1(:,q)\\b';\n";
- mStaticModelFile << " z(q) = z;\n";
- mStaticModelFile << " spparms('autommd',1);\n";
- mStaticModelFile << " y_save=y(y_index);\n";
- mStaticModelFile << " dx= (z-y_save);\n";
- mStaticModelFile << " y(y_index)=y_save+lambda*dx;\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " iter=iter+1;\n";
- mStaticModelFile << " disp(['iter=' num2str(iter,'%d') ' err=' num2str(max_res,'%f')]);\n";
- mStaticModelFile << " end\n";
- mStaticModelFile << " if cvg==0\n";
- mStaticModelFile << " fprintf('Error in steady: Convergence not achieved in block " << i << ", after %d iterations\\n',iter);\n";
- mStaticModelFile << " return;\n";
- mStaticModelFile << " else\n";
- mStaticModelFile << " fprintf('convergence achieved after %d iterations\\n',iter);\n";
- mStaticModelFile << " end\n";
- }
- prev_Simulation_Type=k;
- }
- if(open_par)
- mStaticModelFile << " end;\n";
- mStaticModelFile << " oo_.steady_state = y;\n";
- mStaticModelFile << " if isempty(ys0_)\n";
- mStaticModelFile << " oo_.endo_simul(:,1:M_.maximum_lag) = oo_.steady_state * ones(1,M_.maximum_lag);\n";
- mStaticModelFile << " else\n";
- mStaticModelFile << " options_ =set_default_option(options_,'periods',1);\n";
- mStaticModelFile << " oo_.endo_simul(:,M_.maximum_lag+1:M_.maximum_lag+options_.periods+M_.maximum_lead) = oo_.steady_state * ones(1,options_.periods+M_.maximum_lead);\n";
- mStaticModelFile << " end;\n";
- mStaticModelFile << " disp('Steady State value');\n";
- mStaticModelFile << " disp([strcat(M_.endo_names,' : ') num2str(oo_.steady_state,'%f')]);\n";
- mStaticModelFile << "return;\n";
-
- writeModelStaticEquationsOrdered_M(mStaticModelFile, block_triangular.ModelBlock, static_basename);
- mStaticModelFile.close();
+ }
+ 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])<solve_tolf);\n";
+ mStaticModelFile << " iter=iter+1;\n";
+ mStaticModelFile << " end\n";
+ mStaticModelFile << " if cvg==0\n";
+ mStaticModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
+ mStaticModelFile << " return;\n";
+ mStaticModelFile << " end\n";
+ mStaticModelFile << " end\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);\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)<solve_tolf);\n";
+ mStaticModelFile << " iter=iter+1;\n";
+ mStaticModelFile << " end\n";
+ mStaticModelFile << " if cvg==0\n";
+ mStaticModelFile << " fprintf('Convergence not achieved in block " << i << ", after %d iterations\\n',iter);\n";
+ mStaticModelFile << " return;\n";
+ mStaticModelFile << " end\n";
+ }
+ else if ((k == SOLVE_FOREWARD_COMPLETE || k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_TWO_BOUNDARIES_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ {
+ mStaticModelFile << "end\n";
+ }
+ open_par=false;
+ mStaticModelFile << " y_index=[";
+ for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+ {
+ mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ }
+ mStaticModelFile << " ];\n";
+ mStaticModelFile << " g1=0;g2=0;g3=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, 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)<solve_tolf);\n";
+ mStaticModelFile << " iter=iter+1;\n";
+ mStaticModelFile << " end\n";
+ mStaticModelFile << " if cvg==0\n";
+ mStaticModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
+ mStaticModelFile << " return;\n";
+ mStaticModelFile << " end\n";
+ mStaticModelFile << " end\n";*/
+ mStaticModelFile << " cvg=0;\n";
+ mStaticModelFile << " iter=0;\n";
+ /*mStaticModelFile << " Per_y_=it_*y_size;\n";*/
+ mStaticModelFile << " lambda=1;\n";
+ mStaticModelFile << " stpmx = 100 ;\n";
+ mStaticModelFile << " stpmax = stpmx*max([sqrt(y'*y);size(y_index,2)]);\n";
+ mStaticModelFile << " nn=1:size(y_index,2);\n";
+ mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n";
+ mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n";
+ mStaticModelFile << " max_res=max(abs(r));\n";
+ mStaticModelFile << " cvg=(max_res<solve_tolf);\n";
+ mStaticModelFile << " if (cvg==0),\n";
+ mStaticModelFile << " g = (r'*g1)';\n";
+ mStaticModelFile << " f = 0.5*r'*r;\n";
+ mStaticModelFile << " p = -g1\\r ;\n";
+ mStaticModelFile << " [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,@" << static_basename << "_" << i + 1 << ",nn,y_index,x);\n";
+ mStaticModelFile << " end;\n";
+ mStaticModelFile << " iter=iter+1;\n";
+ mStaticModelFile << " disp(['iter=' num2str(iter,'%d') ' err=' num2str(max_res,'%f')]);\n";
+ mStaticModelFile << " end\n";
+ mStaticModelFile << " if cvg==0\n";
+ mStaticModelFile << " fprintf('Error in steady: Convergence not achieved in block " << i << ", after %d iterations\\n',iter);\n";
+ mStaticModelFile << " return;\n";
+ mStaticModelFile << " else\n";
+ mStaticModelFile << " fprintf('convergence achieved after %d iterations\\n',iter);\n";
+ mStaticModelFile << " end\n";
+ }
+ prev_Simulation_Type=k;
+ }
+ if(open_par)
+ mStaticModelFile << " end;\n";
+ mStaticModelFile << " oo_.steady_state = y;\n";
+ /*mStaticModelFile << " if isempty(ys0_)\n";
+ mStaticModelFile << " oo_.endo_simul(:,1:M_.maximum_lag) = oo_.steady_state * ones(1,M_.maximum_lag);\n";
+ mStaticModelFile << " else\n";
+ mStaticModelFile << " options_ =set_default_option(options_,'periods',1);\n";
+ mStaticModelFile << " oo_.endo_simul(:,M_.maximum_lag+1:M_.maximum_lag+options_.periods+M_.maximum_lead) = oo_.steady_state * ones(1,options_.periods+M_.maximum_lead);\n";
+ mStaticModelFile << " end;\n";*/
+ mStaticModelFile << " disp('Steady State value');\n";
+ mStaticModelFile << " disp([strcat(M_.endo_names,' : ') num2str(oo_.steady_state,'%f')]);\n";
+ mStaticModelFile << " varargout{2}=0;\n";
+ mStaticModelFile << " varargout{1}=oo_.steady_state;\n";
+ mStaticModelFile << "return;\n";
+ writeModelStaticEquationsOrdered_M(mStaticModelFile, block_triangular.ModelBlock, static_basename);
+ mStaticModelFile.close();
}
@@ -2264,6 +2285,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
{
string filename, sp;
ofstream mDynamicModelFile;
+ ostringstream tmp, tmp1, tmp_eq;
int prev_Simulation_Type;
SymbolicGaussElimination SGE;
bool OK;
@@ -2342,6 +2364,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
{
mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
mDynamicModelFile << " global oo_ options_ M_ ;\n";
+ mDynamicModelFile << " g2=[];g3=[];\n";
//Temporary variables declaration
{
OK=true;
@@ -2383,10 +2406,15 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << " y_size=M_.endo_nbr;\n";
mDynamicModelFile << " if(length(varargin)>0)\n";
mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
- mDynamicModelFile << " global it_;\n";
+ //mDynamicModelFile << " global it_;\n";
+ mDynamicModelFile << " it_=varargin{3};\n";
+ //mDynamicModelFile << " g=zeros(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
+ mDynamicModelFile << " g1=spalloc(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),y_size*y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
mDynamicModelFile << " Per_u_=0;\n";
mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " y1=varargin{1};\n";
+ mDynamicModelFile << " y=varargin{1};\n";
+ mDynamicModelFile << " ys=y(it_,:);\n";
+ /*mDynamicModelFile << " y1=varargin{1};\n";
mDynamicModelFile << " cnb_nz_elem=1;\n";
mDynamicModelFile << " for i = -y_kmin:y_kmax\n";
mDynamicModelFile << " nz_elem=find(M_.lead_lag_incidence(:,1+i+y_kmin));\n";
@@ -2397,19 +2425,47 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << " nz_elem_s=nz_elem;\n";
mDynamicModelFile << " end;\n";
mDynamicModelFile << " cnb_nz_elem=cnb_nz_elem+nb_nz_elem;\n";
- mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " end;\n";*/
mDynamicModelFile << " x=varargin{2};\n";
prev_Simulation_Type=-1;
+ tmp.str("");
+ tmp_eq.str("");
for(i = 0;i < block_triangular.ModelBlock->Size;i++)
{
- mDynamicModelFile << " y_index=[";
+
+ k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
+ if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
+ ((prev_Simulation_Type==EVALUATE_FOREWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FOREWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R)
+ || (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)))
+ {
+ mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
+ tmp_eq.str("");
+ mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
+ tmp.str("");
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
+ }
+ //mDynamicModelFile << " y_index=[";
for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
}
- mDynamicModelFile << " ];\n";
- k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if (prev_Simulation_Type==k &&
+ //mDynamicModelFile << " ];\n";
+ if(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R
+ )
+ {
+ if(i==block_triangular.ModelBlock->Size-1)
+ {
+ mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
+ tmp_eq.str("");
+ mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
+ tmp.str("");
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
+ }
+ }
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
skip_head=true;
else
@@ -2421,19 +2477,49 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
case EVALUATE_FOREWARD_R:
case EVALUATE_BACKWARD_R:
if(!skip_head)
- mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_kmin, Per_u_, Per_y_, y_size);\n";
- mDynamicModelFile << " residual(y_index)=ys(y_index)-y(it_, y_index);\n";
+ {
+ tmp1 << " [y, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n";
+ tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
+ }
+ break;
+ case SOLVE_FOREWARD_SIMPLE:
+ case SOLVE_BACKWARD_SIMPLE:
+ 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 << " residual(y_index_eq)=r;\n";
break;
case SOLVE_FOREWARD_COMPLETE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
- mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, it_);\n";
- mDynamicModelFile << " residual(y_index)=r;\n";
+ //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, 1);\n";
+ mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
+ mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
+ tmp.str("");
+ tmp_eq.str("");
+ int tmp_i=variable_table.max_lag+variable_table.max_lead+1;
+ 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, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\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";
+ 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 << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
break;
}
prev_Simulation_Type=k;
}
- mDynamicModelFile << " varagout{1}=residual;\n";
+ if(tmp1.str().length())
+ {
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
+ }
+ mDynamicModelFile << " varargout{1}=residual;\n";
+ mDynamicModelFile << " varargout{2}=g1;\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n";
@@ -2459,7 +2545,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
for(i = 0;i < block_triangular.ModelBlock->Size;i++)
{
k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if (prev_Simulation_Type==k &&
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
skip_head=true;
else
@@ -2490,7 +2576,8 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << " Per_u_=0;\n";
mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, it_);\n";
+ mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
+ mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n";
}
}
if(mode==eSparseDLLMode)
@@ -2524,7 +2611,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << " Per_u_=0;\n";
mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_);\n";
+ mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3);\n";
}
}
if(mode==eSparseDLLMode)
@@ -2585,7 +2672,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
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_);\n";
+ mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\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)<solve_tolf);\n";
mDynamicModelFile << " iter=iter+1;\n";
@@ -2649,9 +2736,9 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
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_);\n";
- mDynamicModelFile << " y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r[it_]/g1;\n";
- mDynamicModelFile << " cvg=((r[it_]*r[it_])<solve_tolf);\n";
+ mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\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)<solve_tolf);\n";
mDynamicModelFile << " iter=iter+1;\n";
mDynamicModelFile << " end\n";
mDynamicModelFile << " if cvg==0\n";
@@ -2728,7 +2815,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << " {\n";
mDynamicModelFile << " Per_u_=(it_-y_kmin)*" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";\n";
mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, r, g1, g2);\n";
+ mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, r, g1, g2, g3);\n";
mDynamicModelFile << "#ifdef PRINT_OUT\n";
mDynamicModelFile << " for(j=0;j<" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";j++)\n";
mDynamicModelFile << " {\n";
@@ -2800,7 +2887,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
}
else
{
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
+ mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2, g3);\n";
mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", 0, false);\n";
}
mDynamicModelFile << " }\n";
@@ -2997,6 +3084,19 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << " y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n";
mDynamicModelFile << " lambda=options_.slowc;\n";
mDynamicModelFile << " correcting_factor=0.01;\n";
+ mDynamicModelFile << " luinc_tol=1e-10;\n";
+ mDynamicModelFile << " max_resa=1e100;\n";
+ int nze, m;
+ for(nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
+ 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=" ";
@@ -3006,7 +3106,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
{
sp="";
}
- mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods);\n";
+ mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\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";
@@ -3014,9 +3114,8 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << sp << " else\n";
mDynamicModelFile << sp << " max_res=max(max(abs(r)));\n";
mDynamicModelFile << sp << " end;\n";
-
mDynamicModelFile << sp << " if(iter>0)\n";
- mDynamicModelFile << sp << " if(~isreal(max_res) | isnan(max_res) | max_resa<max_res)\n";
+ mDynamicModelFile << sp << " if(~isreal(max_res) | isnan(max_res) | (max_resa<max_res && iter>1))\n";
mDynamicModelFile << sp << " if(isnan(max_res))\n";
mDynamicModelFile << sp << " detJ=det(g1aa);\n";
mDynamicModelFile << sp << " if(abs(detJ)<1e-7)\n";
@@ -3037,8 +3136,9 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << sp << " return;\n";
mDynamicModelFile << sp << " end;\n";
mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " elseif(lambda>1e-6)\n";
+ mDynamicModelFile << sp << " elseif(lambda>1e-8)\n";
mDynamicModelFile << sp << " lambda=lambda/2;\n";
+ mDynamicModelFile << sp << " reduced = 1;\n";
mDynamicModelFile << sp << " disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);\n";
mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';\n";
if (!block_triangular.ModelBlock->Block_List[i].is_linear)
@@ -3059,41 +3159,50 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << sp << " g1aa=g1a;\n";
mDynamicModelFile << sp << " ba=b;\n";
mDynamicModelFile << sp << " max_resa=max_res;\n";
-
-
mDynamicModelFile << sp << " if(options_.simulation_method==0),\n";
mDynamicModelFile << sp << " dx = g1a\\b- ya;\n";
mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
mDynamicModelFile << sp << " elseif(options_.simulation_method==2),\n";
- mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,1e-6);\n";
- mDynamicModelFile << sp << " [za,flag1] = gmres(g1a,b," << block_triangular.ModelBlock->Block_List[i].Size << ",1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
- mDynamicModelFile << sp << " dx = za - ya;\n";
- mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
- mDynamicModelFile << sp << " if (flag1>0)\n";
- mDynamicModelFile << sp << " if(flag1==1)\n";
- mDynamicModelFile << sp << " disp(['No convergence inside GMRES after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
- mDynamicModelFile << sp << " elseif(flag1==2)\n";
- mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
- mDynamicModelFile << sp << " elseif(flag1==3)\n";
- mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\n";
+ mDynamicModelFile << sp << " flag1=1;\n";
+ mDynamicModelFile << sp << " while(flag1>0)\n";
+ mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n";
+ mDynamicModelFile << sp << " [za,flag1] = gmres(g1a,b," << block_triangular.ModelBlock->Block_List[i].Size << ",1e-6," << 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 GMRES after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==2)\n";
+ mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==3)\n";
+ mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n";
+ mDynamicModelFile << sp << " reduced = 0;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " dx = za - ya;\n";
+ mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
mDynamicModelFile << sp << " end;\n";
mDynamicModelFile << sp << " end;\n";
mDynamicModelFile << sp << " elseif(options_.simulation_method==3),\n";
- mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,1e-7);\n";
- mDynamicModelFile << sp << " [za,flag1] = bicgstab(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
- mDynamicModelFile << sp << " dx = za - ya;\n";
- mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
- //mDynamicModelFile << sp << " [ya,flag1] = eval(strcat(fullfile(matlabroot,'toolbox','matlab','sparfun','bicgstab'),'(g1a,b,1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1)'));\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
- mDynamicModelFile << sp << " if (flag1>0)\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";
- mDynamicModelFile << sp << " elseif(flag1==2)\n";
- mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
- mDynamicModelFile << sp << " elseif(flag1==3)\n";
- mDynamicModelFile << sp << " disp(['BICGSTAB stagnated. (Two consecutive iterates were the same.)']);\n";
+ 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 << " 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";
+ mDynamicModelFile << sp << " elseif(flag1==2)\n";
+ mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==3)\n";
+ mDynamicModelFile << sp << " disp(['BICGSTAB stagnated. (Two consecutive iterates were the same.)']);\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n";
+ mDynamicModelFile << sp << " reduced = 0;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " dx = za - ya;\n";
+ mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
mDynamicModelFile << sp << " end;\n";
mDynamicModelFile << sp << " end;\n";
mDynamicModelFile << sp << " end;\n";
@@ -3509,7 +3618,171 @@ ModelTree::writeOutput(ostream &output) const
output << ";";
}
output << "]';\n";
+ //In case of sparse model, writes the block structure of the model
+ if(mode==eSparseMode || mode==eSparseDLLMode)
+ {
+ int prev_Simulation_Type=-1;
+ bool skip_the_head;
+ int k=0;
+ int count_lead_lag_incidence = 0;
+ int max_lead, max_lag;
+ 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
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j].Simulation_Type)
+ && (block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
+ ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
+ ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
+ ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R ))
+ skip_the_head=true;
+ else
+ {
+ skip_the_head=false;
+ k++;
+ count_lead_lag_incidence = 0;
+ int Block_size=block_triangular.ModelBlock->Block_List[j].Size;
+ max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ;
+ max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead;
+ bool evaluate=false;
+ ostringstream tmp_s, tmp_s_eq;
+ tmp_s.str("");
+ tmp_s_eq.str("");
+ for(int i=0;i<block_triangular.ModelBlock->Block_List[j].Size;i++)
+ {
+ tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1;
+ tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1;
+ }
+ if (block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
+ ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
+ ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
+ ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R
+ && j+Block_size<(block_triangular.ModelBlock->Size))
+ {
+ bool OK=true;
+ evaluate=true;
+ while(j+Block_size<(block_triangular.ModelBlock->Size) && OK)
+ {
+ if(BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j].Simulation_Type)!=BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j+Block_size].Simulation_Type))
+ OK=false;
+ else
+ {
+ if(max_lag <block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag )
+ max_lag =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag ;
+ if(max_lead<block_triangular.ModelBlock->Block_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";
+ for(int i=0;i<block_triangular.ModelBlock->Block_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;
+ }
+ 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 << ").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";
+ tmp_s.str("");
+ bool done_IM=false;
+ if(!evaluate)
+ {
+ output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n";
+ for(int l=-max_lag;l<max_lead+1;l++)
+ {
+ bool *tmp_IM;
+ tmp_IM=block_triangular.bGet_IM(l);
+ for(int l_var=0;l_var<block_triangular.ModelBlock->Block_List[j].Size;l_var++)
+ {
+ for(int l_equ=0;l_equ<block_triangular.ModelBlock->Block_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;l<max_lead+1;l++)
+ {
+ bool not_increm=true;
+ bool *tmp_IM;
+ tmp_IM=block_triangular.bGet_IM(l);
+ int ii=j;
+ while(ii-j<Block_size)
+ {
+ for(int l_var=0;l_var<block_triangular.ModelBlock->Block_List[ii].Size;l_var++)
+ {
+ for(int l_equ=0;l_equ<block_triangular.ModelBlock->Block_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";
+ done_IM=false;
+ }
+ ii++;
+ }
+ output << tmp_s.str() << "\n";
+ tmp_s.str("");
+ }
+ output << "];\n";
+ }
+ }
+ prev_Simulation_Type=block_triangular.ModelBlock->Block_List[j].Simulation_Type;
+
+ }
+ for(int j=-block_triangular.Model_Max_Lag;j<block_triangular.Model_Max_Lead+1;j++)
+ {
+ bool* IM = block_triangular.bGet_IM(j);
+ if(IM)
+ {
+ bool new_entry=true;
+ output << "M_.block_structure.incidence(" << block_triangular.Model_Max_Lag+j+1 << ").sparse_IM = [";
+ for(int i=0;i<symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
+ {
+ if(IM[i])
+ {
+ if(!new_entry)
+ output << " ; ";
+ else
+ output << " ";
+ output << i/symbol_table.endo_nbr+1 << " " << i % symbol_table.endo_nbr+1;
+ new_entry=false;
+ }
+ }
+ output << "];\n";
+ }
+ }
+ }
// Writing initialization for some other variables
output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr << "];\n";
output << "M_.maximum_lag = " << variable_table.max_lag << ";\n";
@@ -3594,7 +3867,8 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
}
if (IM[eq*symbol_table.endo_nbr+var] && (fabs(val) < cutoff))
{
- 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";
+ 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);
i++;
}
@@ -3742,6 +4016,8 @@ ModelTree::writeDynamicFile(const string &basename) const
break;
case eSparseMode:
writeSparseDynamicFileAndBinFile(basename + "_dynamic", basename, output_type, mode);
+ block_triangular.Free_Block(block_triangular.ModelBlock);
+ block_triangular.Free_IM(block_triangular.First_IM);
break;
case eDLLMode:
writeDynamicCFile(basename + "_dynamic");
@@ -3750,6 +4026,8 @@ ModelTree::writeDynamicFile(const string &basename) const
writeSparseDynamicFileAndBinFile(basename + "_dynamic", basename, output_type, mode);
if (compiler==GCC_COMPILE || compiler==LCC_COMPILE )
writeSparseDLLDynamicHFile(basename + "_dynamic");
+ block_triangular.Free_Block(block_triangular.ModelBlock);
+ block_triangular.Free_IM(block_triangular.First_IM);
break;
}
}
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index 1c498bed3371e7bd54e18c296a65058c779ad4d7..ac254e3e990d099456c6494c32f28b62f16d1a36 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -729,6 +729,15 @@ ParsingDriver::simul()
options_list.clear();
}
+
+void
+ParsingDriver::model_info()
+{
+ mod_file->addStatement(new Model_InfoStatement(options_list));
+ options_list.clear();
+}
+
+
void
ParsingDriver::check()
{
diff --git a/include/BlockTriangular.hh b/include/BlockTriangular.hh
index ceb96dda5e0807e8e9a96d4c1120e8cf09c09bce..3d31c7485318230c0851e0edebab82405b55066b 100644
--- a/include/BlockTriangular.hh
+++ b/include/BlockTriangular.hh
@@ -57,14 +57,14 @@ public:
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);
+ void Free_IM(List_IM* First_IM) const;
void Print_SIM(bool* IM, int n) const;
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);
+ 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, int type, Model_Block * ModelBlock);
- void Free_Block(Model_Block* ModelBlock);
+ void Free_Block(Model_Block* ModelBlock) const;
List_IM *First_IM ;
List_IM *Last_IM ;
simple *Index_Equ_IM;
diff --git a/include/ComputingTasks.hh b/include/ComputingTasks.hh
index aae31cf6d73bd870db6d76ec130f99830bc5ffa3..f497a45bf668aaf7e56acc5e370dbad581681b38 100644
--- a/include/ComputingTasks.hh
+++ b/include/ComputingTasks.hh
@@ -76,6 +76,16 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const;
};
+class Model_InfoStatement : public Statement
+{
+private:
+ const OptionsList options_list;
+public:
+ Model_InfoStatement(const OptionsList &options_list_arg);
+ virtual void checkPass(ModFileStructure &mod_file_struct);
+ virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
class StochSimulStatement : public Statement
{
private:
diff --git a/include/ModelNormalization.hh b/include/ModelNormalization.hh
index dc99b87c3cc3141e6b4ec604559a7c0276cda19c..47460b94a38fbacf7fbf562dba03364037e47a6d 100644
--- a/include/ModelNormalization.hh
+++ b/include/ModelNormalization.hh
@@ -73,6 +73,7 @@ public:
void Gr_to_IM_basic(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation,bool transpose);
const SymbolTable &symbol_table;
void Set_fp_verbose(bool ok);
+ void Free_Equation(int n, Equation_set* Equation);
private:
void IM_to_Gr(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation, Variable_set *Variable );
void Inits(Equation_set *Equation);
@@ -83,7 +84,6 @@ private:
int MeasureMatching(Equation_set *Equation);
void OutputMatching(Equation_set* Equation);
void Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* Index_Var_IM, Equation_set *Equation,bool mixing, bool* IM_s);
- void Free_Equation(int n, Equation_set* Equation);
void Free_Other(Variable_set* Variable);
void Free_All(int n, Equation_set* Equation, Variable_set* Variable);
int eq, eex;
diff --git a/include/ParsingDriver.hh b/include/ParsingDriver.hh
index d13dd5046d1399595a1d32937d3fe85a7290e478..c9da8d05bb4912fdce3c317cfdff6dbf6dda0cf9 100644
--- a/include/ParsingDriver.hh
+++ b/include/ParsingDriver.hh
@@ -46,6 +46,7 @@ using namespace std;
# undef yyFlexLexer
#endif
+
//! The lexer class
/*! Actually it was necessary to subclass the DynareFlexLexer class generated by Flex,
since the prototype for DynareFlexLexer::yylex() was not convenient.
@@ -53,7 +54,7 @@ using namespace std;
class DynareFlex : public DynareFlexLexer
{
public:
- DynareFlex(istream* in = 0, ostream* out = 0);
+ DynareFlex(std::istream* in = 0, ostream* out = 0);
//! The main lexing function
Dynare::parser::token_type lex(Dynare::parser::semantic_type *yylval,
@@ -145,7 +146,7 @@ public:
//! Starts parsing, and constructs the MOD file representation
/*! The returned pointer should be deleted after use */
- ModFile *parse(istream &in, bool debug);
+ ModFile *parse(std::istream &in, bool debug);
//! Reference to the lexer
class DynareFlex *lexer;
@@ -282,6 +283,8 @@ public:
void simul();
//! Writes check command
void check();
+ //! Writes model_info command
+ void model_info();
//! Writes estimated params command
void estimated_params();
//! Writes estimated params init command
diff --git a/macro/MacroDriver.hh b/macro/MacroDriver.hh
index cb8a29e25ea59fe62ad37fe0015e9bafc1abd95b..63bded30be0966d931b4bb305d5f363f7bd1c922 100644
--- a/macro/MacroDriver.hh
+++ b/macro/MacroDriver.hh
@@ -53,12 +53,12 @@ private:
class ScanContext
{
public:
- istream *input;
+ std::istream *input;
struct yy_buffer_state *buffer;
const Macro::parser::location_type yylloc;
const string for_body;
const Macro::parser::location_type for_body_loc;
- ScanContext(istream *input_arg, struct yy_buffer_state *buffer_arg,
+ ScanContext(std::istream *input_arg, struct yy_buffer_state *buffer_arg,
Macro::parser::location_type &yylloc_arg, const string &for_body_arg,
Macro::parser::location_type &for_body_loc_arg) :
input(input_arg), buffer(buffer_arg), yylloc(yylloc_arg), for_body(for_body_arg),
@@ -70,7 +70,7 @@ private:
//! Input stream used for initialization of current scanning context
/*! Kept for deletion at end of current scanning buffer */
- istream *input;
+ std::istream *input;
//! If current context is the body of a loop, contains the string of the loop body. Empty otherwise.
string for_body;
@@ -125,7 +125,7 @@ private:
and initialise a new scanning context with the loop body */
bool iter_loop(MacroDriver &driver, Macro::parser::location_type *yylloc);
public:
- MacroFlex(istream* in = 0, ostream* out = 0);
+ MacroFlex(std::istream* in = 0, ostream* out = 0);
//! The main lexing function
Macro::parser::token_type lex(Macro::parser::semantic_type *yylval,