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

Fix derivatives of STEADY_STATE operator w.r.t. parameters (ticket #128)

time-shift
Sébastien Villemot 2011-01-13 18:08:26 +01:00
parent a6a4b3bc28
commit 9c6d65bc0a
8 changed files with 174 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
preprocessor/CodeInterpreter.hh
View File

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

37
preprocessor/DataTree.cc
View File

@ -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<int> &deriv_id_set)
{
}
int
DataTree::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException)
{

23
preprocessor/DataTree.hh
View File

@ -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<pair<int, int>, 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<pair<expr_t, UnaryOpcode>, pair<int, string> >, 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<expr_t, UnaryOpcode>, pair<pair<int, string>, pair<int, int> > >, 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<pair<pair<pair<expr_t, expr_t>, 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<int> &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

12
preprocessor/DynamicModel.cc
View File

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

3
preprocessor/DynamicModel.hh
View File

@ -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<int> &deriv_id_set);
//! Returns true indicating that this is a dynamic model
virtual bool

107
preprocessor/ExprNode.cc
View File

@ -1213,16 +1213,19 @@ VariableNode::removeTrendLeadLag(map<int, expr_t> 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<VariableNode *>(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<VariableNode *>(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<VariableNode *>(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<VariableNode *>(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<pair<int, pair<expr_t, expr_
return (make_pair(1, (expr_t) NULL));
case oSteadyState:
return (make_pair(1, (expr_t) NULL));
case oSteadyStateParamDeriv:
cerr << "UnaryOpNode::normalizeEquation: oSteadyStateParamDeriv not handled" << endl;
exit(EXIT_FAILURE);
case oSteadyStateParam2ndDeriv:
cerr << "UnaryOpNode::normalizeEquation: oSteadyStateParam2ndDeriv not handled" << endl;
exit(EXIT_FAILURE);
case oExpectation:
cerr << "UnaryOpNode::normalizeEquation: oExpectation not handled" << endl;
exit(EXIT_FAILURE);
@ -1832,6 +1913,12 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_
return (make_pair(0, datatree.AddSqrt(New_expr_t)));
case oSteadyState:
return (make_pair(0, datatree.AddSteadyState(New_expr_t)));
case oSteadyStateParamDeriv:
cerr << "UnaryOpNode::normalizeEquation: oSteadyStateParamDeriv not handled" << endl;
exit(EXIT_FAILURE);
case oSteadyStateParam2ndDeriv:
cerr << "UnaryOpNode::normalizeEquation: oSteadyStateParam2ndDeriv not handled" << endl;
exit(EXIT_FAILURE);
case oExpectation:
cerr << "UnaryOpNode::normalizeEquation: oExpectation not handled" << endl;
exit(EXIT_FAILURE);
@ -1846,7 +1933,7 @@ expr_t
UnaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &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:

9
preprocessor/ExprNode.hh
View File

@ -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<expr_t, int> &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;

6
preprocessor/ModelTree.hh
View File

@ -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<bool> BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector<int> &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