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

Improves the global convergence of Newton method for static model

time-shift
Ferhat MIHOUBI 2010-02-05 18:33:34 +01:00
parent 4a33777ef7
commit f7ac31b58a
5 changed files with 98 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
mex/sources/bytecode/Interpreter.cc
View File

@ -62,11 +62,14 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
minimal_solving_periods = minimal_solving_periods_arg;
}
string
Interpreter::error_location()
{
string tmp;
#ifndef DEBUG_EX
mxArray *M_;
#endif
char * P;
unsigned int R, C;
stringstream Error_loc("in ");
@ -85,17 +88,19 @@ Interpreter::error_location()
Error_loc << "equation " << EQN_equation+1;
break;
case FirstEndoDerivative:
M_ = mexGetVariable("global", "M_");
if (EQN_block_number > 1)
Error_loc << "first order derivative of equation " << EQN_equation+1 << " in block " << EQN_block+1 << " with respect to endogenous variable ";
else
Error_loc << "first order derivative of equation " << EQN_equation+1 << " with respect to endogenous variable ";
#ifndef DEBUG_EX
M_ = mexGetVariable("global", "M_");
C = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "endo_names")));
R = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "endo_names")));
P = (char*) mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "endo_names")));
if (EQN_dvar1 < R)
for (unsigned int i = 0; i < C; i++)
Error_loc << P[2*(EQN_dvar1+i*R)];
#endif
break;
default:
break;
@ -145,11 +150,11 @@ double
Interpreter::log10_1(double a)
{
double r = log(a);
if (fetestexcept(FE_INVALID))
if (fetestexcept(FE_INVALID | FE_DIVBYZERO))
{
if (error_not_printed)
{
mexPrintf("--------------------------------------\n Error: log(X) with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location().c_str());
mexPrintf("--------------------------------------\n Error: log10(X) with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location().c_str());
error_not_printed = false;
r = 0.0000000000000000000000001;
}
@ -170,7 +175,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
double ll;
EQN_block = block_num;
//feclearexcept (FE_ALL_EXCEPT);
while (go_on)
while (go_on && error_not_printed)
{
switch (it_code->first)
{
@ -1237,6 +1242,8 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
max_res_idx = 0;
cvg = false;
iter = 0;
glambda2 = g0 = very_big;
try_at_iteration = 0;
while (!(cvg || (iter > maxit_)))
{
it_code = begining;
@ -1266,11 +1273,13 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, true/*false*/, cvg, iter, true, EQN_block_number);
iter++;
if (iter > prev_iter)
{
g0 = res2;
gp0 = -res2;
try_at_iteration = 0;
}
}
if (!cvg or !result)

110
mex/sources/bytecode/SparseMatrix.cc
View File

@ -933,7 +933,6 @@ SparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
NonZeroElem *first;
int i, j, k;
double yy;
res1 = res2 = max_res = 0;
for (i = 0; i < y_size*(periods+y_kmin); i++)
y[i] = ya[i];
if (symbolic && tbreak)
@ -973,7 +972,6 @@ SparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
int i, k;
double yy;
NonZeroElem *first;
res1 = res2 = max_res = 0;
for (i = 0; i < y_size; i++)
y[i+it_*y_size] = ya[i+it_*y_size];
for (i = Size-1; i >= 0; i--)
@ -1014,7 +1012,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
NR = (int *) mxMalloc(Size*sizeof(int));
error_not_printed = true;
u_count_alloc_save = u_count_alloc;
if (isnan(res1) || isinf(res1))
if (isnan(res1) || isinf(res1) || (res2 > 8*g0 && iter>0))
{
if (iter == 0)
{
@ -1051,47 +1049,61 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
if (fabs(slowc_save) < 1e-8)
{
if (slowc_save>0)
slowc_save = -slowc;
for (j = 0; j < y_size; j++)
{
bool select = false;
for (int i = 0; i < Size; i++)
if (j == index_vara[i])
{
select = true;
break;
}
if (select)
mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
else
mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
}
mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx);
mexEvalString("drawnow;");
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
if (steady_state)
return false;
else
{
for (j = 0; j < y_size; j++)
{
bool select = false;
for (int i = 0; i < Size; i++)
if (j == index_vara[i])
{
select = true;
break;
}
if (select)
mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
else
mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
}
mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx);
mexEvalString("drawnow;");
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
if (steady_state)
return false;
else
{
mexEvalString("st=fclose('all');clear all;");
filename += " stopped";
mexErrMsgTxt(filename.c_str());
}
mexEvalString("st=fclose('all');clear all;");
filename += " stopped";
mexErrMsgTxt(filename.c_str());
}
}
#ifdef GLOBAL_CONVERGENCE
if (iter == 1)
#else
slowc_save /= 2;
#endif
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)))
{
if (try_at_iteration == 0)
{
prev_slowc_save = slowc_save;
slowc_save = max( - gp0 / (2 * (res2 - g0 - gp0)) , 0.1);
}
else
{
double t1 = res2 - gp0 * slowc_save - g0;
double t2 = glambda2 - gp0 * prev_slowc_save - g0;
double a = (1/(slowc_save * slowc_save) * t1 - 1/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save);
double b = (-prev_slowc_save/(slowc_save * slowc_save) * t1 + slowc_save/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save);
prev_slowc_save = slowc_save;
slowc_save = max(min( -b + sqrt(b*b - 3 * a * gp0) / (3 * a), 0.5 * slowc_save), 0.1 * slowc_save);
}
glambda2 = res2;
try_at_iteration ++;
}
else
{
prev_slowc_save = slowc_save;
slowc_save /= 1.1;
}
if (print_it)
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
for (i = 0; i < y_size; i++)
y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size];
mxFree(piv_v);
@ -1109,6 +1121,16 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mxFree(NR);
return (true);
}
if (print_it )
{
mexPrintf("solwc=%f g0=%f res2=%f glambda2=%f\n",slowc_save,g0, res2, glambda2);
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration no %d \n", iter+1);
mexPrintf(" max. error=%.10e \n", double (max_res));
mexPrintf(" sqr. error=%.10e \n", double (res2));
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
}
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
NonZeroElem **bc;
bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc));
@ -1202,7 +1224,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
else
markovitz = fabs(piv_v[j])/piv_abs;
if (fetestexcept(FE_DIVBYZERO))
/*if (fetestexcept(FE_DIVBYZERO))
{
if (error_not_printed)
{
@ -1212,7 +1234,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
}
}*/
//mexPrintf("piv_v[j]=%f NR[j]=%d markovitz=%f markovitz_max=%f\n", piv_v[j], NR[j], markovitz, markovitz_max);
if (markovitz > markovitz_max)
{
@ -1242,7 +1264,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
else
markovitz = fabs(piv_v[j])/piv_abs;
if (fetestexcept(FE_DIVBYZERO))
/*if (fetestexcept(FE_DIVBYZERO))
{
if (error_not_printed)
{
@ -1252,7 +1274,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
}
}*/
if (/*markovitz > markovitz_max &&*/ NR[j] == 1)
{
piv = piv_v[j];

4
mex/sources/bytecode/SparseMatrix.hh
View File

@ -105,7 +105,7 @@ public:
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries);
void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u);
double g0, gp0, glambda2, try_at_iteration;
private:
void Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
@ -173,7 +173,7 @@ protected:
int u_count_alloc, u_count_alloc_save;
double *u, *y, *ya;
double res1, res2, max_res, max_res_idx;
double slowc, slowc_save, markowitz_c;
double slowc, slowc_save, prev_slowc_save, markowitz_c;
int y_kmin, y_kmax, y_size, periods, y_decal;
int *index_vara, *index_equa;
int u_count, tbreak_g;

12
mex/sources/bytecode/bytecode.cc
View File

@ -55,6 +55,14 @@ main(int argc, const char *argv[])
bool steady_state = false;
bool evaluate = false;
printf("argc=%d\n", argc);
fexcept_t *flagp;
flagp = (fexcept_t*) mxMalloc(sizeof(fexcept_t));
if (fegetexceptflag(flagp, FE_ALL_EXCEPT))
mexPrintf("fegetexceptflag failed\n");
if (fesetexceptflag(flagp,FE_INVALID | FE_DIVBYZERO))
mexPrintf("fesetexceptflag failed\n");
mxFree(flagp);
feclearexcept (FE_ALL_EXCEPT);
if (argc < 2)
{
mexPrintf("model filename expected\n");
@ -197,6 +205,10 @@ main(int argc, const char *argv[])
mxFree(ya);
if (direction)
mxFree(direction);
if (steady_yd)
mxFree(steady_yd);
if (steady_xd)
mxFree(steady_xd);
free(params);
}

4
mex/sources/bytecode/testing/simulate_debug.m
View File

@ -1,7 +1,9 @@
function simulate_debug(steady_state)
global M_ oo_ options_;
fid = fopen([M_.fname '_options.txt'],'wt');
fprintf(fid,'%d\n',options_.periods);
if steady_state~=1
fprintf(fid,'%d\n',options_.periods);
end;
fprintf(fid,'%d\n',options_.maxit_);
fprintf(fid,'%6.20f\n',options_.slowc);
fprintf(fid,'%6.20f\n',options_.markowitz);