diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh index 3f1b55523..76530b0f4 100644 --- a/preprocessor/CodeInterpreter.hh +++ b/preprocessor/CodeInterpreter.hh @@ -1,4 +1,3 @@ - /* * 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, diff --git a/preprocessor/DataTree.cc b/preprocessor/DataTree.cc index a967a3ed6..db60f58be 100644 --- a/preprocessor/DataTree.cc +++ b/preprocessor/DataTree.cc @@ -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 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; +} diff --git a/preprocessor/DataTree.hh b/preprocessor/DataTree.hh index f49460ae8..7e5076499 100644 --- a/preprocessor/DataTree.hh +++ b/preprocessor/DataTree.hh @@ -62,7 +62,8 @@ protected: //! Pair( Pair (arg1, UnaryOpCode), Pair(Expectation Info Set, Expectation Info Set Name) ) typedef map, pair >, UnaryOpNode *> unary_op_node_map_t; unary_op_node_map_t unary_op_node_map; - typedef map, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t; + //! Pair( Pair( Pair(arg1, arg2), order of Power Derivative), opCode) + typedef map, int>, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t; binary_op_node_map_t binary_op_node_map; typedef map, 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 diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 5013eee79..07ff6e810 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -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(); } diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index 676f0fac3..32b922687 100644 --- a/preprocessor/ExprNode.cc +++ b/preprocessor/ExprNode.cc @@ -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(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(this); + + expr_t front = datatree.One; + for (int i=0; i &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; ieval(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 &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 > &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 normalizeEquation(int symb_id_endo, vector > > &List_of_Op_RHS) const; virtual expr_t getChainRuleDerivative(int deriv_id, const map &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 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 diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 9488efda7..73e77048f 100644 --- a/preprocessor/StaticModel.cc +++ b/preprocessor/StaticModel.cc @@ -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