v4 parser.src: added function normal cumulative density normcdf(x,mean,std dev) and its derivative

added tests/parser/t_normcdf.mod
               added Trinary Operators
               increased precision of constants written to *.m or *.c files


git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1440 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2007-11-11 15:24:50 +00:00
parent a82c653ab2
commit 2a79c2cecc
10 changed files with 468 additions and 63 deletions

View File

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

View File

@ -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 ';'

View File

@ -270,6 +270,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT,DYNARE_BLOCK>sqrt {return token::SQRT;}
<DYNARE_STATEMENT,DYNARE_BLOCK>max {return token::MAX;}
<DYNARE_STATEMENT,DYNARE_BLOCK>min {return token::MIN;}
<DYNARE_STATEMENT,DYNARE_BLOCK>normcdf {return token::NORMCDF;}
/* options for GSA module by Marco Ratto */
<DYNARE_STATEMENT>identification {return token::IDENTIFICATION;}

View File

@ -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<TrinaryOpNode *>(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<TrinaryOpNode *>(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<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
bool is_matlab) const
{
NodeID this2 = const_cast<TrinaryOpNode *>(this);
map<NodeID, int>::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<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, int> &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
map_idx_type &map_idx) const
{
NodeID this2 = const_cast<TrinaryOpNode *>(this);
map<NodeID, int>::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<TrinaryOpNode *>(this));
if (it != temporary_terms.end())
{
CompileCode.write(&FLDT, sizeof(FLDT));
int var=map_idx[idx];
CompileCode.write(reinterpret_cast<char *>(&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<char *>(&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<TrinaryOpNode *>(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<NodeID> &arguments_arg) :

View File

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

View File

@ -7,6 +7,7 @@ using namespace std;
#include <map>
#include <list>
#include <sstream>
#include <iomanip>
#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<pair<pair<NodeID, NodeID>, int>, NodeID> binary_op_node_map_type;
binary_op_node_map_type binary_op_node_map;
typedef map<pair<pair<pair<NodeID, NodeID>,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

View File

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

View File

@ -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<int> 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<NodeID, int> &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<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, int> &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
{

View File

@ -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

View File

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