stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

- Bugs corrections in deterministic simulation with sparse option 

- Check allowed with sparse option
- New command "MODEL_INFO" providing informations about the block structure of the model
- Memory leak corrections

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1993 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2008-08-25 15:06:36 +00:00
parent 4f3d871b77
commit 3bcff08e87
16 changed files with 887 additions and 468 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

162
preprocessor/BlockTriangular.cc
View File

@ -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;
cout << "alloc ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag = (" << (Lead + Lag + 1) * sizeof(IM_compact) << ")\n";
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].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++)
for(int li = 0;li < Lead + Lag + 1;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].size = tmp_size[Model_Max_Lag - Lag + li];
if(tmp_size[Model_Max_Lag - Lag + li])
{
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,24 +726,36 @@ 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(*prologue+*epilogue<n)
{
if(Do_Normalization)
{
cout << "Normalizing the model ...\n";
if(mixing)
{
double* max_val=(double*)malloc(n*sizeof(double));
//cout << "n=" << n << "\n";
memset(max_val,0,n*sizeof(double));
for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
{
//cout << "iter->first.first=" << iter->first.first << "\n";
if(fabs(iter->second)>max_val[iter->first.first])
max_val[iter->first.first]=fabs(iter->second);
}
@ -723,13 +769,17 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
{
int suppress=0;
SIM0 = (bool*)malloc(n * n * sizeof(bool));
//cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n";
memset(SIM0,0,n*n*sizeof(bool));
SIM00 = (bool*)malloc(n * n * sizeof(bool));
//cout << "Allocate SIM00=" << SIM00 << " size=" << n * n * sizeof(bool) << "\n";
memset(SIM00,0,n*n*sizeof(bool));
//cout << "n*n=" << n*n << "\n";
for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
{
if(fabs(iter->second)>bi)
{
//cout << "iter->first.first*n+iter->first.second=" << iter->first.first*n+iter->first.second << "\n";
SIM0[iter->first.first*n+iter->first.second]=1;
if(!IM_0[iter->first.first*n+iter->first.second])
{
@ -739,12 +789,21 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
else
suppress++;
}
//cout << "n*n=" << n*n << "\n";
for(i = 0;i < n;i++)
for(j = 0;j < n;j++)
{
//cout << "Index_Equ_IM[i].index * n + Index_Var_IM[j].index=" << Index_Equ_IM[i].index * n + Index_Var_IM[j].index << "\n";
SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index];
}
//cout << "free SIM0=" << SIM0 << "\n";
free(SIM0);
if(suppress!=suppressed)
{
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;
@ -764,10 +823,13 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
}
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;
}
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";

17
preprocessor/ComputingTasks.cc
View File

@ -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)
{

5
preprocessor/DynareBison.yy
View File

@ -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 ')' ';'

1
preprocessor/DynareFlex.ll
View File

@ -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;}

1
preprocessor/DynareMain.cc
View File

@ -86,6 +86,7 @@ main(int argc, char** argv)
}
// Do the rest
main2(macro_output, basename, debug, clear_all);
return 0;

4
preprocessor/DynareMain2.cc
View File

@ -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();

1
preprocessor/ModFile.cc
View File

@ -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)

3
preprocessor/ModelBlocks.cc
View File

@ -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);

27
preprocessor/ModelNormalization.cc
View File

@ -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)
{
//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

542
preprocessor/ModelTree.cc
View File

@ -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";
}
//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,19 +761,57 @@ 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;
else
{
output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl;
}
switch(ModelBlock->Block_List[j].Simulation_Type)
{
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FOREWARD_SIMPLE:
output << " g1(1)=";
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++)
{
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;
<< ") " << 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:
@ -789,8 +826,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
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
<< ", 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";
@ -809,21 +846,25 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
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 << ")";
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 << ")";
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)) = ";
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)) = ";
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";
@ -857,7 +898,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
output << " end;\n";
break;
}
}
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
}
output << "return;\n\n\n";
@ -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";
}
//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,6 +962,12 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
}
else
output << "\n\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 << " "
@ -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 << " 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:
@ -2102,11 +2139,12 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
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 (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
@ -2119,7 +2157,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
{
mStaticModelFile << " end\n";
}
mStaticModelFile << " " << static_basename << "_" << i + 1 << "(y, x);\n";
mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n";
}
open_par=false;
}
@ -2196,36 +2234,18 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
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 << " 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 << " 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";
@ -2242,16 +2262,17 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
if(open_par)
mStaticModelFile << " end;\n";
mStaticModelFile << " oo_.steady_state = y;\n";
mStaticModelFile << " if isempty(ys0_)\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 << " 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,19 +3159,16 @@ 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 << " 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 << " 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>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";
@ -3079,15 +3176,20 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
mDynamicModelFile << sp << " elseif(flag1==3)\n";
mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\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 << " 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 << " [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 << " end;\n";
mDynamicModelFile << sp << " end;\n";
mDynamicModelFile << sp << " elseif(options_.simulation_method==3),\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";
@ -3095,6 +3197,13 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons
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";
if(!block_triangular.ModelBlock->Block_List[i].is_linear)
@ -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,6 +3867,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
}
if (IM[eq*symbol_table.endo_nbr+var] && (fabs(val) < cutoff))
{
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;
}
}

9
preprocessor/ParsingDriver.cc
View File

@ -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()
{

6
preprocessor/include/BlockTriangular.hh
View File

@ -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;

10
preprocessor/include/ComputingTasks.hh
View File

@ -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:

2
preprocessor/include/ModelNormalization.hh
View File

@ -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;

7
preprocessor/include/ParsingDriver.hh
View File

@ -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

8
preprocessor/macro/MacroDriver.hh
View File

@ -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,