diff --git a/doc/manual.xml b/doc/manual.xml
index 13264a34a244a1d3195e3d5b44a0cdd97efb71ed..994f9455c58d58e2d9cce02c007c4a9bcb57114b 100644
--- a/doc/manual.xml
+++ b/doc/manual.xml
@@ -769,6 +769,7 @@ end;
<listitem><para>maximum and minimum: <literal>max(<replaceable>a</replaceable>, <replaceable>b</replaceable>)</literal>, <literal>min(<replaceable>a</replaceable>, <replaceable>b</replaceable>)</literal></para></listitem>
<listitem><para>gaussian cumulative distribution function: <literal>normcdf(<replaceable>x</replaceable>, <replaceable>μ</replaceable>, <replaceable>σ</replaceable>)</literal> (note that <literal>normcdf(<replaceable>x</replaceable>)</literal> is equivalent to <literal>normcdf(<replaceable>x</replaceable>, 0, 1)</literal>)</para></listitem>
<listitem><para>gaussian probability density function: <literal>normpdf(<replaceable>x</replaceable>, <replaceable>μ</replaceable>, <replaceable>σ</replaceable>)</literal> (note that <literal>normpdf(<replaceable>x</replaceable>)</literal> is equivalent to <literal>normpdf(<replaceable>x</replaceable>, 0, 1)</literal>)</para></listitem>
+ <listitem><para>gauss error function: <literal>erf(<replaceable>x</replaceable>)</literal></para></listitem>
</itemizedlist>
</para>
</sect3>
diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh
index 2f81aa128b7f4d14758937f7517441a11799bedd..b542b4b52168379150ae27905c16642e1f536a0a 100644
--- a/preprocessor/CodeInterpreter.hh
+++ b/preprocessor/CodeInterpreter.hh
@@ -180,7 +180,8 @@ enum UnaryOpcode
oAtanh,
oSqrt,
oSteadyState,
- oExpectation
+ oExpectation,
+ oErf
};
enum BinaryOpcode
diff --git a/preprocessor/DataTree.cc b/preprocessor/DataTree.cc
index 27e33b89609e9c700d7e6b2fef4cc27c0830783d..6eba15a0f8dbe79fe1d2d9ebc92df5d04ad708ea 100644
--- a/preprocessor/DataTree.cc
+++ b/preprocessor/DataTree.cc
@@ -399,6 +399,15 @@ DataTree::AddSqrt(NodeID iArg1)
return Zero;
}
+NodeID
+DataTree::AddErf(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oErf, iArg1);
+ else
+ return Zero;
+}
+
NodeID
DataTree::AddMax(NodeID iArg1, NodeID iArg2)
{
diff --git a/preprocessor/DataTree.hh b/preprocessor/DataTree.hh
index 8874a00ed44752f8ca65c602befe701918b0201f..8072ba59eca2a0e97d1a42cf1dbed93fcee08bb7 100644
--- a/preprocessor/DataTree.hh
+++ b/preprocessor/DataTree.hh
@@ -174,6 +174,8 @@ public:
NodeID AddAtanh(NodeID iArg1);
//! Adds "sqrt(arg)" to model tree
NodeID AddSqrt(NodeID iArg1);
+ //! Adds "erf(arg)" to model tree
+ NodeID AddErf(NodeID iArg1);
//! Adds "max(arg1,arg2)" to model tree
NodeID AddMax(NodeID iArg1, NodeID iArg2);
//! Adds "min(arg1,arg2)" to model tree
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index c597f05c61c35757520d9884bf50291a56528215..becd2c0cb892ec2569d99e8d1b40947b15206e85 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -134,7 +134,7 @@ class ParsingDriver;
%left TIMES DIVIDE
%left UMINUS UPLUS
%nonassoc POWER
-%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN SINH COSH TANH
+%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN SINH COSH TANH ERF
%token ASINH ACOSH ATANH SQRT NORMCDF NORMPDF STEADY_STATE EXPECTATION
/* GSA analysis */
%token DYNARE_SENSITIVITY MORRIS STAB REDFORM PPRIOR PRIOR_RANGE PPOST ILPTAU GLUE MORRIS_NLIV
@@ -425,6 +425,8 @@ expression : '(' expression ')'
{ $$ = driver.add_normpdf($3, $5, $7); }
| NORMPDF '(' expression ')'
{ $$ = driver.add_normpdf($3); }
+ | ERF '(' expression ')'
+ { $$ = driver.add_erf($3); }
| NAN_CONSTANT
{ $$ = driver.add_nan_constant(); }
| INF_CONSTANT
@@ -579,6 +581,8 @@ hand_side : '(' hand_side ')'
{ $$ = driver.add_normpdf($3, $5, $7); }
| NORMPDF '(' hand_side ')'
{ $$ = driver.add_normpdf($3); }
+ | ERF '(' hand_side ')'
+ { $$ = driver.add_erf($3); }
| STEADY_STATE '(' hand_side ')'
{ $$ = driver.add_steady_state($3); }
;
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 78aeb7daab8a4c6aadf96977a3accd77c8db78f6..026a8158696e0f261810c5abfb9faf9f14b2a92f 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -467,6 +467,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT,DYNARE_BLOCK>min {return token::MIN;}
<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>steady_state {return token::STEADY_STATE;}
<DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
<DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index 0b0ffc03a7a84367136fdefc9ab84356c3360835..fd0466cd3d2ca52fd873981f873ce302612b2ee1 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -1094,6 +1094,18 @@ UnaryOpNode::composeDerivatives(NodeID darg)
return darg;
case oExpectation:
assert(0);
+ case oErf:
+ // x^2
+ t11 = datatree.AddPower(arg, datatree.Two);
+ // exp(x^2)
+ t12 = datatree.AddExp(t11);
+ // sqrt(pi)
+ t11 = datatree.AddSqrt(datatree.Pi);
+ // sqrt(pi)*exp(x^2)
+ t13 = datatree.AddTimes(t11, t12);
+ // 2/(sqrt(pi)*exp(x^2));
+ return datatree.AddDivide(datatree.Two, t13);
+ break;
}
// Suppress GCC warning
exit(EXIT_FAILURE);
@@ -1127,6 +1139,7 @@ UnaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab) c
case oLog:
return cost + 300;
case oLog10:
+ case oErf:
return cost + 16000;
case oCos:
case oSin:
@@ -1182,6 +1195,7 @@ UnaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab) c
case oCosh:
case oSinh:
case oTanh:
+ case oErf:
return cost + 240;
case oAsinh:
return cost + 220;
@@ -1356,6 +1370,9 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
return;
case oExpectation:
assert(0);
+ case oErf:
+ output << "erf";
+ break;
}
bool close_parenthesis = false;
@@ -1434,6 +1451,8 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) throw (EvalException)
return (v);
case oExpectation:
throw EvalException();
+ case oErf:
+ return (erf(v));
}
// Suppress GCC warning
exit(EXIT_FAILURE);
@@ -1539,6 +1558,8 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeI
return (make_pair(1, (NodeID) NULL));
case oExpectation:
assert(0);
+ case oErf:
+ return (make_pair(1, (NodeID) NULL));
}
}
else
@@ -1583,6 +1604,8 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeI
return (make_pair(0, datatree.AddSteadyState(New_NodeID)));
case oExpectation:
assert(0);
+ case oErf:
+ return (make_pair(0, datatree.AddErf(New_NodeID)));
}
}
return (make_pair(1, (NodeID) NULL));
@@ -1638,6 +1661,8 @@ UnaryOpNode::buildSimilarUnaryOpNode(NodeID alt_arg, DataTree &alt_datatree) con
return alt_datatree.AddSteadyState(alt_arg);
case oExpectation:
return alt_datatree.AddExpectation(expectation_information_set, alt_arg);
+ case oErf:
+ return alt_datatree.AddErf(alt_arg);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index a192147e1d7f5e418db07bba72b319a94f8d881c..3d8b63a46d34474e5dc0f434c2d958e7a5d894a8 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -1630,6 +1630,12 @@ ParsingDriver::add_normpdf(NodeID arg)
return add_normpdf(arg, data_tree->Zero, data_tree->One);
}
+NodeID
+ParsingDriver::add_erf(NodeID arg1)
+{
+ return data_tree->AddErf(arg1);
+}
+
NodeID
ParsingDriver::add_steady_state(NodeID arg1)
{
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index 6593161a603b74f7575eb624944ebb5be5cf0596..dda4553d2c9538e9e57558915b82cd1893354b18 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -481,6 +481,8 @@ public:
NodeID add_normpdf(NodeID arg1, NodeID arg2, NodeID arg3);
//! Writes token "normpdf(arg,0,1)" to model tree
NodeID add_normpdf(NodeID arg);
+ //! Writes token "erf(arg)" to model tree
+ NodeID add_erf(NodeID arg);
//! Writes token "steadyState(arg1)" to model tree
NodeID add_steady_state(NodeID arg1);
//! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)