From 9c6d65bc0a4742f0cc30eaa129567d11a5fc12a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 13 Jan 2011 18:08:26 +0100 Subject: [PATCH] Fix derivatives of STEADY_STATE operator w.r.t. parameters (ticket #128) --- preprocessor/CodeInterpreter.hh | 4 +- preprocessor/DataTree.cc | 37 ++++++++++- preprocessor/DataTree.hh | 23 ++++--- preprocessor/DynamicModel.cc | 12 +++- preprocessor/DynamicModel.hh | 3 +- preprocessor/ExprNode.cc | 107 +++++++++++++++++++++++++++++--- preprocessor/ExprNode.hh | 9 ++- preprocessor/ModelTree.hh | 6 +- 8 files changed, 174 insertions(+), 27 deletions(-) diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh index ef9b46d5e..76f4ae105 100644 --- a/preprocessor/CodeInterpreter.hh +++ b/preprocessor/CodeInterpreter.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2007-2010 Dynare Team + * Copyright (C) 2007-2011 Dynare Team * * This file is part of Dynare. * @@ -190,6 +190,8 @@ enum UnaryOpcode oAtanh, oSqrt, oSteadyState, + oSteadyStateParamDeriv, // for the derivative of the STEADY_STATE operator w.r.t. to a parameter + oSteadyStateParam2ndDeriv, // for the 2nd derivative of the STEADY_STATE operator w.r.t. to a parameter oExpectation, oErf }; diff --git a/preprocessor/DataTree.cc b/preprocessor/DataTree.cc index 6b30f1eb6..21f322b98 100644 --- a/preprocessor/DataTree.cc +++ b/preprocessor/DataTree.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2010 Dynare Team + * Copyright (C) 2003-2011 Dynare Team * * This file is part of Dynare. * @@ -440,6 +440,18 @@ DataTree::AddSteadyState(expr_t iArg1) return AddUnaryOp(oSteadyState, iArg1); } +expr_t +DataTree::AddSteadyStateParamDeriv(expr_t iArg1, int param_symb_id) +{ + return AddUnaryOp(oSteadyStateParamDeriv, iArg1, 0, "", param_symb_id); +} + +expr_t +DataTree::AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id) +{ + return AddUnaryOp(oSteadyStateParam2ndDeriv, iArg1, 0, "", param1_symb_id, param2_symb_id); +} + expr_t DataTree::AddExpectation(int iArg1, expr_t iArg2) { @@ -532,6 +544,29 @@ DataTree::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException) throw UnknownDerivIDException(); } +SymbolType +DataTree::getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException) +{ + throw UnknownDerivIDException(); +} + +int +DataTree::getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException) +{ + throw UnknownDerivIDException(); +} + +int +DataTree::getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException) +{ + throw UnknownDerivIDException(); +} + +void +DataTree::addAllParamDerivId(set &deriv_id_set) +{ +} + int DataTree::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException) { diff --git a/preprocessor/DataTree.hh b/preprocessor/DataTree.hh index 7e5076499..596ce0671 100644 --- a/preprocessor/DataTree.hh +++ b/preprocessor/DataTree.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2010 Dynare Team + * Copyright (C) 2003-2011 Dynare Team * * This file is part of Dynare. * @@ -59,8 +59,8 @@ protected: //! Pair (symbol_id, lag) used as key typedef map, VariableNode *> variable_node_map_t; variable_node_map_t variable_node_map; - //! Pair( Pair (arg1, UnaryOpCode), Pair(Expectation Info Set, Expectation Info Set Name) ) - typedef map, pair >, UnaryOpNode *> unary_op_node_map_t; + //! Pair( Pair (arg1, UnaryOpCode), Pair(Pair(Expectation Info Set, Expectation Info Set Name), (param1_symb_id, param2_symb_id)) ) + typedef map, pair, pair > >, UnaryOpNode *> unary_op_node_map_t; unary_op_node_map_t unary_op_node_map; //! Pair( Pair( Pair(arg1, arg2), order of Power Derivative), opCode) typedef map, int>, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t; @@ -88,7 +88,7 @@ private: int node_counter; 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 AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set = 0, const string &arg_exp_info_set_name="", int param1_symb_id = 0, int param2_symb_id = 0); 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); @@ -188,6 +188,10 @@ public: expr_t AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3); //! Adds "steadyState(arg)" to model tree expr_t AddSteadyState(expr_t iArg1); + //! Add derivative of steady state w.r.t. parameter to model tree + expr_t AddSteadyStateParamDeriv(expr_t iArg1, int param_symb_id); + //! Add 2nd derivative of steady state w.r.t. parameter to model tree + expr_t AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id); //! Adds "arg1=arg2" to model tree expr_t AddEqual(expr_t iArg1, expr_t iArg2); //! Adds a model local variable with its value @@ -236,8 +240,13 @@ public: //! Returns the derivation ID, or throws an exception if the derivation ID does not exist virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); + virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException); + virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException); + virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException); //! Returns the column of the dynamic Jacobian associated to a derivation ID virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException); + //! Adds to the set all the deriv IDs corresponding to parameters + virtual void addAllParamDerivId(set &deriv_id_set); //! Returns bool indicating whether DataTree represents a Dynamic Model (returns true in DynamicModel.hh) virtual bool @@ -268,10 +277,10 @@ DataTree::AddPossiblyNegativeConstant(double v) } inline expr_t -DataTree::AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set, const string &arg_exp_info_set_name) +DataTree::AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set, const string &arg_exp_info_set_name, int param1_symb_id, int param2_symb_id) { // If the node already exists in tree, share it - unary_op_node_map_t::iterator it = unary_op_node_map.find(make_pair(make_pair(arg, op_code), make_pair(arg_exp_info_set, arg_exp_info_set_name))); + unary_op_node_map_t::iterator it = unary_op_node_map.find(make_pair(make_pair(arg, op_code), make_pair(make_pair(arg_exp_info_set, arg_exp_info_set_name), make_pair(param1_symb_id, param2_symb_id)))); if (it != unary_op_node_map.end()) return it->second; @@ -290,7 +299,7 @@ DataTree::AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set, cons { } } - return new UnaryOpNode(*this, op_code, arg, arg_exp_info_set, arg_exp_info_set_name); + return new UnaryOpNode(*this, op_code, arg, arg_exp_info_set, arg_exp_info_set_name, param1_symb_id, param2_symb_id); } inline expr_t diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index ec81d2293..9ca30b89d 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2010 Dynare Team + * Copyright (C) 2003-2011 Dynare Team * * This file is part of Dynare. * @@ -3110,6 +3110,14 @@ DynamicModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDExcept return it->second; } +void +DynamicModel::addAllParamDerivId(set &deriv_id_set) +{ + for (size_t i = 0; i < inv_deriv_id_table.size(); i++) + if (symbol_table.getType(inv_deriv_id_table[i].first) == eParameter) + deriv_id_set.insert(i); +} + void DynamicModel::computeDynJacobianCols(bool jacobianExo) { @@ -3310,7 +3318,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const cerr << "ERROR: Can't open file " << filename << " for writing" << endl; exit(EXIT_FAILURE); } - paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_params_derivs(y, x, params, it_)" << endl + paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_params_derivs(y, x, params, it_, ss_param_deriv, ss_param_2nd_deriv)" << endl << "%" << endl << "% Warning : this file is generated automatically by Dynare" << endl << "% from model file (.mod)" << endl << endl; diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh index 4d98a97d4..82be0a345 100644 --- a/preprocessor/DynamicModel.hh +++ b/preprocessor/DynamicModel.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2010 Dynare Team + * Copyright (C) 2003-2011 Dynare Team * * This file is part of Dynare. * @@ -293,6 +293,7 @@ public: virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException); + virtual void addAllParamDerivId(set &deriv_id_set); //! Returns true indicating that this is a dynamic model virtual bool diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index 13c139b6d..ec907ab99 100644 --- a/preprocessor/ExprNode.cc +++ b/preprocessor/ExprNode.cc @@ -1213,16 +1213,19 @@ VariableNode::removeTrendLeadLag(map trend_symbols_map) const } } -UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg) : +UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, const string &expectation_information_set_name_arg, int param1_symb_id_arg, int param2_symb_id_arg) : ExprNode(datatree_arg), arg(arg_arg), expectation_information_set(expectation_information_set_arg), expectation_information_set_name(expectation_information_set_name_arg), + param1_symb_id(param1_symb_id_arg), + param2_symb_id(param2_symb_id_arg), op_code(op_code_arg) { // Add myself to the unary op map datatree.unary_op_node_map[make_pair(make_pair(arg, op_code), - make_pair(expectation_information_set, expectation_information_set_name))] = this; + make_pair(make_pair(expectation_information_set, expectation_information_set_name), + make_pair(param1_symb_id, param2_symb_id)))] = this; } void @@ -1235,12 +1238,15 @@ UnaryOpNode::prepareForDerivation() arg->prepareForDerivation(); - // Non-null derivatives are those of the argument + // Non-null derivatives are those of the argument (except for STEADY_STATE) non_null_derivatives = arg->non_null_derivatives; + if (op_code == oSteadyState || op_code == oSteadyStateParamDeriv + || op_code == oSteadyStateParam2ndDeriv) + datatree.addAllParamDerivId(non_null_derivatives); } expr_t -UnaryOpNode::composeDerivatives(expr_t darg) +UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id) { expr_t t11, t12, t13; @@ -1304,9 +1310,45 @@ UnaryOpNode::composeDerivatives(expr_t darg) return datatree.AddDivide(darg, t11); case oSteadyState: if (datatree.isDynamic()) - return datatree.Zero; + { + if (datatree.getTypeByDerivID(deriv_id) == eParameter) + { + VariableNode *varg = dynamic_cast(arg); + if (varg == NULL) + { + cerr << "UnaryOpNode::writeOutput: STEADY_STATE() should only be used on standalone variables (like STEADY_STATE(y)) to be derivable w.r.t. parameters" << endl; + exit(EXIT_FAILURE); + } + if (datatree.symbol_table.getType(varg->symb_id) == eEndogenous) + return datatree.AddSteadyStateParamDeriv(arg, datatree.getSymbIDByDerivID(deriv_id)); + else + return datatree.Zero; + } + else + return datatree.Zero; + } else return darg; + case oSteadyStateParamDeriv: + assert(datatree.isDynamic()); + if (datatree.getTypeByDerivID(deriv_id) == eParameter) + { + VariableNode *varg = dynamic_cast(arg); + assert(varg != NULL); + assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous); + return datatree.AddSteadyStateParam2ndDeriv(arg, param1_symb_id, datatree.getSymbIDByDerivID(deriv_id)); + } + else + return datatree.Zero; + case oSteadyStateParam2ndDeriv: + assert(datatree.isDynamic()); + if (datatree.getTypeByDerivID(deriv_id) == eParameter) + { + cerr << "3rd derivative of STEADY_STATE node w.r.t. three parameters not implemented" << endl; + exit(EXIT_FAILURE); + } + else + return datatree.Zero; case oExpectation: cerr << "UnaryOpNode::composeDerivatives: not implemented on oExpectation" << endl; exit(EXIT_FAILURE); @@ -1330,7 +1372,7 @@ expr_t UnaryOpNode::computeDerivative(int deriv_id) { expr_t darg = arg->getDerivative(deriv_id); - return composeDerivatives(darg); + return composeDerivatives(darg, deriv_id); } int @@ -1381,6 +1423,8 @@ UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) cons case oSqrt: return cost + 570; case oSteadyState: + case oSteadyStateParamDeriv: + case oSteadyStateParam2ndDeriv: case oExpectation: return cost; } @@ -1419,6 +1463,8 @@ UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) cons case oSqrt: return cost + 90; case oSteadyState: + case oSteadyStateParamDeriv: + case oSteadyStateParam2ndDeriv: case oExpectation: return cost; } @@ -1582,6 +1628,33 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, } arg->writeOutput(output, new_output_type, temporary_terms); return; + case oSteadyStateParamDeriv: + { + VariableNode *varg = dynamic_cast(arg); + assert(varg != NULL); + assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous); + assert(datatree.symbol_table.getType(param1_symb_id) == eParameter); + int tsid_endo = datatree.symbol_table.getTypeSpecificID(varg->symb_id); + int tsid_param = datatree.symbol_table.getTypeSpecificID(param1_symb_id); + assert(IS_MATLAB(output_type)); + output << "ss_param_deriv(" << tsid_endo+1 << "," << tsid_param+1 << ")"; + } + return; + case oSteadyStateParam2ndDeriv: + { + VariableNode *varg = dynamic_cast(arg); + assert(varg != NULL); + assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous); + assert(datatree.symbol_table.getType(param1_symb_id) == eParameter); + assert(datatree.symbol_table.getType(param2_symb_id) == eParameter); + int tsid_endo = datatree.symbol_table.getTypeSpecificID(varg->symb_id); + int tsid_param1 = datatree.symbol_table.getTypeSpecificID(param1_symb_id); + int tsid_param2 = datatree.symbol_table.getTypeSpecificID(param2_symb_id); + assert(IS_MATLAB(output_type)); + output << "ss_param_2nd_deriv(" << tsid_endo+1 << "," << tsid_param1+1 + << "," << tsid_param2+1 << ")"; + } + return; case oExpectation: cerr << "UnaryOpNode::writeOutput: not implemented on oExpectation" << endl; exit(EXIT_FAILURE); @@ -1675,6 +1748,8 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) throw (EvalException, Ev return (sqrt(v)); case oSteadyState: return (v); + case oSteadyStateParamDeriv: + case oSteadyStateParam2ndDeriv: case oExpectation: throw EvalException(); case oErf: @@ -1785,6 +1860,12 @@ UnaryOpNode::normalizeEquation(int var_endo, vector &recursive_variables) { expr_t darg = arg->getChainRuleDerivative(deriv_id, recursive_variables); - return composeDerivatives(darg); + return composeDerivatives(darg, deriv_id); } expr_t @@ -1890,6 +1977,12 @@ UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) con return alt_datatree.AddSqrt(alt_arg); case oSteadyState: return alt_datatree.AddSteadyState(alt_arg); + case oSteadyStateParamDeriv: + cerr << "UnaryOpNode::buildSimilarUnaryOpNode: oSteadyStateParamDeriv can't be translated" << endl; + exit(EXIT_FAILURE); + case oSteadyStateParam2ndDeriv: + cerr << "UnaryOpNode::buildSimilarUnaryOpNode: oSteadyStateParam2ndDeriv can't be translated" << endl; + exit(EXIT_FAILURE); case oExpectation: return alt_datatree.AddExpectation(expectation_information_set, alt_arg); case oErf: diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh index ea25fb1ff..15bab4f0a 100644 --- a/preprocessor/ExprNode.hh +++ b/preprocessor/ExprNode.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2007-2010 Dynare Team + * Copyright (C) 2007-2011 Dynare Team * * This file is part of Dynare. * @@ -442,6 +442,7 @@ public: //! Symbol or variable node class VariableNode : public ExprNode { + friend class UnaryOpNode; private: //! Id from the symbol table const int symb_id; @@ -505,13 +506,15 @@ private: const int expectation_information_set; //! Stores the information set name. Only used for expectation operator const string expectation_information_set_name; + //! Only used for oSteadyStateParamDeriv and oSteadyStateParam2ndDeriv + const int param1_symb_id, param2_symb_id; const UnaryOpcode op_code; virtual expr_t computeDerivative(int deriv_id); virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const; //! Returns the derivative of this node if darg is the derivative of the argument - expr_t composeDerivatives(expr_t darg); + expr_t composeDerivatives(expr_t darg, int deriv_id); public: - UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg); + UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, const string &expectation_information_set_name_arg, int param1_symb_id_arg, int param2_symb_id_arg); virtual void prepareForDerivation(); virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const; diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh index 11ac9bf7c..e05da9332 100644 --- a/preprocessor/ModelTree.hh +++ b/preprocessor/ModelTree.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2010 Dynare Team + * Copyright (C) 2003-2011 Dynare Team * * This file is part of Dynare. * @@ -183,10 +183,6 @@ protected: //! Determine for each block if it is linear or not vector BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector &variable_reordered) const; - virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0; - virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0; - virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0; - //! Determine the simulation type of each block virtual BlockSimulationType getBlockSimulationType(int block_number) const = 0; //! Return the number of blocks