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

- Correction of several bugs in sparse_dll 

- Pound expressions accepted with sparse option


git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2329 ac1d8469-bf42-47a9-8791-bf33cf982152
issue#70
ferhat 2008-12-19 10:24:31 +00:00
parent 72ef44fc09
commit ff0ae482ce
12 changed files with 425 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

183
BlockTriangular.cc
View File

@ -103,12 +103,13 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
void
BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock)
{
int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1;
int *tmp_size, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_exo, tmp_nb_exo, nb_lead_lag_endo;
int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li;
int *tmp_size, *tmp_size_other_endo, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_other_endo, *tmp_exo, tmp_nb_other_endo, tmp_nb_exo, nb_lead_lag_endo;
bool *tmp_variable_evaluated;
bool *Cur_IM;
bool *IM, OK;
ModelBlock->Periods = periods;
int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo;
int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo;
if ((type == PROLOGUE) || (type == EPILOGUE))
{
for(i = 0;i < size;i++)
@ -120,15 +121,24 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
ModelBlock->Block_List[*count_Block].Temporary_terms=new temporary_terms_type ();
ModelBlock->Block_List[*count_Block].Temporary_terms->clear();
tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_other_endo = (int*)malloc(symbol_table.endo_nbr * sizeof(int));
tmp_size_other_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_var = (int*)malloc(sizeof(int));
tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
nb_lead_lag_endo = Lead = Lag = 0;
Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0;
memset(tmp_size_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_other_endo, 0, symbol_table.endo_nbr*sizeof(int));
nb_lead_lag_endo = Lead = Lag = 0;
Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = Lag_Other_Endo = Lead_Other_Endo = 0;
tmp_variable_evaluated = (bool*)malloc(symbol_table.endo_nbr*sizeof(bool));
memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool));
for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
@ -139,6 +149,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
{
if(Cur_IM[i_1 + Index_Var_IM[*count_Equ].index])
{
tmp_variable_evaluated[Index_Var_IM[*count_Equ].index] = true;
nb_lead_lag_endo++;
tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
if(k > Lead)
@ -149,6 +160,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
{
if(Cur_IM[i_1 + Index_Var_IM[*count_Equ].index])
{
tmp_variable_evaluated[Index_Var_IM[*count_Equ].index] = true;
tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
nb_lead_lag_endo++;
if(-k > Lag)
@ -160,6 +172,45 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
Lag_Endo = Lag;
Lead_Endo = Lead;
tmp_nb_other_endo = 0;
for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if(Cur_IM)
{
i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr;
for(j = 0;j < symbol_table.endo_nbr;j++)
{
int ij = Index_Var_IM[j].index;
if(Cur_IM[i_1 + ij])
{
if(!tmp_variable_evaluated[ij])
{
if(!tmp_other_endo[ij])
{
tmp_other_endo[ij] = 1;
tmp_nb_other_endo++;
}
if(k>0 && k>Lead_Other_Endo)
Lead_Other_Endo = k;
else if(k<0 && (-k)>Lag_Other_Endo)
Lag_Other_Endo = -k;
if(k>0 && k>Lead)
Lead = k;
else if(k<0 && (-k)>Lag)
Lag = -k;
tmp_size_other_endo[k+incidencematrix.Model_Max_Lag_Endo]++;
}
}
}
}
}
ModelBlock->Block_List[*count_Block].nb_other_endo = tmp_nb_other_endo;
ModelBlock->Block_List[*count_Block].Other_Endogenous = (int*)malloc(tmp_nb_other_endo * sizeof(int));
tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int));
memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int));
tmp_nb_exo = 0;
@ -214,6 +265,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(sizeof(int));
ModelBlock->Block_List[*count_Block].Equation[0] = Index_Equ_IM[*count_Equ].index;
ModelBlock->Block_List[*count_Block].Variable[0] = Index_Var_IM[*count_Equ].index;
if ((Lead > 0) && (Lag > 0))
ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
else if((Lead > 0) && (Lag == 0))
@ -267,6 +321,16 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_other_endo = tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li];
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_other_endo = tmp_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li];
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l;
IM = incidencematrix.Get_IM(li - Lag, eEndogenous);
if(IM)
@ -290,10 +354,24 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
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;
tmp_variable_evaluated[Index_Var_IM[*count_Equ].index] = true;
l++;
m++;
}
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1;
m = 0;
for(k = 0;k < symbol_table.endo_nbr;k++)
if((!tmp_variable_evaluated[Index_Var_IM[k].index]) && IM[Index_Var_IM[k].index + i_1])
{
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_other_endo[m] = l;
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_other_endo[m] = 0; //j - first_count_equ;
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_other_endo[m] = k ;
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index_other_endo[m] = Index_Equ_IM[*count_Equ].index;
ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index_other_endo[m] = Index_Var_IM[k].index;
l++;
m++;
}
}
}
else
@ -333,6 +411,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
free(tmp_endo);
free(tmp_exo);
free(tmp_var);
free(tmp_size_other_endo);
free(tmp_other_endo);
free(tmp_variable_evaluated);
}
}
else
@ -350,14 +431,21 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
first_count_equ = *count_Equ;
tmp_var = (int*)malloc(size * sizeof(int));
tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_other_endo = (int*)malloc(symbol_table.endo_nbr * sizeof(int));
tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size_other_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_size_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_other_endo, 0, symbol_table.endo_nbr*sizeof(int));
nb_lead_lag_endo = 0;
Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0;
Lag_Endo = Lead_Endo = Lag_Other_Endo = Lead_Other_Endo = Lag_Exo = Lead_Exo = 0;
//Variable by variable looking for all leads and lags its occurence in each equation of the block
tmp_variable_evaluated = (bool*)malloc(symbol_table.endo_nbr*sizeof(bool));
memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool));
for(i = 0;i < size;i++)
{
ModelBlock->Block_List[*count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index;
@ -375,6 +463,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
{
if(Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr])
{
tmp_variable_evaluated[i_1] = true;
tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
if (!OK)
{
@ -396,6 +485,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
if (!OK)
{
tmp_variable_evaluated[i_1] = true;
tmp_endo[incidencematrix.Model_Max_Lag + k]++;
nb_lead_lag_endo++;
OK = true;
@ -409,8 +499,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
}
(*count_Equ)++;
}
if ((Lag > 0) && (Lead > 0))
ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
else if(size > 1)
@ -427,12 +515,45 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
else
ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE;
}
Lag_Endo = Lag;
Lead_Endo = Lead;
tmp_nb_other_endo = 0;
for(i = 0;i < size;i++)
{
for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if(Cur_IM)
{
i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.endo_nbr;
for(j = 0;j < symbol_table.endo_nbr;j++)
if(Cur_IM[i_1 + j])
{
if(!tmp_variable_evaluated[j])
{
tmp_other_endo[j] = 1;
tmp_nb_other_endo++;
}
if(k>0 && k>Lead_Other_Endo)
Lead_Other_Endo = k;
else if(k<0 && (-k)>Lag_Other_Endo)
Lag_Other_Endo = -k;
if(k>0 && k>Lead)
Lead = k;
else if(k<0 && (-k)>Lag)
Lag = -k;
tmp_size_other_endo[k+incidencematrix.Model_Max_Lag_Endo]++;
}
}
}
}
ModelBlock->Block_List[*count_Block].nb_other_endo = tmp_nb_other_endo;
ModelBlock->Block_List[*count_Block].Other_Endogenous = (int*)malloc(tmp_nb_other_endo * sizeof(int));
tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int));
memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int));
memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int));
tmp_nb_exo = 0;
for(i = 0;i < size;i++)
{
@ -481,10 +602,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
ModelBlock->Block_List[*count_Block].Max_Lead = Lead;
ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo;
ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo;
ModelBlock->Block_List[*count_Block].Max_Lag_Other_Endo = Lag_Other_Endo;
ModelBlock->Block_List[*count_Block].Max_Lead_Other_Endo = Lead_Other_Endo;
ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo;
ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo;
ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact));
ls = l = size;
ls = l = li = size;
i1 = 0;
ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo;
for(i = 0;i < Lead + Lag + 1;i++)
@ -499,6 +622,14 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_other_endo = tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_other_endo = tmp_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
}
else
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = 0;
@ -513,6 +644,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
else
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = 0;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l;
memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool));
IM = incidencematrix.Get_IM(i - Lag, eEndogenous);
if(IM)
{
@ -541,16 +673,34 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls;
ls++;
}
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = li;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index;
tmp_variable_evaluated[Index_Var_IM[k].index] = true;
l++;
m++;
li++;
}
}
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = li - 1;
m = 0;
for(j = first_count_equ;j < size + first_count_equ;j++)
{
i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr;
for(k = 0;k < symbol_table.endo_nbr;k++)
if((!tmp_variable_evaluated[Index_Var_IM[k].index]) && IM[Index_Var_IM[k].index + i_1])
{
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_other_endo[m] = l;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_other_endo[m] = j - first_count_equ;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_other_endo[m] = k - first_count_equ;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index_other_endo[m] = Index_Equ_IM[j].index;
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index_other_endo[m] = Index_Var_IM[k].index;
l++;
m++;
}
}
ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1;
}
IM = incidencematrix.Get_IM(i - Lag, eExogenous);
if(IM)
@ -575,10 +725,13 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
}
(*count_Block)++;
free(tmp_size);
free(tmp_size_other_endo);
free(tmp_size_exo);
free(tmp_endo);
free(tmp_other_endo);
free(tmp_exo);
free(tmp_var);
free(tmp_variable_evaluated);
}
}
@ -712,7 +865,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
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";
cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << " " << iter->second << "\n";
}
}
else

3
DynareMain2.cc
View File

@ -35,6 +35,9 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool no_tm
// Run checking pass
mod_file->checkPass();
// Evaluate parameters initialization, initval, endval and pounds
mod_file->evalAllExpressions();
// Do computations
mod_file->computingPass(no_tmp_terms);

22
ExprNode.cc
View File

@ -278,7 +278,10 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case eModelLocalVariable:
case eModFileLocalVariable:
output << datatree.symbol_table.getNameByID(type, symb_id);
if(type==oMatlabDynamicModelSparse || type==oMatlabStaticModelSparse)
datatree.local_variables_table[symb_id]->writeOutput(output, output_type,temporary_terms);
else
output << datatree.symbol_table.getNameByID(type, symb_id);
break;
case eEndogenous:
@ -406,9 +409,14 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
// ModelTree::evaluateJacobian need to have the initval values applied to lead/lagged variables also
/*if (lag != 0)
throw EvalException();*/
/*if(type==eModelLocalVariable)
cout << "eModelLocalVariable = " << symb_id << "\n";*/
eval_context_type::const_iterator it = eval_context.find(make_pair(symb_id, type));
if (it == eval_context.end())
throw EvalException();
{
//cout << "unknonw variable type = " << type << " simb_id = " << symb_id << "\n";
throw EvalException();
}
return it->second;
}
@ -464,10 +472,10 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou
break;
case eModelLocalVariable:
case eModFileLocalVariable:
cerr << "VariableNode::compile: unhandled variable type" << endl;
exit(EXIT_FAILURE);
datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, output_type, temporary_terms, map_idx);
break;
case eUnknownFunction:
cerr << "Impossible case" << endl;
cerr << "Impossible case: eUnknownFuncion" << endl;
exit(EXIT_FAILURE);
}
}
@ -477,6 +485,8 @@ VariableNode::collectEndogenous(set<pair<int, int> > &result) const
{
if (type == eEndogenous)
result.insert(make_pair(symb_id, lag));
else if (type == eModelLocalVariable)
datatree.local_variables_table[symb_id]->collectEndogenous(result);
}
void
@ -484,6 +494,8 @@ VariableNode::collectExogenous(set<pair<int, int> > &result) const
{
if (type == eExogenous)
result.insert(make_pair(symb_id, lag));
else if (type == eModelLocalVariable)
datatree.local_variables_table[symb_id]->collectExogenous(result);
}

116
ModFile.cc
View File

@ -20,7 +20,7 @@
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <typeinfo>
#include "ModFile.hh"
ModFile::ModFile() : expressions_tree(symbol_table, num_constants),
@ -36,6 +36,116 @@ ModFile::~ModFile()
delete (*it);
}
void
ModFile::evalAllExpressions()
{
//Evaluate Parameters
cout << "Evaluating expressions ...";
InitParamStatement *it;
int j=0;
for(vector<Statement *>::const_iterator it1=statements.begin();it1!=statements.end(); it1++)
{
it=dynamic_cast<InitParamStatement *>(*it1);
if(it)
{
try
{
const NodeID expression = it->get_expression();
double val = expression->eval(global_eval_context);
int symb_id = symbol_table.getID(it->get_name());
global_eval_context[make_pair(symb_id, eParameter)] = val;
j++;
}
catch(ExprNode::EvalException &e)
{
cout << "error in evaluation of param\n";
}
}
}
if (j!=symbol_table.parameter_nbr)
{
cout << "Warning: Uninitialized parameters: \n";
for(j=0;j <symbol_table.parameter_nbr; j++)
{
if(global_eval_context.find(make_pair(j, eParameter))==global_eval_context.end())
cout << " " << symbol_table.getNameByID(eParameter, j) << "\n";
}
}
//Evaluate variables
for(InitOrEndValStatement::init_values_type::const_iterator it=init_values.begin(); it!=init_values.end(); it++)
{
try
{
const string &name = it->first;
const NodeID expression = it->second;
SymbolType type = symbol_table.getType(name);
double val = expression->eval(global_eval_context);
int symb_id = symbol_table.getID(name);
global_eval_context[make_pair(symb_id, type)] = val;
}
catch(ExprNode::EvalException &e)
{
cout << "error in evaluation of variable\n";
}
}
if(init_values.size()!=symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
{
cout << "\nWarning: Uninitialized variable: \n";
cout << "Endogenous\n";
for(j=0;j <symbol_table.endo_nbr; j++)
{
if(global_eval_context.find(make_pair(j, eEndogenous))==global_eval_context.end())
cout << " " << symbol_table.getNameByID(eEndogenous, j) << "\n";
}
cout << "Exogenous\n";
for(j=0;j <symbol_table.exo_nbr; j++)
{
if(global_eval_context.find(make_pair(j, eExogenous))==global_eval_context.end())
cout << " " << symbol_table.getNameByID(eExogenous, j) << "\n";
}
cout << "Deterministic exogenous\n";
for(j=0;j <symbol_table.exo_det_nbr; j++)
{
if(global_eval_context.find(make_pair(j, eExogenousDet))==global_eval_context.end())
cout << " " << symbol_table.getNameByID(eExogenousDet, j) << "\n";
}
}
//Evaluate Local variables
for(map<int, NodeID>::const_iterator it = model_tree.local_variables_table.begin(); it !=model_tree.local_variables_table.end(); it++)
{
try
{
const NodeID expression = it->second;
double val = expression->eval(global_eval_context);
//cout << it->first << " " << symbol_table.getNameByID(eModelLocalVariable, it->first) << " = " << val << "\n";
global_eval_context[make_pair(it->first, eModelLocalVariable)] = val;
}
catch(ExprNode::EvalException &e)
{
cout << "error in evaluation of pound\n";
}
}
if(model_tree.local_variables_table.size()!=symbol_table.model_local_variable_nbr+symbol_table.modfile_local_variable_nbr)
{
cout << "Warning: Unitilialized pound: \n";
cout << "Local variable in a model\n";
for(j=0;j <symbol_table.model_local_variable_nbr; j++)
{
if(global_eval_context.find(make_pair(j, eModelLocalVariable))==global_eval_context.end())
cout << " " << symbol_table.getNameByID(eModelLocalVariable, j) << "\n";
}
cout << "Local variable in a model file\n";
for(j=0;j <symbol_table.modfile_local_variable_nbr; j++)
{
if(global_eval_context.find(make_pair(j, eModFileLocalVariable))==global_eval_context.end())
cout << " " << symbol_table.getNameByID(eModFileLocalVariable, j) << "\n";
}
}
cout << "done\n";
}
void
ModFile::addStatement(Statement *st)
{
@ -112,13 +222,13 @@ ModFile::computingPass(bool no_tmp_terms)
if (mod_file_struct.order_option == 3)
model_tree.computeThirdDerivatives = true;
}
//evalAllExpressions();
model_tree.computingPass(global_eval_context, no_tmp_terms);
}
for(vector<Statement *>::iterator it = statements.begin();
it != statements.end(); it++)
(*it)->computingPass();
//evalAllExpressions();
}
void

209
ModelTree.cc
View File

@ -390,7 +390,7 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
ostringstream Uf[symbol_table.endo_nbr];
map<NodeID, int> reference_count;
int prev_Simulation_Type=-1, count_derivates=0;
int jacobian_max_endo_col;
int jacobian_max_endo_col, jacobian_max_exo_col;
ofstream output;
temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//----------------------------------------------------------------------
@ -404,7 +404,7 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
else
tmp_output << " ";
(*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
(*it)->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
}
if(tmp_output.str().length())
@ -553,8 +553,8 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
output << sps << "residual(" << i+1 << ") = (";
goto end;
/*output << sps << "residual(" << i+1 << ") = (";
goto end;*/
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = residual(" << i+1 << ")";
@ -584,8 +584,6 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl;
switch(ModelBlock->Block_List[j].Simulation_Type)
{
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD_R:
@ -619,13 +617,32 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
output << " g1(" << eq+1 << ", "
<< jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr/*ModelBlock->Block_List[j].nb_exo*/ << ") = ";
<< jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr << ") = ";
writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous);
output << "; % variable=" << symbol_table.getNameByID(eExogenous, var)
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
}
jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
if(block_triangular.incidencematrix.Model_Max_Lag_Endo - ModelBlock->Block_List[j].Max_Lag +m >=0)
{
for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_other_endo;i++)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
output << " g1(" << eq+1 << ", "
<< jacobian_max_endo_col+jacobian_max_exo_col+var+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ") = ";
writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
}
}
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE)
{
@ -645,6 +662,8 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE)
output << " end;" << endl;
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
count_derivates++;
@ -719,7 +738,6 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &
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)
@ -997,10 +1015,12 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
lhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms);
output << ";\n";
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = - residual(" << i+1 << ")";
Uf[ModelBlock->Block_List[j].Equation[i]] << "b(" << i+1 << ") = residual(" << i+1 << ")";
goto end;
default:
end:
@ -1034,33 +1054,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
output << " g1(1)=";
writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
<< "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
<< ") " << ModelBlock->Block_List[j].Variable[0]+1
<< ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
break;
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
output << " g2=0;g3=0;\n";
m=ModelBlock->Block_List[j].Max_Lag;
for(i=0;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 eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")";
output << " g1(" << eqr+1 << ", " << varr+1 << ") = ";
writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
<< "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
output << " g2=0;g3=0;\n";
@ -1075,9 +1070,9 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
if(!IM[eqr*ModelBlock->Block_List[j].Size+varr])
{
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1
<< ", " << varr+1 << ")*y( " << var+1 << ")";
IM[eqr*ModelBlock->Block_List[j].Size+varr]=1;
IM[eqr*ModelBlock->Block_List[j].Size+varr]=true;
}
output << " g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + ";
writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
@ -1090,7 +1085,6 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
#endif
}
}
output << " if(~jacobian_eval)\n";
for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
@ -1099,7 +1093,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
output << " condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
#endif
}
output << " end\n";
//output << " end\n";
#ifdef CONDITION
for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -1909,6 +1903,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
u_count_int++;
}
for(j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
int varr=block_triangular.ModelBlock->Block_List[num].Variable[j];
@ -1947,11 +1942,43 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
mStaticModelFile << "%/\n";
mStaticModelFile << "function [varargout] = " << static_basename << "(varargin)\n";
mStaticModelFile << " global oo_ M_ options_ ys0_ ;\n";
bool OK=true;
ostringstream tmp_output;
for(temporary_terms_type::const_iterator it = temporary_terms.begin();
it != temporary_terms.end(); it++)
{
if (OK)
OK=false;
else
tmp_output << " ";
(*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
}
if (tmp_output.str().length()>0)
mStaticModelFile << " global " << tmp_output.str() << " M_ ;\n";
mStaticModelFile << " T_init=0;\n";
tmp_output.str("");
for(temporary_terms_type::const_iterator it = temporary_terms.begin();
it != temporary_terms.end(); it++)
{
tmp_output << " ";
(*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
tmp_output << "=T_init;\n";
}
if (tmp_output.str().length()>0)
mStaticModelFile << tmp_output.str();
mStaticModelFile << " y_kmin=M_.maximum_lag;\n";
mStaticModelFile << " y_kmax=M_.maximum_lead;\n";
mStaticModelFile << " y_size=M_.endo_nbr;\n";
/*tmp_output.str("");
writeModelLocalVariables(tmp_output, oMatlabDynamicModel);
if (tmp_output.str().length()>0)
mStaticModelFile << tmp_output.str() << "\n";*/
mStaticModelFile << " if(length(varargin)>0)\n";
mStaticModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
mStaticModelFile << " %A simple evaluation of the static model\n";
//mStaticModelFile << " global it_;\n";
mStaticModelFile << " y=varargin{1}(:);\n";
mStaticModelFile << " ys=y;\n";
@ -1994,21 +2021,6 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
}
}
/*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 << " y_index_eq=[";
for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
}
mStaticModelFile << " ];\n";
k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;*/
if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))
skip_head=true;
@ -2030,8 +2042,6 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
}
else
ga_index++;
/*mStaticModelFile << " residual(y_index)=ys(y_index)-y(y_index);\n";
mStaticModelFile << " g1(y_index_eq, y_index) = ga(" << ga_index << ", " << ga_index << ");\n";*/
break;
case SOLVE_FORWARD_COMPLETE:
case SOLVE_BACKWARD_COMPLETE:
@ -2055,7 +2065,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
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 << " %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";
@ -2088,32 +2098,6 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
}
open_par=false;
}
/*else if ((k == SOLVE_FORWARD_SIMPLE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
{
if (open_par)
{
mStaticModelFile << " end\n";
}
open_par=false;
mStaticModelFile << " g1=0;\n";
mStaticModelFile << " r=0;\n";
mStaticModelFile << " cvg=0;\n";
mStaticModelFile << " iter=0;\n";
mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n";
mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n";
mStaticModelFile << " y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
mStaticModelFile << " cvg=((r*r)<solve_tolf);\n";
mStaticModelFile << " iter=iter+1;\n";
mStaticModelFile << " end\n";
mStaticModelFile << " if cvg==0\n";
mStaticModelFile << " if(options_.cutoff == 0)\n";
mStaticModelFile << " fprintf('Error in steady: Convergence not achieved in block " << i+1 << ", after %d iterations.\\n Increase \"options_.maxit_\".\\n',iter);\n";
mStaticModelFile << " else\n";
mStaticModelFile << " fprintf('Error in steady: Convergence not achieved in block " << i+1 << ", after %d iterations.\\n Increase \"options_.maxit_\" or set \"cutoff=0\" in model options.\\n',iter);\n";
mStaticModelFile << " end;\n";
mStaticModelFile << " return;\n";
mStaticModelFile << " end\n";
}*/
else if ((k == SOLVE_FORWARD_SIMPLE || k == SOLVE_BACKWARD_SIMPLE) || (k == SOLVE_FORWARD_COMPLETE || k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_TWO_BOUNDARIES_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
{
if (open_par)
@ -2134,46 +2118,11 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
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;
/*mStaticModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n";
mStaticModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n";
mStaticModelFile << " else\n";
mStaticModelFile << " blck_num = 1;\n";
mStaticModelFile << " end;\n";*/
mStaticModelFile << " y = solve_one_boundary('" << static_basename << "_" << i + 1 << "'" <<
", y, x, params, y_index, " << nze <<
", y_kmin, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", 1, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", " << Blck_Num << ", y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 0, 0);\n";
/*mStaticModelFile << " r=0;\n";
mStaticModelFile << " cvg=0;\n";
mStaticModelFile << " iter=0;\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, 0);\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, 0);\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 << " if(options_.cutoff == 0)\n";
mStaticModelFile << " fprintf('Error in steady: Convergence not achieved in block " << i+1 << ", after %d iterations.\\n Increase \"options_.maxit_\".\\n',iter);\n";
mStaticModelFile << " else\n";
mStaticModelFile << " fprintf('Error in steady: Convergence not achieved in block " << i+1 << ", after %d iterations.\\n Increase \"options_.maxit_\" or set \"cutoff=0\" in model options.\\n',iter);\n";
mStaticModelFile << " end;\n";
mStaticModelFile << " return;\n";
mStaticModelFile << " else\n";
mStaticModelFile << " fprintf('convergence achieved after %d iterations\\n',iter);\n";
mStaticModelFile << " end\n";*/
}
prev_Simulation_Type=k;
}
@ -2449,10 +2398,8 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n";
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n";
mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
//mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, g1, g2, g3, y_kmin, periods);\n";
}
//open_par=true;
}
else if ((k == SOLVE_FORWARD_COMPLETE || k == SOLVE_FORWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
{
@ -2461,12 +2408,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
open_par=false;
mDynamicModelFile << " g1=0;\n";
mDynamicModelFile << " r=0;\n";
tmp_eq.str("");
tmp.str("");
for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
mDynamicModelFile << " y_index = [" << tmp.str() << "];\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;
@ -2476,7 +2423,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
mDynamicModelFile << " blck_num = 1;\n";
mDynamicModelFile << " end;\n";
mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" <<
", y, x, params, y_index_eq, " << nze <<
", y, x, params, y_index, " << nze <<
", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n";
@ -2488,12 +2435,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
open_par=false;
mDynamicModelFile << " g1=0;\n";
mDynamicModelFile << " r=0;\n";
tmp_eq.str("");
tmp.str("");
for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
mDynamicModelFile << " y_index = [" << tmp.str() << "];\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;
@ -2503,7 +2450,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
mDynamicModelFile << " blck_num = 1;\n";
mDynamicModelFile << " end;\n";
mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" <<
", y, x, params, y_index_eq, " << nze <<
", y, x, params, y_index, " << nze <<
", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n";
}
@ -3072,7 +3019,9 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
}
catch(ExprNode::EvalException &e)
{
cerr << "ModelTree::evaluateJacobian: evaluation of Jacobian failed!" << endl;
cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getNameByID(eEndogenous, variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ") [" << variable_table.getSymbolID(it->first.second) << "] !" << endl;
Id->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);cout << "\n";
cerr << "ModelTree::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getNameByID(eEndogenous, variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ")!" << endl;
}
int eq=it->first.first;
int var=variable_table.getSymbolID(it->first.second);

13
NumericalInitialization.cc
View File

@ -38,6 +38,19 @@ InitParamStatement::writeOutput(ostream &output, const string &basename) const
output << param_name << " = M_.params( " << id << " );\n";
}
NodeID
InitParamStatement::get_expression() const
{
return(param_value);
}
string
InitParamStatement::get_name() const
{
return(param_name);
}
InitOrEndValStatement::InitOrEndValStatement(const init_values_type &init_values_arg,
const SymbolTable &symbol_table_arg) :
init_values(init_values_arg),

36
ParsingDriver.cc
View File

@ -74,8 +74,8 @@ ParsingDriver::parse(istream &in, bool debug)
/* Deleting filename in DynareFlex::yyterminate() is prematurate,
because if there is an unexpected end of file, the call to
ParsingDriver::error() needs filename */
if (location.begin.filename)
delete location.begin.filename;
/*if (location.begin.filename)
delete location.begin.filename;*/
return mod_file;
}
@ -179,13 +179,6 @@ ParsingDriver::add_model_variable(string *name, string *olag)
NodeID id = model_tree->AddVariable(*name, lag);
/*if (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode)
{
if (type == eEndogenous)
model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
if (type == eExogenous)
model_tree->block_triangular.fill_IM_X(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
}*/
delete name;
delete olag;
return id;
@ -250,6 +243,7 @@ ParsingDriver::dsample(string *arg1, string *arg2)
delete arg2;
}
void
ParsingDriver::init_param(string *name, NodeID rhs)
{
@ -260,7 +254,7 @@ ParsingDriver::init_param(string *name, NodeID rhs)
mod_file->addStatement(new InitParamStatement(*name, rhs, mod_file->symbol_table));
// Update global eval context
try
/*try
{
double val = rhs->eval(mod_file->global_eval_context);
int symb_id = mod_file->symbol_table.getID(*name);
@ -269,7 +263,7 @@ ParsingDriver::init_param(string *name, NodeID rhs)
catch(ExprNode::EvalException &e)
{
}
*/
delete name;
}
@ -283,11 +277,11 @@ ParsingDriver::init_val(string *name, NodeID rhs)
&& type != eExogenous
&& type != eExogenousDet)
error("initval/endval: " + *name + " should be an endogenous or exogenous variable");
init_values.push_back(make_pair(*name, rhs));
//cout << "mod_file->init_values = " << mod_file->init_values << "\n";
mod_file->init_values.push_back(make_pair(*name, rhs));
//cout << "init_val " << *name << " mod_file->init_values.size()=" << mod_file->init_values.size() << "\n";
// Update global evaluation context
try
/*try
{
double val = rhs->eval(mod_file->global_eval_context);
int symb_id = mod_file->symbol_table.getID(*name);
@ -296,7 +290,7 @@ ParsingDriver::init_val(string *name, NodeID rhs)
catch(ExprNode::EvalException &e)
{
}
*/
delete name;
}
@ -379,15 +373,17 @@ ParsingDriver::sparse()
void
ParsingDriver::end_initval()
{
mod_file->addStatement(new InitValStatement(init_values, mod_file->symbol_table));
init_values.clear();
mod_file->addStatement(new InitValStatement(mod_file->init_values, mod_file->symbol_table));
//mod_file->init_values.clear();
//cout << "mod_file->init_values.clear() in end_initval()\n";
}
void
ParsingDriver::end_endval()
{
mod_file->addStatement(new EndValStatement(init_values, mod_file->symbol_table));
init_values.clear();
mod_file->addStatement(new EndValStatement(mod_file->init_values, mod_file->symbol_table));
//mod_file->init_values.clear();
//cout << "mod_file->init_values.clear() in end_endval()\n";
}
void

4
include/DataTree.hh
View File

@ -56,8 +56,6 @@ protected:
//! A counter for filling ExprNode's idx field
int node_counter;
//! Stores local variables value
map<int, NodeID> local_variables_table;
typedef map<int, NodeID> num_const_node_map_type;
num_const_node_map_type num_const_node_map;
@ -81,6 +79,8 @@ public:
//! The variable table
VariableTable variable_table;
NodeID Zero, One, MinusOne;
//! Stores local variables value
map<int, NodeID> local_variables_table;
//! Raised when a local parameter is declared twice
class LocalParameterException

10
include/ExprNode.hh
View File

@ -321,27 +321,31 @@ public:
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
struct IM_compact
{
int size, u_init, u_finish, nb_endo, size_exo;
int size, u_init, u_finish, nb_endo, nb_other_endo, size_exo, size_other_endo;
int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index;
int *u_other_endo, *Var_other_endo, *Equ_other_endo, *Var_Index_other_endo, *Equ_Index_other_endo;
};
//! One block of the model
struct Block
{
int Size, Sized, nb_exo, nb_exo_det;
int Size, Sized, nb_exo, nb_exo_det, nb_other_endo;
BlockType Type;
BlockSimulationType Simulation_Type;
int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo;
int Max_Lag_Endo, Max_Lead_Endo;
int Max_Lag_Other_Endo, Max_Lead_Other_Endo;
int Max_Lag_Exo, Max_Lead_Exo;
bool is_linear;
int *Equation, *Own_Derivative;
int *Variable, *Exogenous;
int *Variable, *Other_Endogenous, *Exogenous;
temporary_terms_type *Temporary_terms;
IM_compact *IM_lead_lag;
int Code_Start, Code_Length;
};
//! The set of all blocks of the model
struct Model_Block
{

5
include/ModFile.hh
View File

@ -26,6 +26,7 @@ using namespace std;
#include "SymbolTable.hh"
#include "NumericalConstants.hh"
#include "NumericalInitialization.hh"
#include "ModelTree.hh"
#include "VariableTable.hh"
#include "Statement.hh"
@ -49,6 +50,8 @@ public:
//! Global evaluation context
/*! Filled using initval blocks and parameters initializations */
eval_context_type global_eval_context;
//! Temporary storage for initval/endval blocks
InitOrEndValStatement::init_values_type init_values;
private:
//! List of statements
@ -59,6 +62,8 @@ private:
public:
//! Add a statement
void addStatement(Statement *st);
//! Evaluate all the statements
void evalAllExpressions();
//! Do some checking and fills mod_file_struct
/*! \todo add check for number of equations and endogenous if ramsey_policy is present */
void checkPass();

2
include/NumericalInitialization.hh
View File

@ -39,6 +39,8 @@ public:
InitParamStatement(const string &param_name_arg, const NodeID param_value_arg,
const SymbolTable &symbol_table_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
NodeID get_expression() const;
string get_name() const;
};
class InitOrEndValStatement : public Statement

2
include/ParsingDriver.hh
View File

@ -125,8 +125,6 @@ private:
SigmaeStatement::row_type sigmae_row;
//! Temporary storage for Sigma_e matrix
SigmaeStatement::matrix_type sigmae_matrix;
//! Temporary storage for initval/endval blocks
InitOrEndValStatement::init_values_type init_values;
//! Temporary storage for histval blocks
HistValStatement::hist_values_type hist_values;
//! Temporary storage for homotopy_setup blocks