Verified Commit e3a3992c authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Derivation engine w.r.t. parameters generalized to any order

Also, no longer compute two times symmetric elements in derivation w.r.t.
parameters at order 2, for consistency with derivation w.r.t. endogenous.
It is therefore now necessary to duplicate them in the output to keep behavior
unchanged.
parent 6fa115ae
......@@ -5438,6 +5438,22 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
hessian1_output << ";" << endl;
i++;
if (param1 != param2)
{
// Treat symmetric elements
hessian1_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
<< RIGHT_ARRAY_SUBSCRIPT(output_type)
<< "=rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",4"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
i++;
}
}
i = 1;
......@@ -5465,6 +5481,24 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
third_derivs_output << ";" << endl;
i++;
if (param1 != param2)
{
// Treat symmetric elements
third_derivs_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
<< RIGHT_ARRAY_SUBSCRIPT(output_type)
<< "=gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
i++;
}
}
i = 1;
......
......@@ -2033,72 +2033,50 @@ ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, Expr
void
ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
{
if (!(paramsDerivsOrder == 1 || paramsDerivsOrder == 2))
return;
assert(paramsDerivsOrder >= 1);
set<int> deriv_id_set;
addAllParamDerivId(deriv_id_set);
// First-order derivatives w.r.t. params
for (int param : deriv_id_set)
{
for (int eq = 0; eq < (int) equations.size(); eq++)
{
expr_t d1 = equations[eq]->getDerivative(param);
if (d1 == Zero)
expr_t d = equations[eq]->getDerivative(param);
if (d == Zero)
continue;
params_derivatives[{ 0, 1 }][{ eq, param }] = d1;
params_derivatives[{ 0, 1 }][{ eq, param }] = d;
}
if (paramsDerivsOrder == 2)
for (const auto &it : params_derivatives[{ 0, 1 }])
for (int endoOrd = 1; endoOrd < (int) derivatives.size(); endoOrd++)
for (const auto &it : derivatives[endoOrd])
{
int eq, param1;
tie(eq, param1) = vectorToTuple<2>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
expr_t d = it.second->getDerivative(param);
if (d == Zero)
continue;
params_derivatives[{ 0, 2 }][{ eq, param1, param }] = d2;
vector<int> indices{it.first};
indices.push_back(param);
params_derivatives[{ endoOrd, 1 }][indices] = d;
}
for (const auto &it : derivatives[1])
{
int eq, var;
tie(eq, var) = vectorToTuple<2>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
params_derivatives[{ 1, 1 }][{ eq, var, param }] = d2;
}
if (paramsDerivsOrder == 2)
{
for (const auto &it : params_derivatives[{ 1, 1 }])
// Higher-order derivatives w.r.t. parameters
for (int endoOrd = 0; endoOrd < (int) derivatives.size(); endoOrd++)
for (int paramOrd = 2; paramOrd <= paramsDerivsOrder; paramOrd++)
for (const auto &it : params_derivatives[{ endoOrd, paramOrd-1 }])
for (int param : deriv_id_set)
{
int eq, var, param1;
tie(eq, var, param1) = vectorToTuple<3>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
if (it.first.back() > param)
continue;
params_derivatives[{ 1, 2 }][{ eq, var, param1, param }] = d2;
}
for (const auto &it : derivatives[2])
{
int eq, var1, var2;
tie(eq, var1, var2) = vectorToTuple<3>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
expr_t d = it.second->getDerivative(param);
if (d == Zero)
continue;
params_derivatives[{ 2, 1 }][{ eq, var1, var2, param }] = d2;
}
}
vector<int> indices{it.first};
indices.push_back(param);
// At this point, indices of both endogenous and parameters are sorted in non-decreasing order
params_derivatives[{ endoOrd, paramOrd }][indices] = d;
}
}
......
......@@ -105,7 +105,8 @@ protected:
derivation order w.r.t. parameters). For e.g., { 1, 2 } corresponds to the jacobian
differentiated twice w.r.t. to parameters.
In inner maps, the vector of integers consists of: the equation index, then
the derivation IDs of endogenous (in non-decreasing order), then the IDs of parameters */
the derivation IDs of endogenous (in non-decreasing order),
then the IDs of parameters (in non-decreasing order)*/
map<pair<int,int>, map<vector<int>, expr_t>> params_derivatives;
//! Storage for temporary terms in block/bytecode mode
......
......@@ -2710,6 +2710,22 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
hessian1_output << ";" << endl;
i++;
if (param1 != param2)
{
// Treat symmetric elements
hessian1_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
<< RIGHT_ARRAY_SUBSCRIPT(output_type)
<< "=rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",4"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
i++;
}
}
i = 1;
......@@ -2737,6 +2753,24 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
third_derivs_output << ";" << endl;
i++;
if (param1 != param2)
{
// Treat symmetric elements
third_derivs_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
<< RIGHT_ARRAY_SUBSCRIPT(output_type)
<< "=gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
i++;
}
}
i = 1;
......
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