diff --git a/parser.src/DataTree.cc b/parser.src/DataTree.cc index e903afb2b..f40d8f89f 100644 --- a/parser.src/DataTree.cc +++ b/parser.src/DataTree.cc @@ -367,6 +367,12 @@ DataTree::AddMin(NodeID iArg1, NodeID iArg2) return AddBinaryOp(iArg1, oMin, iArg2); } +NodeID +DataTree::AddNormcdf(NodeID iArg1, NodeID iArg2, NodeID iArg3) +{ + return AddTrinaryOp(iArg1, oNormcdf, iArg2, iArg3); +} + NodeID DataTree::AddEqual(NodeID iArg1, NodeID iArg2) { diff --git a/parser.src/DynareBison.yy b/parser.src/DynareBison.yy index ab90486d2..0e7f95c5b 100644 --- a/parser.src/DynareBison.yy +++ b/parser.src/DynareBison.yy @@ -70,6 +70,7 @@ class ParsingDriver; %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL %token VALUES VAR VAREXO VAREXO_DET VAROBS %token XLS_SHEET XLS_RANGE +%token NORMCDF %left LESS GREATER LESS_EQUAL GREATER_EQUAL EQUAL_EQUAL EXCLAMATION EXCLAMATION_EQUAL %left COMMA %left PLUS MINUS @@ -300,6 +301,8 @@ expression : '(' expression ')' { $$ = driver.add_min($3 , $5); } | NAME '(' comma_expression ')' { $$ = driver.add_unknown_function($1); } + | NORMCDF '(' expression COMMA expression COMMA expression ')' + { $$ = driver.add_normcdf($3,$5,$7);} ; comma_expression : expression @@ -428,6 +431,8 @@ hand_side : '(' hand_side ')' { $$ = driver.add_max($3 , $5); } | MIN '(' hand_side COMMA hand_side ')' { $$ = driver.add_min($3 , $5); } + | NORMCDF '(' expression COMMA expression COMMA expression ')' + { $$ = driver.add_normcdf($3,$5,$7);} ; pound_expression: '#' NAME EQUAL hand_side ';' diff --git a/parser.src/DynareFlex.ll b/parser.src/DynareFlex.ll index eb4908b9e..6456a749d 100644 --- a/parser.src/DynareFlex.ll +++ b/parser.src/DynareFlex.ll @@ -270,6 +270,7 @@ int sigma_e = 0; sqrt {return token::SQRT;} max {return token::MAX;} min {return token::MIN;} +normcdf {return token::NORMCDF;} /* options for GSA module by Marco Ratto */ identification {return token::IDENTIFICATION;} diff --git a/parser.src/ExprNode.cc b/parser.src/ExprNode.cc index 5e4dff284..4185c3eea 100644 --- a/parser.src/ExprNode.cc +++ b/parser.src/ExprNode.cc @@ -1301,6 +1301,275 @@ BinaryOpNode::collectEndogenous(NodeID &Id) arg2->collectEndogenous(Id); } +TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, + TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) : + ExprNode(datatree_arg), + arg1(arg1_arg), + arg2(arg2_arg), + arg3(arg3_arg), + op_code(op_code_arg) +{ + datatree.trinary_op_node_map[make_pair(make_pair(make_pair(arg1, arg2), arg3), op_code)] = this; + + // Non-null derivatives are the union of those of the arguments + // Compute set union of arg1->non_null_derivatives and arg2->non_null_derivatives + set_union(arg1->non_null_derivatives.begin(), + arg1->non_null_derivatives.end(), + arg2->non_null_derivatives.begin(), + arg2->non_null_derivatives.end(), + inserter(non_null_derivatives_tmp, non_null_derivatives_tmp.begin())); + set_union(non_null_derivatives_tmp.begin(), + non_null_derivatives_tmp.end(), + arg3->non_null_derivatives.begin(), + arg3->non_null_derivatives.end(), + inserter(non_null_derivatives, non_null_derivatives.begin())); +} + +NodeID +TrinaryOpNode::computeDerivative(int varID) +{ + NodeID darg1 = arg1->getDerivative(varID); + NodeID darg2 = arg2->getDerivative(varID); + NodeID darg3 = arg3->getDerivative(varID); + + NodeID t11, t12, t13, t14, t15; + + switch(op_code) + { + case oNormcdf: + // normal pdf is inlined in the tree + NodeID y; + t11 = datatree.AddNumConstant("2"); + t12 = datatree.AddNumConstant("3.141592653589793"); + // 2 * pi + t13 = datatree.AddTimes(t11,t12); + // sqrt(2*pi) + t14 = datatree.AddSqRt(t13); + // x - mu + t12 = datatree.AddMinus(arg1,arg2); + // y = (x-mu)/sigma + y = datatree.AddDivide(t12,arg3); + // (x-mu)^2/sigma^2 + t12 = datatree.AddTimes(y,y); + // -(x-mu)^2/sigma^2 + t13 = datatree.AddUMinus(t12); + // -((x-mu)^2/sigma^2)/2 + t12 = datatree.AddDivide(t13,t11); + // exp(-((x-mu)^2/sigma^2)/2) + t13 = datatree.AddExp(t12); + // derivative of a standardized normal + // t15 = (1/sqrt(2*pi))*exp(-y^2/2) + t15 = datatree.AddDivide(t13,t14); + // derivatives thru x + t11 = datatree.AddDivide(darg1,arg3); + // derivatives thru mu + t12 = datatree.AddDivide(darg2,arg3); + // intermediary sum + t14 = datatree.AddMinus(t11,t12); + // derivatives thru sigma + t11 = datatree.AddDivide(y,arg3); + t12 = datatree.AddTimes(t11,darg3); + //intermediary sum + t11 = datatree.AddMinus(t14,t12); + // total derivative: + // (darg1/sigma - darg2/sigma - darg3*(x-mu)/sigma)* t13 + // where t13 is the derivative of a standardized normal + return datatree.AddTimes(t11, t15); + } + cerr << "Impossible case!" << endl; + exit(-1); +} + +int +TrinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const +{ + temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); + // A temporary term behaves as a variable + if (it != temporary_terms.end()) + return 100; + + switch(op_code) + { + case oNormcdf: + return 100; + } + cerr << "Impossible case!" << endl; + exit(-1); +} + +int +TrinaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab) const +{ + temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); + // For a temporary term, the cost is null + if (it != temporary_terms.end()) + return 0; + + int cost = arg1->cost(temporary_terms, is_matlab); + cost += arg2->cost(temporary_terms, is_matlab); + + if (is_matlab) + // Cost for Matlab files + switch(op_code) + { + case oNormcdf: + return cost+1000; + } + else + // Cost for C files + switch(op_code) + { + case oNormcdf: + return cost+1000; + } + cerr << "Impossible case!" << endl; + exit(-1); +} + +void +TrinaryOpNode::computeTemporaryTerms(map &reference_count, + temporary_terms_type &temporary_terms, + bool is_matlab) const +{ + NodeID this2 = const_cast(this); + map::iterator it = reference_count.find(this2); + if (it == reference_count.end()) + { + // If this node has never been encountered, set its ref count to one, + // and travel through its children + reference_count[this2] = 1; + arg1->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + arg2->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + arg3->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + } + else + { + // If the node has already been encountered, increment its ref count + // and declare it as a temporary term if it is too costly + reference_count[this2]++; + if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)) + temporary_terms.insert(this2); + } +} + +void +TrinaryOpNode::computeTemporaryTerms(map &reference_count, + temporary_terms_type &temporary_terms, + map &first_occurence, + int Curr_block, + Model_Block *ModelBlock, + map_idx_type &map_idx) const +{ + NodeID this2 = const_cast(this); + map::iterator it = reference_count.find(this2); + if (it == reference_count.end()) + { + reference_count[this2] = 1; + first_occurence[this2] = Curr_block; + arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); + arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); + arg3->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); + } + else + { + reference_count[this2]++; + if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C) + { + temporary_terms.insert(this2); + ModelBlock->Block_List[first_occurence[this2]].Temporary_terms->insert(this2); + } + } +} + +double +TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException) +{ + switch(op_code) + { + case oNormcdf: + cerr << "NORMCDF: eval not implemented" << endl; + exit(-1); + } + cerr << "Impossible case!" << endl; + exit(-1); +} + +double +TrinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalException) +{ + double v1 = arg1->eval(eval_context); + double v2 = arg2->eval(eval_context); + double v3 = arg3->eval(eval_context); + + return eval_opcode(v1, op_code, v2, v3); +} + +/*New*/ +void +TrinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const +{ + // If current node is a temporary term + temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); + if (it != temporary_terms.end()) + { + CompileCode.write(&FLDT, sizeof(FLDT)); + int var=map_idx[idx]; + CompileCode.write(reinterpret_cast(&var), sizeof(var)); + return; + } + arg1->compile(CompileCode, lhs_rhs, output_type, temporary_terms, map_idx); + arg2->compile(CompileCode, lhs_rhs, output_type, temporary_terms, map_idx); + arg3->compile(CompileCode, lhs_rhs, output_type, temporary_terms, map_idx); + CompileCode.write(&FBINARY, sizeof(FBINARY)); + TrinaryOpcode op_codel=op_code; + CompileCode.write(reinterpret_cast(&op_codel),sizeof(op_codel)); +} +/*EndNew*/ + +void +TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, + const temporary_terms_type &temporary_terms) const +{ + + if (!OFFSET(output_type)) + { + cerr << "TrinaryOpNode not implemented for C output" << endl; + exit(-1); + } + + // If current node is a temporary term + temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); + if (it != temporary_terms.end()) + { + if (output_type != oCDynamicModelSparseDLL) + output << "T" << idx; + else + output << "T" << idx << "[it_]"; + return; + } + + switch (op_code) + { + case oNormcdf: + output << "pnorm("; + break; + } + arg1->writeOutput(output, output_type, temporary_terms); + output << ","; + arg2->writeOutput(output, output_type, temporary_terms); + output << ","; + arg3->writeOutput(output, output_type, temporary_terms); + output << ")"; +} + +void +TrinaryOpNode::collectEndogenous(NodeID &Id) +{ + arg1->collectEndogenous(Id); + arg2->collectEndogenous(Id); + arg3->collectEndogenous(Id); +} + UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg, int symb_id_arg, const vector &arguments_arg) : diff --git a/parser.src/ParsingDriver.cc b/parser.src/ParsingDriver.cc index 9085fa262..1c35332ba 100644 --- a/parser.src/ParsingDriver.cc +++ b/parser.src/ParsingDriver.cc @@ -1200,6 +1200,12 @@ ParsingDriver::add_min(NodeID arg1, NodeID arg2) return data_tree->AddMin(arg1,arg2); } +NodeID +ParsingDriver::add_normcdf(NodeID arg1, NodeID arg2, NodeID arg3) +{ + return data_tree->AddNormcdf(arg1,arg2,arg3); +} + void ParsingDriver::add_unknown_function_arg(NodeID arg) { diff --git a/parser.src/include/DataTree.hh b/parser.src/include/DataTree.hh index 8150b6790..5a2e69fe9 100644 --- a/parser.src/include/DataTree.hh +++ b/parser.src/include/DataTree.hh @@ -7,6 +7,7 @@ using namespace std; #include #include #include +#include #include "SymbolTable.hh" #include "NumericalConstants.hh" @@ -20,6 +21,7 @@ class DataTree friend class VariableNode; friend class UnaryOpNode; friend class BinaryOpNode; + friend class TrinaryOpNode; friend class UnknownFunctionNode; protected: //! A reference to the symbol table @@ -45,10 +47,13 @@ protected: unary_op_node_map_type unary_op_node_map; typedef map, int>, NodeID> binary_op_node_map_type; binary_op_node_map_type binary_op_node_map; + typedef map,NodeID>, int>, NodeID> trinary_op_node_map_type; + trinary_op_node_map_type trinary_op_node_map; inline NodeID AddPossiblyNegativeConstant(double val); inline NodeID AddUnaryOp(UnaryOpcode op_code, NodeID arg); inline NodeID AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2); + inline NodeID AddTrinaryOp(NodeID arg1, TrinaryOpcode op_code, NodeID arg2, NodeID arg3); public: DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg); virtual ~DataTree(); @@ -126,6 +131,8 @@ public: NodeID AddMaX(NodeID iArg1, NodeID iArg2); //! Adds "min(arg1,arg2)" to model tree NodeID AddMin(NodeID iArg1, NodeID iArg2); + //! Adds "normcdf(arg1,arg2,arg3)" to model tree + NodeID AddNormcdf(NodeID iArg1, NodeID iArg2, NodeID iArg3); //! Adds "arg1=arg2" to model tree NodeID AddEqual(NodeID iArg1, NodeID iArg2); void AddLocalParameter(const string &name, NodeID value) throw (LocalParameterException); @@ -144,7 +151,7 @@ DataTree::AddPossiblyNegativeConstant(double v) neg = true; } ostringstream ost; - ost << v; + ost << setprecision(16) << v; NodeID cnode = AddNumConstant(ost.str()); @@ -201,4 +208,26 @@ DataTree::AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2) return new BinaryOpNode(*this, arg1, op_code, arg2); } +inline NodeID +DataTree::AddTrinaryOp(NodeID arg1, TrinaryOpcode op_code, NodeID arg2, NodeID arg3) +{ + trinary_op_node_map_type::iterator it = trinary_op_node_map.find(make_pair(make_pair(make_pair(arg1, arg2), arg3), op_code)); + if (it != trinary_op_node_map.end()) + return it->second; + + // Try to reduce to a constant + try + { + double argval1 = arg1->eval(eval_context_type()); + double argval2 = arg2->eval(eval_context_type()); + double argval3 = arg3->eval(eval_context_type()); + double val = TrinaryOpNode::eval_opcode(argval1, op_code, argval2, argval3); + return AddPossiblyNegativeConstant(val); + } + catch(ExprNode::EvalException &e) + { + } + return new TrinaryOpNode(*this, arg1, op_code, arg2, arg3); +} + #endif diff --git a/parser.src/include/DynareBison.hh b/parser.src/include/DynareBison.hh index af014dc39..7b1503562 100644 --- a/parser.src/include/DynareBison.hh +++ b/parser.src/include/DynareBison.hh @@ -261,68 +261,69 @@ namespace yy VAROBS = 389, XLS_SHEET = 390, XLS_RANGE = 391, - EXCLAMATION_EQUAL = 392, - EXCLAMATION = 393, - EQUAL_EQUAL = 394, - GREATER_EQUAL = 395, - LESS_EQUAL = 396, - GREATER = 397, - LESS = 398, - COMMA = 399, - MINUS = 400, - PLUS = 401, - DIVIDE = 402, - TIMES = 403, - UMINUS = 404, - POWER = 405, - EXP = 406, - LOG = 407, - LOG10 = 408, - SIN = 409, - COS = 410, - TAN = 411, - ASIN = 412, - ACOS = 413, - ATAN = 414, - SINH = 415, - COSH = 416, - TANH = 417, - ASINH = 418, - ACOSH = 419, - ATANH = 420, - SQRT = 421, - DYNARE_SENSITIVITY = 422, - IDENTIFICATION = 423, - MORRIS = 424, - STAB = 425, - REDFORM = 426, - PPRIOR = 427, - PRIOR_RANGE = 428, - PPOST = 429, - ILPTAU = 430, - GLUE = 431, - MORRIS_NLIV = 432, - MORRIS_NTRA = 433, - NSAM = 434, - LOAD_REDFORM = 435, - LOAD_RMSE = 436, - LOAD_STAB = 437, - ALPHA2_STAB = 438, - KSSTAT = 439, - LOGTRANS_REDFORM = 440, - THRESHOLD_REDFORM = 441, - KSSTAT_REDFORM = 442, - ALPHA2_REDFORM = 443, - NAMENDO = 444, - NAMLAGENDO = 445, - NAMEXO = 446, - RMSE = 447, - LIK_ONLY = 448, - VAR_RMSE = 449, - PFILT_RMSE = 450, - ISTART_RMSE = 451, - ALPHA_RMSE = 452, - ALPHA2_RMSE = 453 + NORMCDF = 392, + EXCLAMATION_EQUAL = 393, + EXCLAMATION = 394, + EQUAL_EQUAL = 395, + GREATER_EQUAL = 396, + LESS_EQUAL = 397, + GREATER = 398, + LESS = 399, + COMMA = 400, + MINUS = 401, + PLUS = 402, + DIVIDE = 403, + TIMES = 404, + UMINUS = 405, + POWER = 406, + EXP = 407, + LOG = 408, + LOG10 = 409, + SIN = 410, + COS = 411, + TAN = 412, + ASIN = 413, + ACOS = 414, + ATAN = 415, + SINH = 416, + COSH = 417, + TANH = 418, + ASINH = 419, + ACOSH = 420, + ATANH = 421, + SQRT = 422, + DYNARE_SENSITIVITY = 423, + IDENTIFICATION = 424, + MORRIS = 425, + STAB = 426, + REDFORM = 427, + PPRIOR = 428, + PRIOR_RANGE = 429, + PPOST = 430, + ILPTAU = 431, + GLUE = 432, + MORRIS_NLIV = 433, + MORRIS_NTRA = 434, + NSAM = 435, + LOAD_REDFORM = 436, + LOAD_RMSE = 437, + LOAD_STAB = 438, + ALPHA2_STAB = 439, + KSSTAT = 440, + LOGTRANS_REDFORM = 441, + THRESHOLD_REDFORM = 442, + KSSTAT_REDFORM = 443, + ALPHA2_REDFORM = 444, + NAMENDO = 445, + NAMLAGENDO = 446, + NAMEXO = 447, + RMSE = 448, + LIK_ONLY = 449, + VAR_RMSE = 450, + PFILT_RMSE = 451, + ISTART_RMSE = 452, + ALPHA_RMSE = 453, + ALPHA2_RMSE = 454 }; }; diff --git a/parser.src/include/ExprNode.hh b/parser.src/include/ExprNode.hh index deddcfc59..ef913c328 100644 --- a/parser.src/include/ExprNode.hh +++ b/parser.src/include/ExprNode.hh @@ -69,6 +69,7 @@ class ExprNode friend class VariableNode; friend class UnaryOpNode; friend class BinaryOpNode; + friend class TrinaryOpNode; private: //! Computes derivative w.r. to variable varID (but doesn't store it in derivatives map) @@ -285,6 +286,41 @@ public: /*EndNew*/ }; +enum TrinaryOpcode + { + oNormcdf + }; + +//! Trinary operator node +class TrinaryOpNode : public ExprNode +{ + friend class ModelTree; +private: + const NodeID arg1, arg2, arg3; + const TrinaryOpcode op_code; + virtual NodeID computeDerivative(int varID); + virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const; + set non_null_derivatives_tmp; +public: + TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, + TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg); + virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; + virtual void computeTemporaryTerms(map &reference_count, temporary_terms_type &temporary_terms, bool is_matlab) const; + virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; + virtual void computeTemporaryTerms(map &reference_count, + temporary_terms_type &temporary_terms, + map &first_occurence, + int Curr_block, + Model_Block *ModelBlock, + map_idx_type &map_idx) const; + virtual void collectEndogenous(NodeID &Id); + static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException); + virtual double eval(const eval_context_type &eval_context) const throw (EvalException); + /*New*/ + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + /*EndNew*/ +}; + //! Unknown function node class UnknownFunctionNode : public ExprNode { diff --git a/parser.src/include/ParsingDriver.hh b/parser.src/include/ParsingDriver.hh index 30f365ce8..dfa359592 100644 --- a/parser.src/include/ParsingDriver.hh +++ b/parser.src/include/ParsingDriver.hh @@ -370,6 +370,8 @@ public: NodeID add_max(NodeID arg1, NodeID arg2); //! Writes token "min(arg1,arg2)" to model tree NodeID add_min(NodeID arg1, NodeID arg2); + //! Writes token "normcdf(arg1,arg2,arg3)" to model tree + NodeID add_normcdf(NodeID arg1, NodeID arg2, NodeID arg3); //! Adds an unknwon function argument void add_unknown_function_arg(NodeID arg); //! Adds an unknown function call node diff --git a/tests/parser/t_normcdf.mod b/tests/parser/t_normcdf.mod new file mode 100644 index 000000000..5665abc51 --- /dev/null +++ b/tests/parser/t_normcdf.mod @@ -0,0 +1,50 @@ +var y1, y2, y3, x1, x2, x3; + +model; +x1 = 1.96; +x2 = 1; +x3 = 0.5; +y1 = normcdf(x1,0,1); +y2 = normcdf(-x1,-x2,1); +y3 = normcdf(x1/2,0,x3); +end; + +initval; +y1 = 0; +y2 = 0; +y3 = 0; +x1 = 0; +x2 = 0; +x3 = 1; +end; + +steady; + +if abs(oo_.steady_state(1) - pnorm(1.96,0,1)) > 1e-12; + error('Error 1 in t_normcdf') +end; +if abs(oo_.steady_state(2) - pnorm(-1.96,-1,1)) > 1e-12; + error('Error 2 in t_normcdf') +end; +if abs(oo_.steady_state(3) - pnorm(1.96/2,0,1/2)) > 1e-12; + error('Error 3 in t_normcdf') +end; + +[junk,JJ] = t_normcdf_dynamic(oo_.steady_state,[]); + +if abs(JJ(4,4) + dnorm(1.96,0,1)) > 1e-12; + error('Error 4 in t_normcdf') +end; +if abs(JJ(5,4) - dnorm(-1.96,-1,1)) > 1e-12; + error('Error 5 in t_normcdf') +end; +if abs(JJ(5,5) + dnorm(-1.96,-1,1)) > 1e-12; + error('Error 6 in t_normcdf') +end; +if abs(JJ(6,4) + dnorm(1.96/2,0,1/2)/2) > 1e-12; + error('Error 7 in t_normcdf') +end; +if abs(JJ(6,6) - (1/2)*((1.96/2)/(1/2)^2)*dnorm(1.96/2,0,1/2)) > 1e-12; + error('Error 8 in t_normcdf') +end; +