diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index 399f87619..b5e290812 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -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) diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 13a53dea8..21545482a 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -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]; diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index f2cf06dc5..c4d51afb8 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -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, int>, int> &IM); void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map, 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; diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 59292040b..58f7fc33f 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -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); } diff --git a/mex/sources/bytecode/testing/simulate_debug.m b/mex/sources/bytecode/testing/simulate_debug.m index f38aeca98..a220b484f 100644 --- a/mex/sources/bytecode/testing/simulate_debug.m +++ b/mex/sources/bytecode/testing/simulate_debug.m @@ -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);