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