Add erfc() primitive

Closes: #85

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

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)
 {

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

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

src/DynareFlex.ll

@@ -929,6 +929,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT,DYNARE_BLOCK>normcdf {return token::NORMCDF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>normpdf {return token::NORMPDF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>erf {return token::ERF;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>erfc {return token::ERFC;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>steady_state {return token::STEADY_STATE;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}

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

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)
 {

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)