From 6d958b6e8dfdfe71360b2afd9fb0107e87978b0b Mon Sep 17 00:00:00 2001 From: Ferhat MIHOUBI Date: Sat, 6 Feb 2010 15:07:56 +0100 Subject: [PATCH] Improvement of the floating point error messages (the equation and the location of the error are displayed) --- mex/sources/bytecode/Interpreter.cc | 1006 ++++++++++++++++++++++++-- mex/sources/bytecode/Interpreter.hh | 13 +- mex/sources/bytecode/SparseMatrix.cc | 141 ++-- 3 files changed, 1057 insertions(+), 103 deletions(-) diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index c2fd07043..df1e6c993 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -16,9 +16,6 @@ * You should have received a copy of the GNU General Public License * along with Dynare. If not, see . */ -//#define _GLIBCXX_USE_C99_FENV_TR1 1 -//#include - #include #include #include "Interpreter.hh" @@ -62,17 +59,93 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub minimal_solving_periods = minimal_solving_periods_arg; } +string +Interpreter::remove_white(string str) +{ + string temp; + for (unsigned int i = 0; i < str.length(); i++) + if (str[i] != ' ') + temp += str[i]; + return(temp); +} string -Interpreter::error_location() +Interpreter::add_underscore_to_fpe(string str) { - string tmp; + string temp; + int pos1 = -1, pos2 = -1; + string tmp_n(str.length(), ' '); + for (unsigned int i = 0; i < str.length(); i++) + { + if (str[i] != '$' && str[i] != '£') + temp += str[i]; + else + { + if (str[i] == '$') + pos1 = temp.length(); + else + pos2 = temp.length(); + if (pos1>=0 && pos2 >=0) + { + tmp_n.erase(pos1, pos2-pos1+1); + tmp_n.insert(pos1, pos2-pos1, '~'); + pos1 = pos2 = -1; + } + } + } + return(temp + "\n " + tmp_n ); +} + + +string +Interpreter::get_variable(SymbolType variable_type, int variable_num) +{ + ostringstream res; #ifndef DEBUG_EX mxArray *M_; -#endif char * P; unsigned int R, C; - stringstream Error_loc("in "); + M_ = mexGetVariable("global", "M_"); + switch(variable_type) + { + case eEndogenous: + 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++) + res << P[2*(EQN_dvar1+i*R)]; + break; + case eExogenous: + case eExogenousDet: + C = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_names"))); + R = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_names"))); + P = (char*) mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_names"))); + if (EQN_dvar1 < R) + for (unsigned int i = 0; i < C; i++) + res << P[2*(EQN_dvar1+i*R)]; + break; + case eParameter: + C = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names"))); + R = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names"))); + P = (char*) mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names"))); + if (EQN_dvar1 < R) + for (unsigned int i = 0; i < C; i++) + res << P[2*(EQN_dvar1+i*R)]; + break; + default: + break; + } +#endif + return(remove_white(res.str())); +} + + +string +Interpreter::error_location(bool evaluate) +{ + string tmp; + stringstream Error_loc(""); switch(EQN_type) { case TemporaryTerm: @@ -91,32 +164,25 @@ Interpreter::error_location() 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 + Error_loc << "first order derivative of equation " << EQN_equation+1 << " with respect to endogenous variable " << get_variable(eEndogenous, EQN_dvar1); break; default: + return("???"); break; } + Error_loc << "\n " << add_underscore_to_fpe(print_expression(it_code_expr, evaluate)); return(Error_loc.str()); } double -Interpreter::pow1(double a, double b) +Interpreter::pow1(double a, double b, bool evaluate) { double r = pow_(a, b); - if (isnan(r) || isinf(r)) + if (isnan(r)) { if (error_not_printed) { - mexPrintf("--------------------------------------\n Error: X^a with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location().c_str()); + mexPrintf("--------------------------------------\nError: X^a with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location(evaluate).c_str()); error_not_printed = false; r = 0.0000000000000000000000001; } @@ -127,14 +193,14 @@ Interpreter::pow1(double a, double b) } double -Interpreter::log1(double a) +Interpreter::log1(double a, bool evaluate) { double r = log(a); - if (isnan(r) || isinf(r)) + if (isnan(r)) { if (error_not_printed) { - mexPrintf("--------------------------------------\n Error: log(X) with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location().c_str()); + mexPrintf("--------------------------------------\nError: log(X) with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location(evaluate).c_str()); error_not_printed = false; r = 0.0000000000000000000000001; } @@ -145,14 +211,14 @@ Interpreter::log1(double a) } double -Interpreter::log10_1(double a) +Interpreter::log10_1(double a, bool evaluate) { double r = log(a); - if (isnan(r) || isinf(r)) + if (isnan(r)) { if (error_not_printed) { - mexPrintf("--------------------------------------\n Error: log10(X) with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location().c_str()); + mexPrintf("--------------------------------------\nError: log10(X) with X=%5.15f\n in %s \n--------------------------------------\n", a,error_location(evaluate).c_str()); error_not_printed = false; r = 0.0000000000000000000000001; } @@ -162,6 +228,816 @@ Interpreter::log10_1(double a) return r; } + +string +Interpreter::print_expression(it_code_type it_code, bool evaluate) +{ + int var, lag = 0, op; + stack Stack; + stack Stackf; + ostringstream tmp_out, tmp_out2; + string v1, v2; + double v1f, v2f; + bool go_on = true; + double ll; + ExpressionType equation_type; + unsigned int equation_num; + unsigned int dvar1, dvar2, dvar3; + int lag1, lag2, lag3; + size_t found; + + while (go_on) + { + switch (it_code->first) + { + case FNUMEXPR: + switch (((FNUMEXPR_ *) it_code->second)->get_expression_type()) + { + case TemporaryTerm: + equation_type = TemporaryTerm; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + break; + case ModelEquation: + equation_type = ModelEquation; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + break; + case FirstEndoDerivative: + equation_type = FirstEndoDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + break; + case FirstExoDerivative: + equation_type = FirstExoDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + break; + case FirstExodetDerivative: + equation_type = FirstExodetDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + break; + case FirstParamDerivative: + equation_type = FirstParamDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + break; + case SecondEndoDerivative: + equation_type = SecondEndoDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); + break; + case SecondExoDerivative: + equation_type = SecondExoDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); + break; + case SecondExodetDerivative: + equation_type = SecondExodetDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); + break; + case SecondParamDerivative: + equation_type = SecondExodetDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + break; + case ThirdEndoDerivative: + equation_type = ThirdEndoDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); + dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3(); + break; + case ThirdExoDerivative: + equation_type = ThirdExoDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); + dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3(); + break; + case ThirdExodetDerivative: + equation_type = ThirdExodetDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); + dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3(); + break; + case ThirdParamDerivative: + equation_type = ThirdExodetDerivative; + equation_num = ((FNUMEXPR_ *) it_code->second)->get_equation(); + dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); + dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); + break; + } + break; + case FLDV: + //load a variable in the processor + switch (((FLDV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FLDV_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << get_variable(eParameter, var); + Stack.push(tmp_out.str()); + Stackf.push(params[var]); + break; + case eEndogenous: + var = ((FLDV_ *) it_code->second)->get_pos(); + lag = ((FLDV_ *) it_code->second)->get_lead_lag(); + tmp_out.str(""); + if(lag != 0) + tmp_out << get_variable(eEndogenous, var) << "(" << lag << ")"; + else + tmp_out << get_variable(eEndogenous, var); + Stack.push(tmp_out.str()); + if (evaluate) + Stackf.push(ya[(it_+lag)*y_size+var]); + else + Stackf.push(y[(it_+lag)*y_size+var]); + break; + case eExogenous: + var = ((FLDV_ *) it_code->second)->get_pos(); + lag = ((FLDV_ *) it_code->second)->get_lead_lag(); + tmp_out.str(""); + if(lag != 0) + tmp_out << get_variable(eExogenous, var) << "(" << lag << ")"; + else + tmp_out << get_variable(eExogenous, var); + Stack.push(tmp_out.str()); + Stackf.push(x[it_+lag+var*nb_row_x]); + break; + case eExogenousDet: + var = ((FLDV_ *) it_code->second)->get_pos(); + lag = ((FLDV_ *) it_code->second)->get_lead_lag(); + tmp_out.str(""); + if(lag != 0) + tmp_out << get_variable(eExogenousDet, var) << "(" << lag << ")"; + else + tmp_out << get_variable(eExogenousDet, var); + Stack.push(tmp_out.str()); + Stackf.push(x[it_+lag+var*nb_row_xd]); + break; + case eModelLocalVariable: + break; + default: + mexPrintf("FLDV: Unknown variable type\n"); + } + break; + case FLDSV: + case FLDVS: + //load a variable in the processor + switch (((FLDSV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FLDSV_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << get_variable(eParameter, var); + Stack.push(tmp_out.str()); + Stackf.push(params[var]); + break; + case eEndogenous: + var = ((FLDSV_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << get_variable(eEndogenous, var); + Stack.push(tmp_out.str()); + if (it_code->first == FLDSV) + { + if (evaluate) + Stackf.push(ya[var]); + else + Stackf.push(y[var]); + } + else + Stackf.push(steady_y[var]); + break; + case eExogenous: + var = ((FLDSV_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << get_variable(eExogenous, var); + Stack.push(tmp_out.str()); + break; + Stackf.push(x[var]); + case eExogenousDet: + var = ((FLDSV_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << get_variable(eExogenousDet, var); + Stack.push(tmp_out.str()); + Stackf.push(x[var]); + break; + case eModelLocalVariable: + break; + default: + mexPrintf("FLDSV: Unknown variable type\n"); + } + break; + case FLDT: + //load a temporary variable in the processor + var = ((FLDT_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "T" << var+1; + Stack.push(tmp_out.str()); + Stackf.push(T[var*(periods+y_kmin+y_kmax)+it_]); + break; + case FLDST: + //load a temporary variable in the processor + var = ((FLDT_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "T" << var+1; + Stack.push(tmp_out.str()); + Stackf.push(T[var]); + break; + case FLDU: + //load u variable in the processor + var = ((FLDU_ *) it_code->second)->get_pos(); + var += Per_u_; + tmp_out.str(""); + tmp_out << "u[" << var+1 << "]"; + Stack.push(tmp_out.str()); + Stackf.push(u[var]); + break; + case FLDSU: + //load u variable in the processor + var = ((FLDSU_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "u[" << var+1 << "]"; + Stack.push(tmp_out.str()); + Stackf.push(u[var]); + break; + case FLDR: + var = ((FLDR_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "residual[" << var+1 << "]"; + Stack.push(tmp_out.str()); + Stackf.push(r[var]); + break; + case FLDZ: + //load 0 in the processor + tmp_out.str(""); + tmp_out << 0; + Stack.push(tmp_out.str()); + Stackf.push(0.0); + break; + case FLDC: + //load a numerical constant in the processor + ll = ((FLDC_ *) it_code->second)->get_value(); + tmp_out.str(""); + tmp_out << ll; + Stack.push(tmp_out.str()); + Stackf.push(ll); + break; + case FSTPV: + //load a variable in the processor + go_on = false; + switch (((FSTPV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FSTPV_ *) it_code->second)->get_pos(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eParameter, var) << " = " << tmp_out2.str(); + Stack.pop(); + params[var] = Stackf.top(); + Stackf.pop(); + break; + case eEndogenous: + var = ((FSTPV_ *) it_code->second)->get_pos(); + lag = ((FSTPV_ *) it_code->second)->get_lead_lag(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eEndogenous, var); + if (lag != 0) + tmp_out << "(" << lag << ")"; + tmp_out << " = " << tmp_out2.str(); + Stack.pop(); + y[(it_+lag)*y_size+var] = Stackf.top(); + Stackf.pop(); + break; + case eExogenous: + var = ((FSTPV_ *) it_code->second)->get_pos(); + lag = ((FSTPV_ *) it_code->second)->get_lead_lag(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eExogenous, var); + if (lag != 0) + tmp_out << "(" << lag << ")"; + tmp_out << " = " << tmp_out2.str(); + Stack.pop(); + x[it_+lag+var*nb_row_x] = Stackf.top(); + Stackf.pop(); + break; + case eExogenousDet: + var = ((FSTPV_ *) it_code->second)->get_pos(); + lag = ((FSTPV_ *) it_code->second)->get_lead_lag(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eExogenousDet, var); + if (lag != 0) + tmp_out << "(" << lag << ")"; + tmp_out << " = " << tmp_out2.str(); + Stack.pop(); + x[it_+lag+var*nb_row_xd] = Stackf.top(); + Stackf.pop(); + break; + default: + mexPrintf("FSTPV: Unknown variable type\n"); + } + break; + case FSTPSV: + go_on = false; + //load a variable in the processor + switch (((FSTPSV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FSTPSV_ *) it_code->second)->get_pos(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eParameter, var); + tmp_out << " = " << tmp_out2.str(); + Stack.pop(); + params[var] = Stackf.top(); + Stackf.pop(); + break; + case eEndogenous: + var = ((FSTPSV_ *) it_code->second)->get_pos(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eEndogenous, var); + tmp_out << " = " << tmp_out2.str(); + Stack.pop(); + y[var] = Stackf.top(); + Stackf.pop(); + break; + case eExogenous: + case eExogenousDet: + var = ((FSTPSV_ *) it_code->second)->get_pos(); + tmp_out2.str(""); + tmp_out2 << Stack.top(); + tmp_out.str(""); + tmp_out << get_variable(eExogenous, var); + tmp_out << " = " << tmp_out2.str(); + Stack.pop(); + x[var] = Stackf.top(); + Stackf.pop(); + break; + default: + mexPrintf("FSTPSV: Unknown variable type\n"); + } + break; + case FSTPT: + go_on = false; + //store in a temporary variable from the processor + var = ((FSTPT_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "T" << var+1 << " = " << Stack.top(); + Stack.pop(); + T[var*(periods+y_kmin+y_kmax)+it_] = Stackf.top(); + Stackf.pop(); + break; + case FSTPST: + go_on = false; + //store in a temporary variable from the processor + var = ((FSTPT_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "T" << var+1 << " = " << Stack.top(); + Stack.pop(); + T[var] = Stackf.top(); + Stackf.pop(); + break; + case FSTPU: + go_on = false; + //store in u variable from the processor + var = ((FSTPU_ *) it_code->second)->get_pos(); + var += Per_u_; + tmp_out.str(""); + tmp_out << "u[" << var+1 << "] = " << Stack.top(); + Stack.pop(); + u[var] = Stackf.top(); + Stackf.pop(); + break; + case FSTPSU: + go_on = false; + //store in u variable from the processor + var = ((FSTPU_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "u[" << var+1 << "] = " << Stack.top(); + Stack.pop(); + u[var] = Stackf.top(); + Stackf.pop(); + break; + case FSTPR: + go_on = false; + //store in residual variable from the processor + var = ((FSTPR_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "residual[" << var+1 << "] = " << Stack.top(); + Stack.pop(); + r[var] = Stackf.top(); + Stackf.pop(); + break; + case FSTPG: + go_on = false; + //store in derivative (g) variable from the processor + var = ((FSTPG_ *) it_code->second)->get_pos(); + tmp_out.str(""); + tmp_out << "g1[" << var+1 << "] = " << Stack.top(); + Stack.pop(); + g1[var] = Stackf.top(); + Stackf.pop(); + break; + case FBINARY: + op = ((FBINARY_ *) it_code->second)->get_op_type(); + v2 = Stack.top(); + Stack.pop(); + v1 = Stack.top(); + Stack.pop(); + v2f = Stackf.top(); + Stackf.pop(); + v1f = Stackf.top(); + Stackf.pop(); + switch (op) + { + case oPlus: + Stackf.push(v1f + v2f); + tmp_out.str(""); + tmp_out << v1 << " + " << v2; + Stack.push(tmp_out.str()); + break; + case oMinus: + Stackf.push(v1f - v2f); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " - "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oTimes: + Stackf.push(v1f * v2f); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " * "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oDivide: + double r; + r = v1f / v2f; + Stackf.push(r); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + if (isinf(r)) + tmp_out << "$"; + tmp_out << " / "; + if (isinf(r)) + tmp_out << "£"; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oLess: + Stackf.push(double (v1f < v2f)); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " < "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oGreater: + Stackf.push(double (v1f > v2f)); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " > "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oLessEqual: + Stackf.push(double (v1f <= v2f)); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " <= "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oGreaterEqual: + Stackf.push(double (v1f >= v2f)); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " >= "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oEqualEqual: + Stackf.push(double (v1f == v2f)); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " == "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oDifferent: + Stackf.push(double (v1f != v2f)); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + tmp_out << " != "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oPower: + r = pow(v1f, v2f); + Stackf.push(r); + tmp_out.str(""); + found = v1.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v1; + if (found != string::npos) + tmp_out << ")"; + if(isnan(r)) + tmp_out << "$ ^ £"; + else + tmp_out << " ^ "; + found = v2.find(" "); + if (found != string::npos) + tmp_out << "("; + tmp_out << v2; + if (found != string::npos) + tmp_out << ")"; + Stack.push(tmp_out.str()); + break; + case oMax: + Stackf.push(max(v1f, v2f)); + tmp_out.str(""); + tmp_out << "max(" << v1 << ", " << v2 << ")"; + Stack.push(tmp_out.str()); + break; + case oMin: + Stackf.push(min(v1f, v2f)); + tmp_out.str(""); + tmp_out << "min(" << v1 << ", " << v2 << ")"; + Stack.push(tmp_out.str()); + break; + case oEqual: + default: + /*throw EvalException();*/ + ; + } + break; + case FUNARY: + op = ((FUNARY_ *) it_code->second)->get_op_type(); + v1 = Stack.top(); + Stack.pop(); + v1f = Stackf.top(); + Stackf.pop(); + double r; + switch (op) + { + case oUminus: + Stackf.push(-v1f); + tmp_out.str(""); + tmp_out << " -" << v1; + Stack.push(tmp_out.str()); + break; + case oExp: + Stackf.push(exp(v1f)); + tmp_out.str(""); + tmp_out << "exp(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oLog: + r = log(v1f); + Stackf.push(r); + tmp_out.str(""); + if (isnan(r)) + tmp_out << "$log£(" << v1 << ")"; + else + tmp_out << "log(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oLog10: + r = log10(v1f); + Stackf.push(r); + tmp_out.str(""); + if (isnan(r)) + tmp_out << "$log10£(" << v1 << ")"; + else + tmp_out << "log10(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oCos: + Stackf.push(cos(v1f)); + tmp_out.str(""); + tmp_out << "cos(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oSin: + Stackf.push(sin(v1f)); + tmp_out.str(""); + tmp_out << "sin(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oTan: + Stackf.push(tan(v1f)); + tmp_out.str(""); + tmp_out << "tan(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oAcos: + Stackf.push(acos(v1f)); + tmp_out.str(""); + tmp_out << "acos(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oAsin: + Stackf.push(asin(v1f)); + tmp_out.str(""); + tmp_out << "asin(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oAtan: + Stackf.push(atan(v1f)); + tmp_out.str(""); + tmp_out << "atan(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oCosh: + Stackf.push(cosh(v1f)); + tmp_out.str(""); + tmp_out << "cosh(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oSinh: + Stackf.push(sinh(v1f)); + tmp_out.str(""); + tmp_out << "sinh(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oTanh: + Stackf.push(tanh(v1f)); + tmp_out.str(""); + tmp_out << "tanh(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oAcosh: + Stackf.push(acosh(v1f)); + tmp_out.str(""); + tmp_out << "acosh(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oAsinh: + Stackf.push(asinh(v1f)); + tmp_out.str(""); + tmp_out << "asinh(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oAtanh: + Stackf.push(atanh(v1f)); + tmp_out.str(""); + tmp_out << "atanh(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + case oSqrt: + Stackf.push(sqrt(v1f)); + tmp_out.str(""); + tmp_out << "sqrt(" << v1 << ")"; + Stack.push(tmp_out.str()); + break; + default: + ; + } + break; + case FCUML: + case FENDBLOCK: + case FOK: + case FENDEQU: + go_on = false; + break; + default: + mexPrintf("Unknown opcode %d!! FENDEQU=%d\n", it_code->first, FENDEQU); + mexErrMsgTxt("End of bytecode"); + break; + } + it_code++; + } + return(tmp_out.str()); +} + void Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) { @@ -170,13 +1046,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) double v1, v2; bool go_on = true; double ll; + EQN_block = block_num; //feclearexcept (FE_ALL_EXCEPT); - while (go_on && error_not_printed) + while (go_on) { + //tmp_it_code = it_code; switch (it_code->first) { case FNUMEXPR: + it_code_expr = it_code; switch (((FNUMEXPR_ *) it_code->second)->get_expression_type()) { case TemporaryTerm: @@ -211,7 +1090,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); break; case SecondEndoDerivative: - EQN_type = FirstEndoDerivative; + EQN_type = SecondEndoDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); @@ -219,7 +1098,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); break; case SecondExoDerivative: - EQN_type = FirstExoDerivative; + EQN_type = SecondExoDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); @@ -227,7 +1106,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); break; case SecondExodetDerivative: - EQN_type = FirstExodetDerivative; + EQN_type = SecondExodetDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); @@ -235,13 +1114,13 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2(); break; case SecondParamDerivative: - EQN_type = FirstParamDerivative; + EQN_type = SecondParamDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); break; case ThirdEndoDerivative: - EQN_type = FirstEndoDerivative; + EQN_type = ThirdEndoDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); @@ -251,7 +1130,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3(); break; case ThirdExoDerivative: - EQN_type = FirstExoDerivative; + EQN_type = ThirdExoDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); @@ -261,7 +1140,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3(); break; case ThirdExodetDerivative: - EQN_type = FirstExodetDerivative; + EQN_type = ThirdExodetDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1(); @@ -271,7 +1150,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3(); break; case ThirdParamDerivative: - EQN_type = FirstParamDerivative; + EQN_type = ThirdParamDerivative; EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation(); EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1(); EQN_dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); @@ -296,9 +1175,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) if (evaluate) Stack.push(ya[(it_+lag)*y_size+var]); else - { - Stack.push(y[(it_+lag)*y_size+var]); - } + Stack.push(y[(it_+lag)*y_size+var]); #ifdef DEBUG tmp_out << " y[" << it_+lag << ", " << var << "](" << y[(it_+lag)*y_size+var] << ")"; #endif @@ -656,11 +1533,12 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) { if (error_not_printed) { - mexPrintf("--------------------------------------\n Error: Divide by zero with %5.15f/%5.15f\n in %s \n--------------------------------------\n", v1, v2,error_location().c_str()); + mexPrintf("--------------------------------------\n Error: Divide by zero with %5.15f/%5.15f\n in %s \n--------------------------------------\n", v1, v2,error_location(evaluate).c_str()); error_not_printed = false; r = 1e70; } res1 = NAN; + go_on = false; } Stack.push(r); #ifdef DEBUG @@ -704,7 +1582,11 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) #endif break; case oPower: - Stack.push(pow1(v1, v2)); + r = pow1(v1, v2, evaluate); + if (isnan(res1)) + go_on = false; + Stack.push(r); + #ifdef DEBUG tmp_out << " |" << v1 << "^" << v2 << "|"; #endif @@ -747,13 +1629,18 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) #endif break; case oLog: - Stack.push(log1(v1)); + Stack.push(log1(v1, evaluate)); + if (isnan(res1)) + go_on = false; + #ifdef DEBUG tmp_out << " |log(" << v1 << ")|"; #endif break; case oLog10: - Stack.push(log10_1(v1)); + Stack.push(log10_1(v1, evaluate)); + if (isnan(res1)) + go_on = false; #ifdef DEBUG tmp_out << " |log10(" << v1 << ")|"; #endif @@ -1306,6 +2193,8 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, { cvg = false; iter = 0; + glambda2 = g0 = very_big; + try_at_iteration = 0; Per_y_ = it_*y_size; while (!(cvg || (iter > maxit_))) { @@ -1341,6 +2230,12 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, int prev_iter = iter; result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number); iter++; + if (iter > prev_iter) + { + g0 = res2; + gp0 = -res2; + try_at_iteration = 0; + } } if (!cvg) { @@ -1384,6 +2279,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; @@ -1412,8 +2309,15 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = false; 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); iter++; + if (iter > prev_iter) + { + g0 = res2; + gp0 = -res2; + try_at_iteration = 0; + } } if (!cvg or !result) { @@ -1447,6 +2351,8 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, { cvg = false; iter = 0; + glambda2 = g0 = very_big; + try_at_iteration = 0; Per_y_ = it_*y_size; while (!(cvg || (iter > maxit_))) { @@ -1479,8 +2385,15 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = false; if (cvg) continue; + int prev_iter = iter; result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number); iter++; + if (iter > prev_iter) + { + g0 = res2; + gp0 = -res2; + try_at_iteration = 0; + } } if (!cvg) { @@ -1530,6 +2443,8 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, if (!is_linear) { cvg = false; + glambda2 = g0 = very_big; + try_at_iteration = 0; int u_count_saved = u_count; while (!(cvg || (iter > maxit_))) { @@ -1571,8 +2486,15 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, else cvg = (max_res < solve_tolf); u_count = u_count_saved; + int prev_iter = iter; simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number); iter++; + if (iter > prev_iter) + { + g0 = res2; + gp0 = -res2; + try_at_iteration = 0; + } } if (!cvg) { diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index 90dc1a013..fe74981c8 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/mex/sources/bytecode/Interpreter.hh @@ -51,12 +51,17 @@ private: unsigned int EQN_equation, EQN_block, EQN_block_number; unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3; int EQN_lag1, EQN_lag2, EQN_lag3; + it_code_type it_code_expr; protected: - double pow1(double a, double b); - double log1(double a); - double log10_1(double a); - string error_location(); + double pow1(double a, double b, bool evaluate); + double log1(double a, bool evaluate); + double log10_1(double a, bool evaluate); + string remove_white(string str); + string add_underscore_to_fpe(string str); + string get_variable(SymbolType variable_type, int variable_num); + string error_location(bool evaluate); void compute_block_time(int Per_u_, bool evaluate, int block_num); + string print_expression(it_code_type it_code, bool evaluate); void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num, const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0); bool simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num, diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index d43821fde..17e81b424 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -1012,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) || (res2 > 8*g0 && iter>0)) + if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter>0)) { if (iter == 0) { @@ -1205,7 +1205,6 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, } double markovitz = 0, markovitz_max = -9e70; int NR_max = 0; - //mexPrintf("l=%d\n",l); if (!one) { for (j = 0; j < l; j++) @@ -1224,18 +1223,6 @@ 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 (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(1) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max); - error_not_printed = false; - markovitz = 1e70; - } - 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) { piv = piv_v[j]; @@ -1264,32 +1251,20 @@ 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 (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(2) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max); - error_not_printed = false; - markovitz = 1e70; - } - feclearexcept (FE_ALL_EXCEPT); - res1 = NAN; - }*/ - if (/*markovitz > markovitz_max &&*/ NR[j] == 1) + if (NR[j] == 1) { piv = piv_v[j]; pivj = pivj_v[j]; //Line number pivk = pivk_v[j]; //positi markovitz_max = markovitz; NR_max = NR[j]; - //mexPrintf("forced piv = %f NR_max =%d\n",piv, NR_max); } } } - if (fabs(piv) < eps) + /*if (fabs(piv) < eps) mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max); if (NR_max == 0) - mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max); + mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max);*/ pivot[i] = pivj; pivotk[i] = pivk; pivotv[i] = piv; @@ -1544,7 +1519,84 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax mexEvalString("drawnow;"); time00 = clock(); } - if (isnan(res1) || isinf(res1)) + if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter>0)) + { + if (iter == 0) + { + 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("res1=%5.10f\n", res1); + mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); + mexPrintf("Change them!\n"); + mexEvalString("drawnow;"); + filename += " stopped"; + mexErrMsgTxt(filename.c_str()); + } + if (fabs(slowc_save) < 1e-8) + { + 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;"); + filename += " stopped"; + mexErrMsgTxt(filename.c_str()); + } + 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*(periods+y_kmin); i++) + y[i] = ya[i]+slowc_save*direction[i]; + iter--; + return 0; + } + /*if (isnan(res1) || isinf(res1)) { if (iter == 0) { @@ -1577,7 +1629,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax y[i] = ya[i]+slowc_save*direction[i]; iter--; return (0); - } + }*/ u_count += u_count_init; if (alt_symbolic && alt_symbolic_count < alt_symbolic_count_max) { @@ -1699,7 +1751,6 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax } double markovitz = 0, markovitz_max = -9e70; int NR_max = 0; - //mexPrintf("l=%d\n",l); if (!one) { for (j = 0; j < l; j++) @@ -1718,18 +1769,6 @@ SparseMatrix::simulate_NG1(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 (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(1) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max); - error_not_printed = false; - markovitz = 1e70; - } - 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) { piv = piv_v[j]; @@ -1758,25 +1797,13 @@ SparseMatrix::simulate_NG1(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 (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(2) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max); - error_not_printed = false; - markovitz = 1e70; - } - feclearexcept (FE_ALL_EXCEPT); - res1 = NAN; - }*/ - if (/*markovitz > markovitz_max &&*/ NR[j] == 1) + if (NR[j] == 1) { piv = piv_v[j]; pivj = pivj_v[j]; //Line number pivk = pivk_v[j]; //positi markovitz_max = markovitz; NR_max = NR[j]; - //mexPrintf("forced piv = %f NR_max =%d\n",piv, NR_max); } } }