stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added normpdf as an internally supported function and updated manual.

issue#70
Houtan Bastani 2010-03-11 09:43:16 +01:00
parent 1849436c0a
commit 50258dae49
8 changed files with 87 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CodeInterpreter.hh
View File

@ -203,7 +203,8 @@ enum BinaryOpcode
enum TrinaryOpcode
{
oNormcdf
oNormcdf,
oNormpdf
};
struct Block_contain_type

6
DataTree.cc
View File

@ -417,6 +417,12 @@ DataTree::AddNormcdf(NodeID iArg1, NodeID iArg2, NodeID iArg3)
return AddTrinaryOp(iArg1, oNormcdf, iArg2, iArg3);
}
NodeID
DataTree::AddNormpdf(NodeID iArg1, NodeID iArg2, NodeID iArg3)
{
return AddTrinaryOp(iArg1, oNormpdf, iArg2, iArg3);
}
NodeID
DataTree::AddSteadyState(NodeID iArg1)
{

2
DataTree.hh
View File

@ -180,6 +180,8 @@ public:
NodeID AddMin(NodeID iArg1, NodeID iArg2);
//! Adds "normcdf(arg1,arg2,arg3)" to model tree
NodeID AddNormcdf(NodeID iArg1, NodeID iArg2, NodeID iArg3);
//! Adds "normpdf(arg1,arg2,arg3)" to model tree
NodeID AddNormpdf(NodeID iArg1, NodeID iArg2, NodeID iArg3);
//! Adds "steadyState(arg)" to model tree
NodeID AddSteadyState(NodeID iArg1);
//! Adds "arg1=arg2" to model tree

10
DynareBison.yy
View File

@ -135,7 +135,7 @@ class ParsingDriver;
%left UMINUS UPLUS
%nonassoc POWER
%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN SINH COSH TANH
%token ASINH ACOSH ATANH SQRT NORMCDF STEADY_STATE EXPECTATION
%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
%token MORRIS_NTRA NSAM LOAD_REDFORM LOAD_RMSE LOAD_STAB ALPHA2_STAB KSSTAT LOGTRANS_REDFORM THRESHOLD_REDFORM
@ -421,6 +421,10 @@ expression : '(' expression ')'
{ $$ = driver.add_normcdf($3, $5, $7); }
| NORMCDF '(' expression ')'
{ $$ = driver.add_normcdf($3); }
| NORMPDF '(' expression COMMA expression COMMA expression ')'
{ $$ = driver.add_normpdf($3, $5, $7); }
| NORMPDF '(' expression ')'
{ $$ = driver.add_normpdf($3); }
| NAN_CONSTANT
{ $$ = driver.add_nan_constant(); }
| INF_CONSTANT
@ -571,6 +575,10 @@ hand_side : '(' hand_side ')'
{ $$ = driver.add_normcdf($3, $5, $7); }
| NORMCDF '(' hand_side ')'
{ $$ = driver.add_normcdf($3); }
| NORMPDF '(' hand_side COMMA hand_side COMMA hand_side ')'
{ $$ = driver.add_normpdf($3, $5, $7); }
| NORMPDF '(' hand_side ')'
{ $$ = driver.add_normpdf($3); }
| STEADY_STATE '(' hand_side ')'
{ $$ = driver.add_steady_state($3); }
;

1
DynareFlex.ll
View File

@ -466,6 +466,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT,DYNARE_BLOCK>max {return token::MAX;}
<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>steady_state {return token::STEADY_STATE;}
<DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
<DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}

51
ExprNode.cc
View File

@ -2866,6 +2866,26 @@ TrinaryOpNode::composeDerivatives(NodeID darg1, NodeID darg2, NodeID darg3)
// (darg1/sigma - darg2/sigma - darg3*(x-mu)/sigma^2) * t15
// where t15 is the derivative of a standardized normal
return datatree.AddTimes(t11, t15);
case oNormpdf:
// (x - mu)
t11 = datatree.AddMinus(arg1, arg2);
// (x - mu)/sigma
t12 = datatree.AddDivide(t11, arg3);
// darg3 * (x - mu)/sigma
t11 = datatree.AddTimes(darg3, t12);
// darg2 - darg1
t13 = datatree.AddMinus(darg2, darg1);
// darg2 - darg1 + darg3 * (x - mu)/sigma
t14 = datatree.AddPlus(t13, t11);
// ((x - mu)/sigma) * (darg2 - darg1 + darg3 * (x - mu)/sigma)
t11 = datatree.AddTimes(t12, t14);
// ((x - mu)/sigma) * (darg2 - darg1 + darg3 * (x - mu)/sigma) - darg3
t12 = datatree.AddMinus(t11, darg3);
// this / sigma
t11 = datatree.AddDivide(this, arg3);
// total derivative:
// (this / sigma) * (((x - mu)/sigma) * (darg2 - darg1 + darg3 * (x - mu)/sigma) - darg3)
return datatree.AddTimes(t11, t12);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
@ -2891,6 +2911,7 @@ TrinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_
switch (op_code)
{
case oNormcdf:
case oNormpdf:
return 100;
}
// Suppress GCC warning
@ -2913,6 +2934,7 @@ TrinaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab)
switch (op_code)
{
case oNormcdf:
case oNormpdf:
return cost+1000;
}
else
@ -2920,6 +2942,7 @@ TrinaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab)
switch (op_code)
{
case oNormcdf:
case oNormpdf:
return cost+1000;
}
// Suppress GCC warning
@ -2988,6 +3011,8 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
{
case oNormcdf:
return (0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
case oNormpdf:
return (1/(v3*sqrt(2*M_PI)*exp(pow((v1-v2)/v3,2)/2)));
}
// Suppress GCC warning
exit(EXIT_FAILURE);
@ -3082,6 +3107,30 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
arg3->writeOutput(output, output_type, temporary_terms);
output << ")";
}
case oNormpdf:
if (IS_C(output_type))
{
//(1/(v3*sqrt(2*M_PI)*exp(pow((v1-v2)/v3,2)/2)))
output << "(1/(";
arg3->writeOutput(output, output_type, temporary_terms);
output << "*sqrt(2*M_PI)*exp(pow((";
arg1->writeOutput(output, output_type, temporary_terms);
output << "-";
arg2->writeOutput(output, output_type, temporary_terms);
output << ")/";
arg3->writeOutput(output, output_type, temporary_terms);
output << ",2)/2)))";
}
else
{
output << "normpdf(";
arg1->writeOutput(output, output_type, temporary_terms);
output << ",";
arg2->writeOutput(output, output_type, temporary_terms);
output << ",";
arg3->writeOutput(output, output_type, temporary_terms);
output << ")";
}
break;
}
}
@ -3138,6 +3187,8 @@ TrinaryOpNode::buildSimilarTrinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, NodeI
{
case oNormcdf:
return alt_datatree.AddNormcdf(alt_arg1, alt_arg2, alt_arg3);
case oNormpdf:
return alt_datatree.AddNormpdf(alt_arg1, alt_arg2, alt_arg3);
}
// Suppress GCC warning
exit(EXIT_FAILURE);

12
ParsingDriver.cc
View File

@ -1618,6 +1618,18 @@ ParsingDriver::add_normcdf(NodeID arg)
return add_normcdf(arg, data_tree->Zero, data_tree->One);
}
NodeID
ParsingDriver::add_normpdf(NodeID arg1, NodeID arg2, NodeID arg3)
{
return data_tree->AddNormpdf(arg1, arg2, arg3);
}
NodeID
ParsingDriver::add_normpdf(NodeID arg)
{
return add_normpdf(arg, data_tree->Zero, data_tree->One);
}
NodeID
ParsingDriver::add_steady_state(NodeID arg1)
{

4
ParsingDriver.hh
View File

@ -477,6 +477,10 @@ public:
NodeID add_normcdf(NodeID arg1, NodeID arg2, NodeID arg3);
//! Writes token "normcdf(arg,0,1)" to model tree
NodeID add_normcdf(NodeID arg);
//! Writes token "normpdf(arg1,arg2,arg3)" to model tree
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 "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)