Commit f849b115 authored by Ferhat Mihoubi's avatar Ferhat Mihoubi
Browse files

The functions erf, normpdf and normcdf work fine with bytecode option

parent 4bde0636
......@@ -242,8 +242,8 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
stack<string> Stack;
stack<double> Stackf;
ostringstream tmp_out, tmp_out2;
string v1, v2;
double v1f, v2f;
string v1, v2, v3;
double v1f, v2f, v3f;
bool go_on = true;
double ll;
ExpressionType equation_type;
......@@ -1024,10 +1024,46 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
tmp_out << "sqrt(" << v1 << ")";
Stack.push(tmp_out.str());
break;
case oErf:
Stackf.push(sqrt(v1f));
tmp_out.str("");
tmp_out << "erf(" << v1 << ")";
Stack.push(tmp_out.str());
break;
default:
;
}
break;
case FTRINARY:
op = ((FTRINARY_ *) it_code->second)->get_op_type();
v3 = Stack.top();
Stack.pop();
v2 = Stack.top();
Stack.pop();
v1 = Stack.top();
Stack.pop();
v3f = Stackf.top();
Stackf.pop();
v2f = Stackf.top();
Stackf.pop();
v1f = Stackf.top();
Stackf.pop();
switch (op)
{
case oNormcdf:
Stackf.push(0.5*(1+erf((v1f-v2f)/v3f/M_SQRT2)));
tmp_out.str("");
tmp_out << "normcdf(" << v1 << ", " << v2 << ", " << v3 << ")";
Stack.push(tmp_out.str());
break;
case oNormpdf:
Stackf.push(1/(v3f*sqrt(2*M_PI)*exp(pow((v1f-v2f)/v3f,2)/2)));
tmp_out.str("");
tmp_out << "normpdf(" << v1 << ", " << v2 << ", " << v3 << ")";
Stack.push(tmp_out.str());
break;
}
break;
case FCUML:
case FENDBLOCK:
case FOK:
......@@ -1049,7 +1085,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
int var, lag = 0, op;
ostringstream tmp_out;
double v1, v2;
double v1, v2, v3;
bool go_on = true;
double ll;
......@@ -1727,12 +1763,44 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
Stack.push(sqrt(v1));
#ifdef DEBUG
tmp_out << " |sqrt(" << v1 << ")|";
#endif
break;
case oErf:
Stack.push(erf(v1));
#ifdef DEBUG
tmp_out << " |erf(" << v1 << ")|";
#endif
break;
default:
;
}
break;
case FTRINARY:
op = ((FTRINARY_ *) it_code->second)->get_op_type();
v3 = Stack.top();
Stack.pop();
v2 = Stack.top();
Stack.pop();
v1 = Stack.top();
Stack.pop();
switch (op)
{
case oNormcdf:
//mexPrintf("normcdf(v1=%f, v2=%f, v3=%f)=%f\n", v1, v2, v3, 0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
Stack.push(0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
#ifdef DEBUG
tmp_out << " |normcdf(" << v1 << ", " << v2 << ", " << v3 << ")|";
#endif
break;
case oNormpdf:
Stack.push(1/(v3*sqrt(2*M_PI)*exp(pow((v1-v2)/v3,2)/2)));
#ifdef DEBUG
tmp_out << " |normpdf(" << v1 << ", " << v2 << ", " << v3 << ")|";
#endif
break;
}
break;
case FCUML:
v1 = Stack.top();
Stack.pop();
......
......@@ -86,6 +86,7 @@ enum Tags
FUNARY, //!< A Unary operator - 14
FBINARY, //!< A binary operator - 15
FTRINARY, //!< A trinary operator - 15'
FCUML, //!< Cumulates the result - 16
......@@ -610,6 +611,22 @@ public:
};
};
class FTRINARY_ : public TagWithOneArgument<uint8_t>
{
public:
inline FTRINARY_() : TagWithOneArgument<uint8_t>::TagWithOneArgument(FTRINARY)
{
};
inline FTRINARY_(const int op_type_arg) : TagWithOneArgument<uint8_t>::TagWithOneArgument(FTRINARY, op_type_arg)
{
};
inline uint8_t
get_op_type()
{
return arg1;
};
};
class FOK_ : public TagWithOneArgument<int>
{
public:
......@@ -1178,6 +1195,13 @@ public:
tags_liste.push_back(make_pair(FBINARY, code));
code += sizeof(FBINARY_);
break;
case FTRINARY:
# ifdef DEBUGL
mexPrintf("FTRINARY\n");
# endif
tags_liste.push_back(make_pair(FTRINARY, code));
code += sizeof(FTRINARY_);
break;
case FOK:
# ifdef DEBUGL
mexPrintf("FOK\n");
......
......@@ -237,6 +237,18 @@ ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<
return dynamic_cast<VariableNode *>(substexpr);
}
bool
ExprNode::isNumConstNodeEqualTo(double value) const
{
return false;
}
bool
ExprNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
return false;
}
NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
ExprNode(datatree_arg),
id(id_arg)
......@@ -385,6 +397,21 @@ VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg)
&& (lag == 0 || (type != eModelLocalVariable && type != eModFileLocalVariable)));
}
bool
NumConstNode::isNumConstNodeEqualTo(double value) const
{
if (datatree.num_constants.getDouble(id) == value)
return true;
else
return false;
}
bool
NumConstNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
return false;
}
void
VariableNode::prepareForDerivation()
{
......@@ -999,6 +1026,21 @@ VariableNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
return const_cast<VariableNode *>(this);
}
bool
VariableNode::isNumConstNodeEqualTo(double value) const
{
return false;
}
bool
VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
if (type == type_arg && datatree.symbol_table.getTypeSpecificID(symb_id) == variable_id && lag == lag_arg)
return true;
else
return false;
}
UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg) :
ExprNode(datatree_arg),
arg(arg_arg),
......@@ -1820,6 +1862,18 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
}
}
bool
UnaryOpNode::isNumConstNodeEqualTo(double value) const
{
return false;
}
bool
UnaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
return false;
}
BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
BinaryOpcode op_code_arg, const NodeID arg2_arg) :
ExprNode(datatree_arg),
......@@ -2810,6 +2864,18 @@ BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
}
bool
BinaryOpNode::isNumConstNodeEqualTo(double value) const
{
return false;
}
bool
BinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
return false;
}
TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) :
ExprNode(datatree_arg),
......@@ -3077,8 +3143,8 @@ TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms
arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg3->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FBINARY_ fbinary(op_code);
fbinary.write(CompileCode);
FTRINARY_ ftrinary(op_code);
ftrinary.write(CompileCode);
}
void
......@@ -3306,6 +3372,18 @@ TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOp
return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
}
bool
TrinaryOpNode::isNumConstNodeEqualTo(double value) const
{
return false;
}
bool
TrinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
return false;
}
ExternalFunctionNode::ExternalFunctionNode(DataTree &datatree_arg,
int symb_id_arg,
const vector<NodeID> &arguments_arg) :
......@@ -3632,6 +3710,18 @@ ExternalFunctionNode::getIndxInTefTerms(int the_symb_id, deriv_node_temp_terms_t
throw UnknownFunctionNameAndArgs();
}
bool
ExternalFunctionNode::isNumConstNodeEqualTo(double value) const
{
return false;
}
bool
ExternalFunctionNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
return false;
}
FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatree_arg,
int top_level_symb_id_arg,
const vector<NodeID> &arguments_arg,
......
......@@ -333,6 +333,20 @@ public:
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const = 0;
//! Return true if the nodeID is a numerical constant equal to value and false otherwise
/*!
\param[in] value of the numerical constante
\param[out] the boolean equal to true if NodeId is a constant equal to value
*/
virtual bool isNumConstNodeEqualTo(double value) const = 0;
//! Return true if the nodeID is a variable withe a type equal to type_arg, a specific variable id aqual to varfiable_id and a lag equal to lag_arg and false otherwise
/*!
\param[in] the type (type_arg), specifique variable id (variable_id and the lag (lag_arg)
\param[out] the boolean equal to true if NodeId is the variable
*/
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const = 0;
};
//! Object used to compare two nodes (using their indexes)
......@@ -378,6 +392,8 @@ public:
virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
virtual bool isNumConstNodeEqualTo(double value) const;
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
};
//! Symbol or variable node
......@@ -421,6 +437,8 @@ public:
virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
virtual bool isNumConstNodeEqualTo(double value) const;
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
};
//! Unary operator node
......@@ -482,6 +500,8 @@ public:
virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
virtual bool isNumConstNodeEqualTo(double value) const;
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
};
//! Binary operator node
......@@ -548,6 +568,8 @@ public:
virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
virtual bool isNumConstNodeEqualTo(double value) const;
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
};
//! Trinary operator node
......@@ -596,6 +618,8 @@ public:
virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
virtual bool isNumConstNodeEqualTo(double value) const;
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
};
//! External function node
......@@ -649,6 +673,8 @@ public:
virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
virtual NodeID buildSimilarExternalFunctionNode(vector<NodeID> &alt_args, DataTree &alt_datatree) const;
virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
virtual bool isNumConstNodeEqualTo(double value) const;
virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
};
class FirstDerivExternalFunctionNode : public ExternalFunctionNode
......
......@@ -414,10 +414,7 @@ t_equation_type_and_normalized_equation
ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs)
{
NodeID lhs, rhs;
ostringstream tmp_output;
BinaryOpNode *eq_node;
ostringstream tmp_s;
temporary_terms_type temporary_terms;
EquationType Equation_Simulation_Type;
t_equation_type_and_normalized_equation V_Equation_Simulation_Type(equations.size());
for (unsigned int i = 0; i < equations.size(); i++)
......@@ -429,10 +426,6 @@ ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
Equation_Simulation_Type = E_SOLVE;
tmp_s.str("");
tmp_output.str("");
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
tmp_s << "y(it_, " << Index_Var_IM[i]+1 << ")";
map<pair<int, pair<int, int> >, NodeID>::iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
pair<bool, NodeID> res;
if (derivative != first_order_endo_derivatives.end())
......@@ -441,9 +434,7 @@ ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair
derivative->second->collectEndogenous(result);
set<pair<int, int> >::const_iterator d_endo_variable = result.find(make_pair(var, 0));
//Determine whether the equation could be evaluated rather than to be solved
ostringstream tt("");
derivative->second->writeOutput(tt, oMatlabDynamicModelSparse, temporary_terms);
if (tmp_output.str() == tmp_s.str() and tt.str() == "1")
if (lhs->isVariableNodeEqualTo(eEndogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
{
Equation_Simulation_Type = E_EVALUATE;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment