Commit 571b5d08 authored by Sébastien Villemot's avatar Sébastien Villemot

Computation of temporary terms generalized to any derivation order

parent 67ac4bf8
Pipeline #400 passed with stage
in 1 minute and 29 seconds
......@@ -257,46 +257,6 @@ enum class PriorDistributions
weibull = 8
};
enum class NodeTreeReference
{
residuals,
firstDeriv,
secondDeriv,
thirdDeriv,
residualsParamsDeriv,
jacobianParamsDeriv,
residualsParamsSecondDeriv,
jacobianParamsSecondDeriv,
hessianParamsDeriv
};
/*! Lists elements of the NodeTreeReference enum that come “before” the argument.
Used in AbstractExternalFunctionNode::computeTemporaryTerms */
inline auto
nodeTreeReferencesBefore(NodeTreeReference tr)
{
vector<NodeTreeReference> v;
// Should be same order as the one appearing in ModelTree::computeTemporaryTerms()
for (auto tr2 : { NodeTreeReference::residuals, NodeTreeReference::firstDeriv, NodeTreeReference::secondDeriv, NodeTreeReference::thirdDeriv })
if (tr == tr2)
return v;
else
v.push_back(tr2);
v.clear();
// Should be same order as the one appearing in ModelTree::computeParamsDerivativesTemporaryTerms()
for (auto tr2 : { NodeTreeReference::residualsParamsDeriv, NodeTreeReference::jacobianParamsDeriv, NodeTreeReference::residualsParamsSecondDeriv,
NodeTreeReference::jacobianParamsSecondDeriv, NodeTreeReference::hessianParamsDeriv})
if (tr == tr2)
return v;
else
v.push_back(tr2);
cerr << "nodeTreeReferencesBelow: impossible case" << endl;
exit(EXIT_FAILURE);
}
struct Block_contain_type
{
int Equation, Variable, Own_Derivative;
......
......@@ -5349,8 +5349,8 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
deriv_node_temp_terms_t tef_terms;
writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
for (auto it : { make_pair(0,1), make_pair(1,1), make_pair(0,2), make_pair(1,2), make_pair(2,1) })
writeTemporaryTerms(params_derivs_temporary_terms.find(it)->second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
for (const auto &it : params_derivs_temporary_terms)
writeTemporaryTerms(it.second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
for (const auto & residuals_params_derivative : params_derivatives.find({ 0, 1 })->second)
{
......@@ -6553,8 +6553,8 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
temporary_terms_t temp_term_union;
string concat = "all";
for (auto it : { make_pair(0,1), make_pair(1,1), make_pair(0,2), make_pair(1,2), make_pair(2,1) })
writeJsonTemporaryTerms(params_derivs_temporary_terms.find(it)->second, temp_term_union, model_output, tef_terms, concat);
for (const auto &it : params_derivs_temporary_terms)
writeJsonTemporaryTerms(it.second, temp_term_union, model_output, tef_terms, concat);
jacobian_output << "\"deriv_wrt_params\": {"
<< " \"neqs\": " << equations.size()
......
......@@ -86,7 +86,7 @@ ExprNode::cost(const temporary_terms_t &temp_terms_map, bool is_matlab) const
}
int
ExprNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
ExprNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const
{
// For a terminal node, the cost is null
return 0;
......@@ -146,9 +146,10 @@ ExprNode::collectExogenous(set<pair<int, int>> &result) const
}
void
ExprNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
ExprNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
// Nothing to do for a terminal node
}
......@@ -2169,7 +2170,7 @@ UnaryOpNode::computeDerivative(int deriv_id)
}
int
UnaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
UnaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const
{
// For a temporary term, the cost is null
for (const auto & it : temp_terms_map)
......@@ -2295,17 +2296,18 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
}
void
UnaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
UnaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
expr_t this2 = const_cast<UnaryOpNode *>(this);
auto it = reference_count.find(this2);
if (it == reference_count.end())
{
reference_count[this2] = { 1, tr };
arg->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
reference_count[this2] = { 1, derivOrder };
arg->computeTemporaryTerms(derivOrder, temp_terms_map, reference_count, is_matlab);
}
else
{
......@@ -4057,7 +4059,7 @@ BinaryOpNode::precedenceJson(const temporary_terms_t &temporary_terms) const
}
int
BinaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
BinaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const
{
// For a temporary term, the cost is null
for (const auto & it : temp_terms_map)
......@@ -4142,9 +4144,10 @@ BinaryOpNode::cost(int cost, bool is_matlab) const
}
void
BinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
BinaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
expr_t this2 = const_cast<BinaryOpNode *>(this);
auto it = reference_count.find(this2);
......@@ -4152,9 +4155,9 @@ BinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &r
{
// If this node has never been encountered, set its ref count to one,
// and travel through its children
reference_count[this2] = { 1, tr };
arg1->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
arg2->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
reference_count[this2] = { 1, derivOrder };
arg1->computeTemporaryTerms(derivOrder, temp_terms_map, reference_count, is_matlab);
arg2->computeTemporaryTerms(derivOrder, temp_terms_map, reference_count, is_matlab);
}
else
{
......@@ -5964,7 +5967,7 @@ TrinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_
}
int
TrinaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
TrinaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const
{
// For a temporary term, the cost is null
for (const auto & it : temp_terms_map)
......@@ -6016,9 +6019,10 @@ TrinaryOpNode::cost(int cost, bool is_matlab) const
}
void
TrinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
TrinaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
expr_t this2 = const_cast<TrinaryOpNode *>(this);
auto it = reference_count.find(this2);
......@@ -6026,10 +6030,10 @@ TrinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &
{
// If this node has never been encountered, set its ref count to one,
// and travel through its children
reference_count[this2] = { 1, tr };
arg1->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
arg2->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
arg3->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
reference_count[this2] = { 1, derivOrder };
arg1->computeTemporaryTerms(derivOrder, temp_terms_map, reference_count, is_matlab);
arg2->computeTemporaryTerms(derivOrder, temp_terms_map, reference_count, is_matlab);
arg3->computeTemporaryTerms(derivOrder, temp_terms_map, reference_count, is_matlab);
}
else
{
......@@ -7118,9 +7122,10 @@ AbstractExternalFunctionNode::getIndxInTefTerms(int the_symb_id, const deriv_nod
}
void
AbstractExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
AbstractExternalFunctionNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
/* All external function nodes are declared as temporary terms.
......@@ -7133,18 +7138,17 @@ AbstractExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTr
corresponding to the same external function call is present in that
previous level. */
for (auto tr2 : nodeTreeReferencesBefore(tr))
for (auto &tt : temp_terms_map)
{
auto it = find_if(temp_terms_map[tr2].cbegin(), temp_terms_map[tr2].cend(),
sameTefTermPredicate());
if (it != temp_terms_map[tr2].cend())
auto it = find_if(tt.second.cbegin(), tt.second.cend(), sameTefTermPredicate());
if (it != tt.second.cend())
{
temp_terms_map[tr2].insert(const_cast<AbstractExternalFunctionNode *>(this));
tt.second.insert(const_cast<AbstractExternalFunctionNode *>(this));
return;
}
}
temp_terms_map[tr].insert(const_cast<AbstractExternalFunctionNode *>(this));
temp_terms_map[derivOrder].insert(const_cast<AbstractExternalFunctionNode *>(this));
}
bool
......@@ -8460,9 +8464,10 @@ VarExpectationNode::VarExpectationNode(DataTree &datatree_arg,
}
void
VarExpectationNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
VarExpectationNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
cerr << "VarExpectationNode::computeTemporaryTerms not implemented." << endl;
exit(EXIT_FAILURE);
......@@ -8917,11 +8922,12 @@ PacExpectationNode::PacExpectationNode(DataTree &datatree_arg,
}
void
PacExpectationNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
bool is_matlab, NodeTreeReference tr) const
PacExpectationNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const
{
temp_terms_map[tr].insert(const_cast<PacExpectationNode *>(this));
temp_terms_map[derivOrder].insert(const_cast<PacExpectationNode *>(this));
}
void
......
This diff is collapsed.
......@@ -1315,17 +1315,12 @@ ModelTree::computeDerivatives(int order, const set<int> &vars)
void
ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
{
map<expr_t, pair<int, NodeTreeReference>> reference_count;
temporary_terms.clear();
temporary_terms_mlv.clear();
temporary_terms_derivatives.clear();
temporary_terms_derivatives.resize(4);
/* Collect all model local variables appearing in equations (and only those,
because printing unused model local variables can lead to a crash,
see Dynare/dynare#101).
Then store them in a dedicated structure (temporary_terms_mlv), that will
be treated as the rest of temporary terms. */
temporary_terms_mlv.clear();
set<int> used_local_vars;
for (auto & equation : equations)
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
......@@ -1335,59 +1330,44 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
temporary_terms_mlv[v] = local_variables_table.find(used_local_var)->second;
}
map<NodeTreeReference, temporary_terms_t> temp_terms_map;
temp_terms_map[NodeTreeReference::residuals] = temporary_terms_derivatives[0];
temp_terms_map[NodeTreeReference::firstDeriv] = temporary_terms_derivatives[1];
temp_terms_map[NodeTreeReference::secondDeriv] = temporary_terms_derivatives[2];
temp_terms_map[NodeTreeReference::thirdDeriv] = temporary_terms_derivatives[3];
// Compute the temporary terms in equations and derivatives
map<pair<int, int>, temporary_terms_t> temp_terms_map;
if (!no_tmp_terms)
{
map<expr_t, pair<int, pair<int, int>>> reference_count;
for (auto & equation : equations)
equation->computeTemporaryTerms(reference_count,
equation->computeTemporaryTerms({ 0, 0 },
temp_terms_map,
is_matlab, NodeTreeReference::residuals);
for (auto & first_derivative : derivatives[1])
first_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
is_matlab, NodeTreeReference::firstDeriv);
for (auto & second_derivative : derivatives[2])
second_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
is_matlab, NodeTreeReference::secondDeriv);
for (auto & third_derivative : derivatives[3])
third_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
is_matlab, NodeTreeReference::thirdDeriv);
reference_count,
is_matlab);
for (int order = 1; order < (int) derivatives.size(); order++)
for (const auto &it : derivatives[order])
it.second->computeTemporaryTerms({ 0, order },
temp_terms_map,
reference_count,
is_matlab);
}
for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
it != temp_terms_map.end(); it++)
temporary_terms.insert(it->second.begin(), it->second.end());
// Fill the (now obsolete) temporary_terms structure
temporary_terms.clear();
for (const auto &it : temp_terms_map)
temporary_terms.insert(it.second.begin(), it.second.end());
temporary_terms_derivatives[0] = temp_terms_map[NodeTreeReference::residuals];
temporary_terms_derivatives[1] = temp_terms_map[NodeTreeReference::firstDeriv];
temporary_terms_derivatives[2] = temp_terms_map[NodeTreeReference::secondDeriv];
temporary_terms_derivatives[3] = temp_terms_map[NodeTreeReference::thirdDeriv];
// Fill the new structure
temporary_terms_derivatives.clear();
temporary_terms_derivatives.resize(derivatives.size());
for (int order = 0; order < (int) derivatives.size(); order++)
temporary_terms_derivatives[order] = move(temp_terms_map[{ 0, order }]);
// Compute indices in MATLAB/Julia vector
int idx = 0;
for (auto &it : temporary_terms_mlv)
temporary_terms_idxs[it.first] = idx++;
for (auto it : temporary_terms_derivatives[0])
temporary_terms_idxs[it] = idx++;
for (auto it : temporary_terms_derivatives[1])
temporary_terms_idxs[it] = idx++;
for (auto it : temporary_terms_derivatives[2])
temporary_terms_idxs[it] = idx++;
for (auto it : temporary_terms_derivatives[3])
temporary_terms_idxs[it] = idx++;
for (int order = 0; order < (int) derivatives.size(); order++)
for (const auto &it : temporary_terms_derivatives[order])
temporary_terms_idxs[it] = idx++;
}
void
......@@ -2078,66 +2058,24 @@ ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
void
ModelTree::computeParamsDerivativesTemporaryTerms()
{
map<expr_t, pair<int, NodeTreeReference >> reference_count;
map<NodeTreeReference, temporary_terms_t> temp_terms_map;
temp_terms_map[NodeTreeReference::residualsParamsDeriv] = params_derivs_temporary_terms[{ 0, 1 }];
temp_terms_map[NodeTreeReference::jacobianParamsDeriv] = params_derivs_temporary_terms[{ 1, 1 }];
temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv] = params_derivs_temporary_terms[{ 0, 2 }];
temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv] = params_derivs_temporary_terms[{ 1, 2 }];
temp_terms_map[NodeTreeReference::hessianParamsDeriv] = params_derivs_temporary_terms[{ 2, 1}];
map<expr_t, pair<int, pair<int, int>>> reference_count;
/* The temp terms should be constructed in the same order as the for loops in
{Static,Dynamic}Model::write{Json,}ParamsDerivativesFile() */
for (const auto &residuals_params_derivative : params_derivatives[{ 0, 1 }])
residuals_params_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::residualsParamsDeriv);
for (const auto &jacobian_params_derivative : params_derivatives[{ 1, 1 }])
jacobian_params_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::jacobianParamsDeriv);
for (const auto &it : params_derivatives[{ 0, 2 }])
it.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::residualsParamsSecondDeriv);
for (const auto &it : params_derivatives[{ 1, 2 }])
it.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::jacobianParamsSecondDeriv);
for (const auto &it : params_derivatives[{ 2, 1 }])
it.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::hessianParamsDeriv);
params_derivs_temporary_terms[{ 0, 1 }] = temp_terms_map[NodeTreeReference::residualsParamsDeriv];
params_derivs_temporary_terms[{ 1, 1 }] = temp_terms_map[NodeTreeReference::jacobianParamsDeriv];
params_derivs_temporary_terms[{ 0, 2 }] = temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv];
params_derivs_temporary_terms[{ 1, 2 }] = temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv];
params_derivs_temporary_terms[{ 2, 1 }] = temp_terms_map[NodeTreeReference::hessianParamsDeriv];
params_derivs_temporary_terms.clear();
for (const auto &it : params_derivatives)
for (const auto &it2 : it.second)
it2.second->computeTemporaryTerms(it.first,
params_derivs_temporary_terms,
reference_count,
true);
int idx = 0;
for (auto &it : temporary_terms_mlv)
params_derivs_temporary_terms_idxs[it.first] = idx++;
for (auto tt : params_derivs_temporary_terms[{ 0, 1 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
for (auto tt : params_derivs_temporary_terms[{ 1, 1 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
for (auto tt : params_derivs_temporary_terms[{ 0, 2 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
for (auto tt : params_derivs_temporary_terms[{ 1, 2 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
for (auto tt : params_derivs_temporary_terms[{ 2, 1 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
for (const auto &it : params_derivs_temporary_terms)
for (const auto &tt : it.second)
params_derivs_temporary_terms_idxs[tt] = idx++;
}
bool
......
......@@ -120,12 +120,14 @@ protected:
/*! Index 0 is temp. terms of residuals, index 1 for first derivatives, ... */
vector<temporary_terms_t> temporary_terms_derivatives;
//! Stores, for each temporary term, its index in the MATLAB/Julia vector
temporary_terms_idxs_t temporary_terms_idxs;
//! Temporary terms for parameter derivatives, under a disaggregated form
/*! The pair of integers is to be interpreted as in param_derivatives */
map<pair<int,int>, temporary_terms_t> params_derivs_temporary_terms;
//! Stores, for each temporary term in param. derivs, its index in the MATLAB/Julia vector
temporary_terms_idxs_t params_derivs_temporary_terms_idxs;
//! Trend variables and their growth factors
......
......@@ -2615,8 +2615,8 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
deriv_node_temp_terms_t tef_terms;
writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
for (auto it : { make_pair(0,1), make_pair(1,1), make_pair(0,2), make_pair(1,2), make_pair(2,1) })
writeTemporaryTerms(params_derivs_temporary_terms.find(it)->second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
for (const auto &it : params_derivs_temporary_terms)
writeTemporaryTerms(it.second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
for (const auto & residuals_params_derivative : params_derivatives.find({ 0, 1 })->second)
{
......@@ -3082,8 +3082,8 @@ StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
temporary_terms_t temp_term_union;
string concat = "all";
for (auto it : { make_pair(0,1), make_pair(1,1), make_pair(0,2), make_pair(1,2), make_pair(2,1) })
writeJsonTemporaryTerms(params_derivs_temporary_terms.find(it)->second, temp_term_union, model_output, tef_terms, concat);
for (const auto &it : params_derivs_temporary_terms)
writeJsonTemporaryTerms(it.second, temp_term_union, model_output, tef_terms, concat);
jacobian_output << "\"deriv_wrt_params\": {"
<< " \"neqs\": " << equations.size()
......
Markdown is supported
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