From d15b2110a0dcf3359eec089ff622170672946f10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 7 Dec 2021 15:19:40 +0100 Subject: [PATCH] Add erfc() primitive Closes: #85 --- src/CodeInterpreter.hh | 1 + src/DataTree.cc | 9 +++++++++ src/DataTree.hh | 2 ++ src/DynareBison.yy | 6 +++++- src/DynareFlex.ll | 1 + src/ExprNode.cc | 28 ++++++++++++++++++++++++++-- src/ParsingDriver.cc | 6 ++++++ src/ParsingDriver.hh | 2 ++ 8 files changed, 52 insertions(+), 3 deletions(-) diff --git a/src/CodeInterpreter.hh b/src/CodeInterpreter.hh index 832c6e85..ad914979 100644 --- a/src/CodeInterpreter.hh +++ b/src/CodeInterpreter.hh @@ -186,6 +186,7 @@ enum class UnaryOpcode steadyStateParam2ndDeriv, // for the 2nd derivative of the STEADY_STATE operator w.r.t. to a parameter expectation, erf, + erfc, diff, adl }; diff --git a/src/DataTree.cc b/src/DataTree.cc index f24c3030..d460ad08 100644 --- a/src/DataTree.cc +++ b/src/DataTree.cc @@ -619,6 +619,15 @@ DataTree::AddErf(expr_t iArg1) return AddUnaryOp(UnaryOpcode::erf, iArg1); } +expr_t +DataTree::AddErfc(expr_t iArg1) +{ + if (iArg1 == Zero) + return One; + + return AddUnaryOp(UnaryOpcode::erfc, iArg1); +} + expr_t DataTree::AddMax(expr_t iArg1, expr_t iArg2) { diff --git a/src/DataTree.hh b/src/DataTree.hh index 62da5681..e7baa641 100644 --- a/src/DataTree.hh +++ b/src/DataTree.hh @@ -238,6 +238,8 @@ public: expr_t AddSign(expr_t iArg1); //! Adds "erf(arg)" to model tree expr_t AddErf(expr_t iArg1); + //! Adds "erfc(arg)" to model tree + expr_t AddErfc(expr_t iArg1); //! Adds "max(arg1,arg2)" to model tree expr_t AddMax(expr_t iArg1, expr_t iArg2); //! Adds "min(arg1,arg2)" to model tree diff --git a/src/DynareBison.yy b/src/DynareBison.yy index 1e0e2782..46d14253 100644 --- a/src/DynareBison.yy +++ b/src/DynareBison.yy @@ -136,7 +136,7 @@ class ParsingDriver; %left TIMES DIVIDE %precedence UNARY %nonassoc POWER -%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN ERF DIFF ADL AUXILIARY_MODEL_NAME +%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN ERF ERFC DIFF ADL AUXILIARY_MODEL_NAME %token SQRT CBRT NORMCDF NORMPDF STEADY_STATE EXPECTATION /* GSA analysis */ %token DYNARE_SENSITIVITY MORRIS STAB REDFORM PPRIOR PRIOR_RANGE PPOST ILPTAU MORRIS_NLIV @@ -760,6 +760,8 @@ expression : '(' expression ')' { $$ = driver.add_normpdf($3); } | ERF '(' expression ')' { $$ = driver.add_erf($3); } + | ERFC '(' expression ')' + { $$ = driver.add_erfc($3); } | NAN_CONSTANT { $$ = driver.add_nan_constant(); } | INF_CONSTANT @@ -1068,6 +1070,8 @@ hand_side : '(' hand_side ')' { $$ = driver.add_normpdf($3); } | ERF '(' hand_side ')' { $$ = driver.add_erf($3); } + | ERFC '(' hand_side ')' + { $$ = driver.add_erfc($3); } | STEADY_STATE '(' hand_side ')' { $$ = driver.add_steady_state($3); } ; diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll index f9e81bb2..438dfcdc 100644 --- a/src/DynareFlex.ll +++ b/src/DynareFlex.ll @@ -929,6 +929,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]) normcdf {return token::NORMCDF;} normpdf {return token::NORMPDF;} erf {return token::ERF;} +erfc {return token::ERFC;} steady_state {return token::STEADY_STATE;} expectation {return token::EXPECTATION;} var_expectation {return token::VAR_EXPECTATION;} diff --git a/src/ExprNode.cc b/src/ExprNode.cc index ffa4166d..6f2a95af 100644 --- a/src/ExprNode.cc +++ b/src/ExprNode.cc @@ -1992,7 +1992,7 @@ UnaryOpNode::prepareForDerivation() expr_t UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id) { - expr_t t11, t12, t13, t14; + expr_t t11, t12, t13, t14, t15; switch (op_code) { @@ -2107,6 +2107,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id) cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::expectation" << endl; exit(EXIT_FAILURE); case UnaryOpcode::erf: + case UnaryOpcode::erfc: // x^2 t11 = datatree.AddPower(arg, datatree.Two); // exp(x^2) @@ -2118,7 +2119,11 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id) // 2/(sqrt(pi)*exp(x^2)); t14 = datatree.AddDivide(datatree.Two, t13); // (2/(sqrt(pi)*exp(x^2)))*dx; - return datatree.AddTimes(t14, darg); + t15 = datatree.AddTimes(t14, darg); + if (op_code == UnaryOpcode::erf) + return t15; + else // erfc + return datatree.AddUMinus(t15); case UnaryOpcode::diff: cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::diff" << endl; exit(EXIT_FAILURE); @@ -2179,6 +2184,7 @@ UnaryOpNode::cost(int cost, bool is_matlab) const return cost + 300; case UnaryOpcode::log10: case UnaryOpcode::erf: + case UnaryOpcode::erfc: return cost + 16000; case UnaryOpcode::cos: case UnaryOpcode::sin: @@ -2246,6 +2252,7 @@ UnaryOpNode::cost(int cost, bool is_matlab) const case UnaryOpcode::sinh: case UnaryOpcode::tanh: case UnaryOpcode::erf: + case UnaryOpcode::erfc: return cost + 240; case UnaryOpcode::asinh: return cost + 220; @@ -2404,6 +2411,9 @@ UnaryOpNode::writeJsonAST(ostream &output) const case UnaryOpcode::erf: output << "erf"; break; + case UnaryOpcode::erfc: + output << "erfc"; + break; } output << R"(", "arg" : )"; arg->writeJsonAST(output); @@ -2555,6 +2565,9 @@ UnaryOpNode::writeJsonOutput(ostream &output, case UnaryOpcode::erf: output << "erf"; break; + case UnaryOpcode::erfc: + output << "erfc"; + break; } bool close_parenthesis = false; @@ -2770,6 +2783,9 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, case UnaryOpcode::erf: output << "erf"; break; + case UnaryOpcode::erfc: + output << "erfc"; + break; case UnaryOpcode::diff: output << "diff"; break; @@ -2891,6 +2907,8 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) noexcept(false) exit(EXIT_FAILURE); case UnaryOpcode::erf: return erf(v); + case UnaryOpcode::erfc: + return erfc(v); case UnaryOpcode::diff: cerr << "UnaryOpNode::eval_opcode: not implemented on UnaryOpcode::diff" << endl; exit(EXIT_FAILURE); @@ -3095,6 +3113,8 @@ UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) con return alt_datatree.AddExpectation(expectation_information_set, alt_arg); case UnaryOpcode::erf: return alt_datatree.AddErf(alt_arg); + case UnaryOpcode::erfc: + return alt_datatree.AddErfc(alt_arg); case UnaryOpcode::diff: return alt_datatree.AddDiff(alt_arg); case UnaryOpcode::adl: @@ -3265,6 +3285,7 @@ UnaryOpNode::createAuxVarForUnaryOpNode() const case UnaryOpcode::abs: case UnaryOpcode::sign: case UnaryOpcode::erf: + case UnaryOpcode::erfc: return true; default: return false; @@ -3462,6 +3483,9 @@ UnaryOpNode::substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_ case UnaryOpcode::erf: unary_op = "erf"; break; + case UnaryOpcode::erfc: + unary_op = "erfc"; + break; default: cerr << "UnaryOpNode::substituteUnaryOpNodes: Shouldn't arrive here" << endl; exit(EXIT_FAILURE); diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc index 2220b8bd..23accb33 100644 --- a/src/ParsingDriver.cc +++ b/src/ParsingDriver.cc @@ -2834,6 +2834,12 @@ ParsingDriver::add_erf(expr_t arg1) return data_tree->AddErf(arg1); } +expr_t +ParsingDriver::add_erfc(expr_t arg1) +{ + return data_tree->AddErfc(arg1); +} + expr_t ParsingDriver::add_steady_state(expr_t arg1) { diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh index a4aee28b..35ba70ab 100644 --- a/src/ParsingDriver.hh +++ b/src/ParsingDriver.hh @@ -802,6 +802,8 @@ public: expr_t add_normpdf(expr_t arg); //! Writes token "erf(arg)" to model tree expr_t add_erf(expr_t arg); + //! Writes token "erfc(arg)" to model tree + expr_t add_erfc(expr_t arg); //! Writes token "steadyState(arg1)" to model tree expr_t add_steady_state(expr_t arg1); //! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)