Skip to content

GitLab

Commit 35450a29 authored by Houtan Bastani's avatar Houtan Bastani
Browse files

extend domain of derivative of x^p to x=0 under certain conditions (see ticket 78)

parent d9b9360a
Hide whitespace changes
Inline Side-by-side
preprocessor/CodeInterpreter.hh
View file @ 35450a29
/*
* Copyright (C) 2007-2010 Dynare Team
*
......@@ -200,6 +199,7 @@ enum BinaryOpcode
oTimes,
oDivide,
oPower,
oPowerDeriv, // for the derivative of the power function (see trac ticket #78)
oEqual,
oMax,
oMin,
......
preprocessor/DataTree.cc
View file @ 35450a29
......@@ -240,6 +240,13 @@ DataTree::AddPower(expr_t iArg1, expr_t iArg2)
return Zero;
}
expr_t
DataTree::AddPowerDeriv(expr_t iArg1, expr_t iArg2, int powerDerivOrder)
{
assert(powerDerivOrder > 0);
return AddBinaryOp(iArg1, oPowerDeriv, iArg2, powerDerivOrder);
}
expr_t
DataTree::AddExp(expr_t iArg1)
{
......@@ -607,3 +614,50 @@ DataTree::minLagForSymbol(int symb_id) const
r = it->first.second;
return r;
}
void
DataTree::writePowerDerivCHeader(ostream &output) const
{
if (isBinaryOpUsed(oPowerDeriv))
output << "double getPowerDeriv(double, double, int);" << endl;
}
void
DataTree::writePowerDeriv(ostream &output, bool use_dll) const
{
if (!isBinaryOpUsed(oPowerDeriv))
return;
if (use_dll)
output << "/*" << endl
<< " * The k-th derivative of x^p" << endl
<< " */" << endl
<< "double getPowerDeriv(double x, double p, int k)" << endl
<< "{" << endl
<< " if ( fabs(x) < " << NEAR_ZERO << " && p > 0 && k >= p && fabs(p-nearbyint(p)) < " << NEAR_ZERO << " )" << endl
<< " return 0.0;" << endl
<< " else" << endl
<< " {" << endl
<< " int i = 0;" << endl
<< " double dxp = pow(x, p-k);" << endl
<< " for (; i<k; i++)" << endl
<< " dxp *= p--;" << endl
<< " return dxp;" << endl
<< " }" << endl
<< "}" << endl;
else
output << endl
<< "%" << endl
<< "% The k-th derivative of x^p" << endl
<< "%" << endl
<< "function dxp=getPowerDeriv(x,p,k)" << endl
<< "if (abs(x) < " << NEAR_ZERO << ") && (p > 0) && (k >= p) && (abs(p - round(p)) < " << NEAR_ZERO << ")" << endl
<< " dxp = 0;" << endl
<< "else" << endl
<< " dxp = x^(p-k);" << endl
<< " for i=0:k-1" << endl
<< " dxp = dxp*p;" << endl
<< " p = p-1;" << endl
<< " end" << endl
<< "end" << endl;
}
preprocessor/DataTree.hh
View file @ 35450a29
......@@ -62,7 +62,8 @@ protected:
//! Pair( Pair (arg1, UnaryOpCode), Pair(Expectation Info Set, Expectation Info Set Name) )
typedef map<pair<pair<expr_t, UnaryOpcode>, pair<int, string> >, UnaryOpNode *> unary_op_node_map_t;
unary_op_node_map_t unary_op_node_map;
typedef map<pair<pair<expr_t, expr_t>, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t;
//! Pair( Pair( Pair(arg1, arg2), order of Power Derivative), opCode)
typedef map<pair<pair<pair<expr_t, expr_t>, int>, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t;
binary_op_node_map_t binary_op_node_map;
typedef map<pair<pair<pair<expr_t, expr_t>, expr_t>, TrinaryOpcode>, TrinaryOpNode *> trinary_op_node_map_t;
trinary_op_node_map_t trinary_op_node_map;
......@@ -88,7 +89,7 @@ private:
inline expr_t AddPossiblyNegativeConstant(double val);
inline expr_t AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set = 0, const string &arg_exp_info_set_name="");
inline expr_t AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2);
inline expr_t AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2, int powerDerivOrder = 0);
inline expr_t AddTrinaryOp(expr_t arg1, TrinaryOpcode op_code, expr_t arg2, expr_t arg3);
public:
......@@ -137,6 +138,8 @@ public:
expr_t AddDifferent(expr_t iArg1, expr_t iArg2);
//! Adds "arg1^arg2" to model tree
expr_t AddPower(expr_t iArg1, expr_t iArg2);
//! Adds "getPowerDeriv(arg1, arg2, powerDerivOrder)" to model tree
expr_t AddPowerDeriv(expr_t iArg1, expr_t iArg2, int powerDerivOrder);
//! Adds "E(arg1)(arg2)" to model tree
expr_t AddExpectation(int iArg1, expr_t iArg2);
//! Adds "E(arg1)(arg2)" to model tree
......@@ -212,6 +215,10 @@ public:
//! Returns the minimum lag (as a negative number) of the given symbol in the whole data tree (and not only in the equations !!)
/*! Returns 0 if the symbol is not used */
int minLagForSymbol(int symb_id) const;
//! Write the Header for getPowerDeriv when use_dll is used
void writePowerDerivCHeader(ostream &output) const;
//! Write getPowerDeriv
void writePowerDeriv(ostream &output, bool use_dll) const;
//! Thrown when trying to access an unknown variable by deriv_id
class UnknownDerivIDException
{
......@@ -287,9 +294,9 @@ DataTree::AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set, cons
}
inline expr_t
DataTree::AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2)
DataTree::AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2, int powerDerivOrder)
{
binary_op_node_map_t::iterator it = binary_op_node_map.find(make_pair(make_pair(arg1, arg2), op_code));
binary_op_node_map_t::iterator it = binary_op_node_map.find(make_pair(make_pair(make_pair(arg1, arg2), powerDerivOrder), op_code));
if (it != binary_op_node_map.end())
return it->second;
......@@ -298,13 +305,13 @@ DataTree::AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2)
{
double argval1 = arg1->eval(eval_context_t());
double argval2 = arg2->eval(eval_context_t());
double val = BinaryOpNode::eval_opcode(argval1, op_code, argval2);
double val = BinaryOpNode::eval_opcode(argval1, op_code, argval2, powerDerivOrder);
return AddPossiblyNegativeConstant(val);
}
catch (ExprNode::EvalException &e)
{
}
return new BinaryOpNode(*this, arg1, op_code, arg2);
return new BinaryOpNode(*this, arg1, op_code, arg2, powerDerivOrder);
}
inline expr_t
......
preprocessor/DynamicModel.cc
View file @ 35450a29
......@@ -795,6 +795,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
default:
break;
}
writePowerDeriv(output, false);
output.close();
}
}
......@@ -1046,6 +1047,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
fendblock.write(code_file, instruction_number);
FEND_ fend;
fend.write(code_file, instruction_number);
writePowerDeriv(code_file, false);
code_file.close();
}
......@@ -1527,7 +1529,7 @@ DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
mDynamicModelFile << "global oo_;" << endl << endl;
writeDynamicModel(mDynamicModelFile, false);
writePowerDeriv(mDynamicModelFile, false);
mDynamicModelFile.close();
}
......@@ -1556,6 +1558,9 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order)
<< "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
<< "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
// Write function definition if oPowerDeriv is used
writePowerDerivCHeader(mDynamicModelFile);
// Writing the function body
writeDynamicModel(mDynamicModelFile, true);
......@@ -1625,6 +1630,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order)
<< " /* Call the C subroutines. */" << endl
<< " Dynamic(y, x, nb_row_x, params, steady_state, it_, residual, g1, v2, v3);" << endl
<< "}" << endl;
writePowerDeriv(mDynamicModelFile, true);
mDynamicModelFile.close();
}
......@@ -2028,6 +2034,8 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " oo_.endo_simul = y';\n";
mDynamicModelFile << "return;\n";
writePowerDeriv(mDynamicModelFile, false);
mDynamicModelFile.close();
writeModelEquationsOrdered_M(dynamic_basename);
......@@ -3434,6 +3442,8 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
paramsDerivsFile << "end" << endl;
writePowerDeriv(paramsDerivsFile, false);
paramsDerivsFile.close();
}
......
preprocessor/ExprNode.cc
View file @ 35450a29
......@@ -2102,9 +2102,22 @@ BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
ExprNode(datatree_arg),
arg1(arg1_arg),
arg2(arg2_arg),
op_code(op_code_arg)
op_code(op_code_arg),
powerDerivOrder(0)
{
datatree.binary_op_node_map[make_pair(make_pair(arg1, arg2), op_code)] = this;
datatree.binary_op_node_map[make_pair(make_pair(make_pair(arg1, arg2), powerDerivOrder), op_code)] = this;
}
BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
BinaryOpcode op_code_arg, const expr_t arg2_arg, int powerDerivOrder_arg) :
ExprNode(datatree_arg),
arg1(arg1_arg),
arg2(arg2_arg),
op_code(op_code_arg),
powerDerivOrder(powerDerivOrder_arg)
{
assert(powerDerivOrder >= 0);
datatree.binary_op_node_map[make_pair(make_pair(make_pair(arg1, arg2), powerDerivOrder), op_code)] = this;
}
void
......@@ -2162,17 +2175,18 @@ BinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2)
return datatree.Zero;
case oPower:
if (darg2 == datatree.Zero)
{
if (darg1 == datatree.Zero)
return datatree.Zero;
else
if(darg1 == datatree.Zero)
return datatree.Zero;
else
if (dynamic_cast<NumConstNode *>(arg2) != NULL)
{
t11 = datatree.AddMinus(arg2, datatree.One);
t12 = datatree.AddPower(arg1, t11);
t13 = datatree.AddTimes(arg2, t12);
return datatree.AddTimes(darg1, t13);
}
}
else
return datatree.AddTimes(darg1, datatree.AddPowerDeriv(arg1, arg2, powerDerivOrder + 1));
else
{
t11 = datatree.AddLog(arg1);
......@@ -2182,6 +2196,38 @@ BinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2)
t15 = datatree.AddPlus(t12, t14);
return datatree.AddTimes(t15, this);
}
case oPowerDeriv:
if (darg2 == datatree.Zero)
return datatree.AddTimes(darg1, datatree.AddPowerDeriv(arg1, arg2, powerDerivOrder + 1));
else
{
t11 = datatree.AddTimes(darg2, datatree.AddLog(arg1));
t12 = datatree.AddMinus(arg2, datatree.AddPossiblyNegativeConstant(powerDerivOrder));
t13 = datatree.AddTimes(darg1, t12);
t14 = datatree.AddDivide(t13, arg1);
t15 = datatree.AddPlus(t11, t14);
expr_t f = datatree.AddPower(arg1, t12);
expr_t first_part = datatree.AddTimes(f, t15);
for (int i=0; i<powerDerivOrder; i++)
first_part = datatree.AddTimes(first_part, datatree.AddMinus(arg2, datatree.AddPossiblyNegativeConstant(i)));
t13 = datatree.Zero;
for (int i=0; i<powerDerivOrder; i++)
{
t11 = datatree.One;
for (int j=0; j<powerDerivOrder; j++)
if (i != j)
{
t12 = datatree.AddMinus(arg2, datatree.AddPossiblyNegativeConstant(j));
t11 = datatree.AddTimes(t11, t12);
}
t13 = datatree.AddPlus(t13, t11);
}
t13 = datatree.AddTimes(darg2, t13);
t14 = datatree.AddTimes(f, t13);
return datatree.AddPlus(first_part, t14);
}
case oMax:
t11 = datatree.AddGreater(arg1, arg2);
t12 = datatree.AddTimes(t11, darg1);
......@@ -2201,6 +2247,23 @@ BinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2)
exit(EXIT_FAILURE);
}
expr_t
BinaryOpNode::unpackPowerDeriv() const
{
if (op_code != oPowerDeriv)
return const_cast<BinaryOpNode *>(this);
expr_t front = datatree.One;
for (int i=0; i<powerDerivOrder; i++)
front = datatree.AddTimes(front,
datatree.AddMinus(arg2,
datatree.AddPossiblyNegativeConstant(i)));
expr_t tmp = datatree.AddPower(arg1,
datatree.AddMinus(arg2,
datatree.AddPossiblyNegativeConstant(powerDerivOrder)));
return datatree.AddTimes(front, tmp);
}
expr_t
BinaryOpNode::computeDerivative(int deriv_id)
{
......@@ -2236,6 +2299,7 @@ BinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_t
case oDivide:
return 4;
case oPower:
case oPowerDeriv:
if (IS_C(output_type))
// In C, power operator is of the form pow(a, b)
return 100;
......@@ -2281,6 +2345,7 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con
case oDivide:
return cost + 990;
case oPower:
case oPowerDeriv:
return cost + 1160;
case oEqual:
return cost;
......@@ -2306,6 +2371,7 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con
case oDivide:
return cost + 15;
case oPower:
case oPowerDeriv:
return cost + 520;
case oEqual:
return cost;
......@@ -2371,7 +2437,7 @@ BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
}
double
BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException, EvalExternalFunctionException)
BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2, int derivOrder) throw (EvalException, EvalExternalFunctionException)
{
switch (op_code)
{
......@@ -2385,6 +2451,18 @@ BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (Eva
return (v1 / v2);
case oPower:
return (pow(v1, v2));
case oPowerDeriv:
if (fabs(v1) < NEAR_ZERO && v2 > 0 &&
derivOrder >= v2 &&
fabs(v2-nearbyint(v2)) < NEAR_ZERO)
return 0.0;
else
{
double dxp = pow(v1, v2-derivOrder);
for (int i=0; i<derivOrder; i++)
dxp *= v2--;
return dxp;
}
case oMax:
if (v1 < v2)
return v2;
......@@ -2420,7 +2498,7 @@ BinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
double v1 = arg1->eval(eval_context);
double v2 = arg2->eval(eval_context);
return eval_opcode(v1, op_code, v2);
return eval_opcode(v1, op_code, v2, powerDerivOrder);
}
void
......@@ -2482,6 +2560,22 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
return;
}
// Treat derivative of Power
if (op_code == oPowerDeriv)
{
if (IS_LATEX(output_type))
unpackPowerDeriv()->writeOutput(output, output_type, temporary_terms);
else
{
output << "getPowerDeriv(";
arg1->writeOutput(output, output_type, temporary_terms);
output << ",";
arg2->writeOutput(output, output_type, temporary_terms);
output << "," << powerDerivOrder << ")";
}
return;
}
// Treat special case of power operator in C, and case of max and min operators
if ((op_code == oPower && IS_C(output_type)) || op_code == oMax || op_code == oMin)
{
......@@ -2848,6 +2942,8 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr
return (make_pair(1, (expr_t) NULL));
}
break;
case oPowerDeriv:
exit(EXIT_FAILURE);
case oEqual:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
{
......@@ -2968,6 +3064,8 @@ BinaryOpNode::buildSimilarBinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, DataTre
return alt_datatree.AddEqualEqual(alt_arg1, alt_arg2);
case oDifferent:
return alt_datatree.AddDifferent(alt_arg1, alt_arg2);
case oPowerDeriv:
return alt_datatree.AddBinaryOp(alt_arg1, oPowerDeriv, alt_arg2);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
......
preprocessor/ExprNode.hh
View file @ 35450a29
......@@ -107,6 +107,8 @@ enum ExprNodeOutputType
#define MIN_COST_C (40*4)
#define MIN_COST(is_matlab) ((is_matlab) ? MIN_COST_MATLAB : MIN_COST_C)
#define NEAR_ZERO (1e-12)
//! Base class for expression nodes
class ExprNode
{
......@@ -578,9 +580,12 @@ private:
virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
//! Returns the derivative of this node if darg1 and darg2 are the derivatives of the arguments
expr_t composeDerivatives(expr_t darg1, expr_t darg2);
const int powerDerivOrder;
public:
BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
BinaryOpcode op_code_arg, const expr_t arg2_arg);
BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
BinaryOpcode op_code_arg, const expr_t arg2_arg, int powerDerivOrder);
virtual void prepareForDerivation();
virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
......@@ -600,7 +605,7 @@ public:
int equation) const;
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException, EvalExternalFunctionException );
static double eval_opcode(double v1, BinaryOpcode op_code, double v2, int derivOrder) throw (EvalException, EvalExternalFunctionException );
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
......@@ -622,6 +627,11 @@ public:
{
return (op_code);
};
int
get_power_deriv_order() const
{
return powerDerivOrder;
}
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
......@@ -644,6 +654,8 @@ public:
virtual expr_t detrend(int symb_id, expr_t trend) const;
virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
//! Function to write out the oPowerNode in expr_t terms as opposed to writing out the function itself
expr_t unpackPowerDeriv() const;
};
//! Trinary operator node
......
preprocessor/StaticModel.cc
View file @ 35450a29
......@@ -407,6 +407,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
default:
break;
}
writePowerDeriv(output, false);
output.close();
}
}
......@@ -574,6 +575,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
fendblock.write(code_file, instruction_number);
FEND_ fend;
fend.write(code_file, instruction_number);
writePowerDeriv(code_file, false);
code_file.close();
}
......@@ -1242,6 +1244,7 @@ StaticModel::writeStaticMFile(const string &func_name) const
output << " g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl;
output << "end;" << endl; // Close the if nargout >= 3 statement
writePowerDeriv(output, false);
output.close();
}
......@@ -1325,8 +1328,8 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
}
output << " end" << endl
<< "end" << endl;
writePowerDeriv(output, false);
output.close();
}
void
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment