From e89dd9cb8bebdc194a0045a4b9ca3926db86d74d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 15 Nov 2018 16:39:53 +0100
Subject: [PATCH] Make storage of derivatives ready for k-order
This is the first step towards computing k-order derivatives.
Ref Dynare/dynare#217, #10
---
src/DynamicModel.cc | 287 +++++++++++++++++++++-----------------------
src/ModelTree.cc | 268 ++++++++++++++++++-----------------------
src/ModelTree.hh | 114 +++++++-----------
src/StaticModel.cc | 272 ++++++++++++++++++++---------------------
4 files changed, 433 insertions(+), 508 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 912f261e..2a1094dd 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -245,8 +245,8 @@ DynamicModel::operator=(const DynamicModel &m)
void
DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
{
- auto it = first_derivatives.find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
- if (it != first_derivatives.end())
+ auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
+ if (it != derivatives[1].end())
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
else
{
@@ -274,7 +274,6 @@ DynamicModel::computeTemporaryTermsOrdered()
map<expr_t, pair<int, int>> first_occurence;
map<expr_t, int> reference_count;
BinaryOpNode *eq_node;
- first_derivatives_t::const_iterator it;
first_chain_rule_derivatives_t::const_iterator it_chr;
ostringstream tmp_s;
v_temporary_terms.clear();
@@ -1016,10 +1015,10 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
map<pair< int, pair<int, int>>, expr_t> first_derivatives_reordered_endo;
map<pair< pair<int, SymbolType>, pair<int, int>>, expr_t> first_derivatives_reordered_exo;
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
- unsigned int eq = first_derivative.first.first;
+ int deriv_id = first_derivative.first[1];
+ unsigned int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
@@ -1103,24 +1102,24 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
fjmp_if_eval.write(code_file, instruction_number);
int prev_instruction_number = instruction_number;
- vector<vector<pair<pair<int, int>, int >>> derivatives;
- derivatives.resize(symbol_table.endo_nbr());
+ vector<vector<pair<pair<int, int>, int >>> my_derivatives;
+ my_derivatives.resize(symbol_table.endo_nbr());
count_u = symbol_table.endo_nbr();
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
+ int deriv_id = first_derivative.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
expr_t d1 = first_derivative.second;
- unsigned int eq = first_derivative.first.first;
+ unsigned int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
fnumexpr.write(code_file, instruction_number);
- if (!derivatives[eq].size())
- derivatives[eq].clear();
- derivatives[eq].emplace_back(make_pair(var, lag), count_u);
+ if (!my_derivatives[eq].size())
+ my_derivatives[eq].clear();
+ my_derivatives[eq].emplace_back(make_pair(var, lag), count_u);
d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FSTPU_ fstpu(count_u);
@@ -1132,10 +1131,10 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
{
FLDR_ fldr(i);
fldr.write(code_file, instruction_number);
- if (derivatives[i].size())
+ if (my_derivatives[i].size())
{
- for (vector<pair<pair<int, int>, int>>::const_iterator it = derivatives[i].begin();
- it != derivatives[i].end(); it++)
+ for (vector<pair<pair<int, int>, int>>::const_iterator it = my_derivatives[i].begin();
+ it != my_derivatives[i].end(); it++)
{
FLDU_ fldu(it->second);
fldu.write(code_file, instruction_number);
@@ -1143,7 +1142,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
fldv.write(code_file, instruction_number);
FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
fbinary.write(code_file, instruction_number);
- if (it != derivatives[i].begin())
+ if (it != my_derivatives[i].begin())
{
FBINARY_ fbinary{static_cast<int>(BinaryOpcode::plus)};
fbinary.write(code_file, instruction_number);
@@ -1722,7 +1721,7 @@ DynamicModel::writeDynamicCFile(const string &basename, const int order) const
string filename_mex = basename + "/model/src/dynamic_mex.c";
ofstream mDynamicModelFile, mDynamicMexFile;
- int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size();
+ int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
mDynamicModelFile.open(filename, ios::out | ios::binary);
if (!mDynamicModelFile.is_open())
@@ -1835,7 +1834,7 @@ DynamicModel::writeDynamicCFile(const string &basename, const int order) const
<< " if (nlhs >= 3)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v2. */" << endl
- << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3
+ << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3
<< ", mxREAL);" << endl
<< " double *v2 = mxGetPr(plhs[2]);" << endl
<< " dynamic_g2_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
@@ -1845,7 +1844,7 @@ DynamicModel::writeDynamicCFile(const string &basename, const int order) const
<< " if (nlhs >= 4)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v3. */" << endl
- << " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl
+ << " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 3 << ", mxREAL);" << endl
<< " double *v3 = mxGetPr(plhs[3]);" << endl
<< " dynamic_g3_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
<< " dynamic_g3(y, x, nb_row_x, params, steady_state, it_, T, v3);" << endl
@@ -1893,9 +1892,9 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const
void
DynamicModel::setNonZeroHessianEquations(map<int, string> &eqs)
{
- for (const auto &it : second_derivatives)
+ for (const auto &it : derivatives[2])
{
- int eq = get<0>(it.first);
+ int eq = it.first[0];
if (nonzero_hessian_eqs.find(eq) == nonzero_hessian_eqs.end())
{
nonzero_hessian_eqs[eq] = "";
@@ -2460,7 +2459,7 @@ DynamicModel::writeDynamicMatlabCompatLayer(const string &basename) const
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
- int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size();
+ int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
output << "function [residual, g1, g2, g3] = dynamic(y, x, params, steady_state, it_)" << endl
<< " T = NaN(" << ntt << ", 1);" << endl
@@ -2514,11 +2513,11 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
model_tt_output, output_type, tef_terms);
- writeTemporaryTerms(temporary_terms_res,
+ writeTemporaryTerms(temporary_terms_derivatives[0],
temp_term_union,
temporary_terms_idxs,
model_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_res.begin(), temporary_terms_res.end());
+ temp_term_union.insert(temporary_terms_derivatives[0].begin(), temporary_terms_derivatives[0].end());
writeModelEquations(model_output, output_type, temp_term_union);
@@ -2526,18 +2525,18 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
// Writing Jacobian
- if (!first_derivatives.empty())
+ if (!derivatives[1].empty())
{
- writeTemporaryTerms(temporary_terms_g1,
+ writeTemporaryTerms(temporary_terms_derivatives[1],
temp_term_union,
temporary_terms_idxs,
jacobian_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
+ temp_term_union.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
int eq, var;
- tie(eq, var) = first_derivative.first;
+ tie(eq, var) = vectorToTuple<2>(first_derivative.first);
expr_t d1 = first_derivative.second;
jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
@@ -2549,19 +2548,19 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
}
// Writing Hessian
- if (!second_derivatives.empty())
+ if (!derivatives[2].empty())
{
- writeTemporaryTerms(temporary_terms_g2,
+ writeTemporaryTerms(temporary_terms_derivatives[2],
temp_term_union,
temporary_terms_idxs,
hessian_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
+ temp_term_union.insert(temporary_terms_derivatives[2].begin(), temporary_terms_derivatives[2].end());
int k = 0; // Keep the line of a 2nd derivative in v2
- for (const auto & second_derivative : second_derivatives)
+ for (const auto & second_derivative : derivatives[2])
{
int eq, var1, var2;
- tie(eq, var1, var2) = second_derivative.first;
+ tie(eq, var1, var2) = vectorToTuple<3>(second_derivative.first);
expr_t d2 = second_derivative.second;
int id1 = getDynJacobianCol(var1);
@@ -2618,19 +2617,19 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
}
// Writing third derivatives
- if (!third_derivatives.empty())
+ if (!derivatives[3].empty())
{
- writeTemporaryTerms(temporary_terms_g3,
+ writeTemporaryTerms(temporary_terms_derivatives[3],
temp_term_union,
temporary_terms_idxs,
third_derivatives_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
+ temp_term_union.insert(temporary_terms_derivatives[3].begin(), temporary_terms_derivatives[3].end());
int k = 0; // Keep the line of a 3rd derivative in v3
- for (const auto & third_derivative : third_derivatives)
+ for (const auto & third_derivative : derivatives[3])
{
int eq, var1, var2, var3;
- tie(eq, var1, var2, var3) = third_derivative.first;
+ tie(eq, var1, var2, var3) = vectorToTuple<4>(third_derivative.first);
expr_t d3 = third_derivative.second;
int id1 = getDynJacobianCol(var1);
@@ -2714,7 +2713,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
init_output << "residual = zeros(" << nrows << ", 1);";
writeDynamicModelHelper(basename, "dynamic_resid", "residual",
"dynamic_resid_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(),
"", init_output, end_output,
model_output, model_tt_output);
@@ -2723,7 +2722,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
init_output << "g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");";
writeDynamicModelHelper(basename, "dynamic_g1", "g1",
"dynamic_g1_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
"dynamic_resid_tt",
init_output, end_output,
jacobian_output, jacobian_tt_output);
@@ -2731,17 +2730,17 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
init_output.str(string());
init_output.clear();
- if (second_derivatives.size())
+ if (derivatives[2].size())
{
- init_output << "v2 = zeros(" << NNZDerivatives[1] << ",3);";
+ init_output << "v2 = zeros(" << NNZDerivatives[2] << ",3);";
end_output << "g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << nrows << "," << hessianColsNbr << ");";
}
else
init_output << "g2 = sparse([],[],[]," << nrows << "," << hessianColsNbr << ");";
writeDynamicModelHelper(basename, "dynamic_g2", "g2",
"dynamic_g2_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size()
- + temporary_terms_g2.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()
+ + temporary_terms_derivatives[2].size(),
"dynamic_g1_tt",
init_output, end_output,
hessian_output, hessian_tt_output);
@@ -2752,17 +2751,17 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
end_output.str(string());
end_output.clear();
int ncols = hessianColsNbr * dynJacobianColsNbr;
- if (third_derivatives.size())
+ if (derivatives[3].size())
{
- init_output << "v3 = zeros(" << NNZDerivatives[2] << ",3);";
+ init_output << "v3 = zeros(" << NNZDerivatives[3] << ",3);";
end_output << "g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");";
}
else
init_output << "g3 = sparse([],[],[]," << nrows << "," << ncols << ");";
writeDynamicModelHelper(basename, "dynamic_g3", "g3",
"dynamic_g3_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size()
- + temporary_terms_g2.size() + temporary_terms_g3.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()
+ + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size(),
"dynamic_g2_tt",
init_output, end_output,
third_derivatives_output, third_derivatives_tt_output);
@@ -2878,10 +2877,10 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
// Write the number of temporary terms
output << "tmp_nbr = zeros(Int,4)" << endl
- << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_res.size() << "# Number of temporary terms for the residuals" << endl
- << "tmp_nbr[2] = " << temporary_terms_g1.size() << "# Number of temporary terms for g1 (jacobian)" << endl
- << "tmp_nbr[3] = " << temporary_terms_g2.size() << "# Number of temporary terms for g2 (hessian)" << endl
- << "tmp_nbr[4] = " << temporary_terms_g3.size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
+ << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl
+ << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl
+ << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl
+ << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
// dynamicResidTT!
output << "function dynamicResidTT!(T::Vector{Float64}," << endl
@@ -2895,7 +2894,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
output << "function dynamicResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl
<< " y::Vector{Float64}, x::Matrix{Float64}, "
<< "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
- << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_res.size() << endl
+ << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl
<< " @assert length(residual) == " << nrows << endl
<< " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
@@ -2920,7 +2919,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
<< " y::Vector{Float64}, x::Matrix{Float64}, "
<< "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
<< " @assert length(T) >= "
- << temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() << endl
+ << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
<< " @assert size(g1) == (" << nrows << ", " << dynJacobianColsNbr << ")" << endl
<< " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
@@ -2945,7 +2944,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
output << "function dynamicG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
<< " y::Vector{Float64}, x::Matrix{Float64}, "
<< "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
- << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() << endl
+ << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl
<< " @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl
<< " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
@@ -2972,7 +2971,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
<< " y::Vector{Float64}, x::Matrix{Float64}, "
<< "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
<< " @assert length(T) >= "
- << temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size() << endl
+ << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl
<< " @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
<< " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
@@ -3117,10 +3116,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
<< modstruct << "ndynamic = " << npred+nboth+nfwrd << ";" << endl;
if (!julia)
output << modstruct << "dynamic_tmp_nbr = zeros(4,1); % Number of temporaries used for the dynamic model" << endl
- << modstruct << "dynamic_tmp_nbr(1) = " << temporary_terms_res.size() << "; % Number of temporaries used for the evaluation of the residuals" << endl
- << modstruct << "dynamic_tmp_nbr(2) = " << temporary_terms_g1.size() << "; % Number of temporaries used for the evaluation of g1 (jacobian)" << endl
- << modstruct << "dynamic_tmp_nbr(3) = " << temporary_terms_g2.size() << "; % Number of temporaries used for the evaluation of g2 (hessian)" << endl
- << modstruct << "dynamic_tmp_nbr(4) = " << temporary_terms_g3.size() << "; % Number of temporaries used for the evaluation of g3 (third order derivatives)" << endl;
+ << modstruct << "dynamic_tmp_nbr(1) = " << temporary_terms_derivatives[0].size() << "; % Number of temporaries used for the evaluation of the residuals" << endl
+ << modstruct << "dynamic_tmp_nbr(2) = " << temporary_terms_derivatives[1].size() << "; % Number of temporaries used for the evaluation of g1 (jacobian)" << endl
+ << modstruct << "dynamic_tmp_nbr(3) = " << temporary_terms_derivatives[2].size() << "; % Number of temporaries used for the evaluation of g2 (hessian)" << endl
+ << modstruct << "dynamic_tmp_nbr(4) = " << temporary_terms_derivatives[3].size() << "; % Number of temporaries used for the evaluation of g3 (third order derivatives)" << endl;
// Write equation tags
if (julia)
@@ -3413,12 +3412,12 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
state_equ.push_back(equation_reordered[variable_inv_reordered[*it - 1]]+1);
map<pair< int, pair<int, int>>, int> lag_row_incidence;
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
+ int deriv_id = first_derivative.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
- int eq = first_derivative.first.first;
+ int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
@@ -3625,14 +3624,14 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
// Write number of non-zero derivatives
// Use -1 if the derivatives have not been computed
output << modstruct << (julia ? "nnzderivatives" : "NNZDerivatives")
- << " = [" << NNZDerivatives[0] << "; ";
+ << " = [" << NNZDerivatives[1] << "; ";
if (order > 1)
- output << NNZDerivatives[1] << "; ";
+ output << NNZDerivatives[2] << "; ";
else
output << "-1; ";
if (order > 2)
- output << NNZDerivatives[2];
+ output << NNZDerivatives[3];
else
output << "-1";
output << "];" << endl;
@@ -3646,13 +3645,13 @@ map<pair<int, pair<int, int >>, expr_t>
DynamicModel::collect_first_order_derivatives_endogenous()
{
map<pair<int, pair<int, int >>, expr_t> endo_derivatives;
- for (auto & first_derivative : first_derivatives)
+ for (auto & first_derivative : derivatives[1])
{
- if (getTypeByDerivID(first_derivative.first.second) == SymbolType::endogenous)
+ if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous)
{
- int eq = first_derivative.first.first;
- int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
- int lag = getLagByDerivID(first_derivative.first.second);
+ int eq = first_derivative.first[0];
+ int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
+ int lag = getLagByDerivID(first_derivative.first[1]);
endo_derivatives[{ eq, { var, lag } }] = first_derivative.second;
}
}
@@ -4718,7 +4717,7 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivat
int eqr = it_l.second.first;
int varr = it_l.second.second;
if (Deriv_type == 0)
- first_chain_rule_derivatives[{ eqr, { varr, lag } }] = first_derivatives[{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
+ first_chain_rule_derivatives[{ eqr, { varr, lag } }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
else if (Deriv_type == 1)
first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
@@ -4763,16 +4762,16 @@ DynamicModel::collect_block_first_order_derivatives()
exo_max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
exo_det_max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
- for (auto & first_derivative : first_derivatives)
+ for (auto & first_derivative : derivatives[1])
{
- int eq = first_derivative.first.first;
- int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
- int lag = getLagByDerivID(first_derivative.first.second);
+ int eq = first_derivative.first[0];
+ int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
+ int lag = getLagByDerivID(first_derivative.first[1]);
int block_eq = equation_2_block[eq];
int block_var = 0;
derivative_t tmp_derivative;
lag_var_t lag_var;
- switch (getTypeByDerivID(first_derivative.first.second))
+ switch (getTypeByDerivID(first_derivative.first[1]))
{
case SymbolType::endogenous:
block_var = variable_2_block[var];
@@ -4783,7 +4782,7 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > endo_max_leadlag_block[block_eq].second)
endo_max_leadlag_block[block_eq] = { endo_max_leadlag_block[block_eq].first, lag };
tmp_derivative = derivative_endo[block_eq];
- tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
+ tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
derivative_endo[block_eq] = tmp_derivative;
}
else
@@ -4807,7 +4806,7 @@ DynamicModel::collect_block_first_order_derivatives()
}
}
}
- tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
+ tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
derivative_other_endo[block_eq] = tmp_derivative;
lag_var = other_endo_block[block_eq];
if (lag_var.find(lag) == lag_var.end())
@@ -4836,7 +4835,7 @@ DynamicModel::collect_block_first_order_derivatives()
}
}
}
- tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
+ tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
derivative_exo[block_eq] = tmp_derivative;
lag_var = exo_block[block_eq];
if (lag_var.find(lag) == lag_var.end())
@@ -4864,7 +4863,7 @@ DynamicModel::collect_block_first_order_derivatives()
}
}
}
- tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
+ tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
derivative_exo_det[block_eq] = tmp_derivative;
lag_var = exo_det_block[block_eq];
if (lag_var.find(lag) == lag_var.end())
@@ -5378,11 +5377,7 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context
void
DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) const
{
- if (!residuals_params_derivatives.size()
- && !residuals_params_second_derivatives.size()
- && !jacobian_params_derivatives.size()
- && !jacobian_params_second_derivatives.size()
- && !hessian_params_derivatives.size())
+ if (!params_derivatives.size())
return;
ExprNodeOutputType output_type = (julia ? ExprNodeOutputType::juliaDynamicModel : ExprNodeOutputType::matlabDynamicModel);
@@ -5398,10 +5393,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
writeTemporaryTerms(params_derivs_temporary_terms, {}, params_derivs_temporary_terms_idxs, model_output, output_type, tef_terms);
- for (const auto & residuals_params_derivative : residuals_params_derivatives)
+ for (const auto & residuals_params_derivative : params_derivatives.find({ 0, 1 })->second)
{
int eq, param;
- tie(eq, param) = residuals_params_derivative.first;
+ tie(eq, param) = vectorToTuple<2>(residuals_params_derivative.first);
expr_t d1 = residuals_params_derivative.second;
int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
@@ -5412,10 +5407,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
jacobian_output << ";" << endl;
}
- for (const auto & jacobian_params_derivative : jacobian_params_derivatives)
+ for (const auto & jacobian_params_derivative : params_derivatives.find({ 1, 1 })->second)
{
int eq, var, param;
- tie(eq, var, param) = jacobian_params_derivative.first;
+ tie(eq, var, param) = vectorToTuple<3>(jacobian_params_derivative.first);
expr_t d2 = jacobian_params_derivative.second;
int var_col = getDynJacobianCol(var) + 1;
@@ -5428,10 +5423,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
}
int i = 1;
- for (const auto &it : residuals_params_second_derivatives)
+ for (const auto &it : params_derivatives.find({ 0, 2 })->second)
{
int eq, param1, param2;
- tie(eq, param1, param2) = it.first;
+ tie(eq, param1, param2) = vectorToTuple<3>(it.first);
expr_t d2 = it.second;
int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
@@ -5452,10 +5447,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
}
i = 1;
- for (const auto &it : jacobian_params_second_derivatives)
+ for (const auto &it : params_derivatives.find({ 1, 2 })->second)
{
int eq, var, param1, param2;
- tie(eq, var, param1, param2) = it.first;
+ tie(eq, var, param1, param2) = vectorToTuple<4>(it.first);
expr_t d2 = it.second;
int var_col = getDynJacobianCol(var) + 1;
@@ -5479,10 +5474,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
}
i = 1;
- for (const auto &it : hessian_params_derivatives)
+ for (const auto &it : params_derivatives.find({ 2, 1 })->second)
{
int eq, var1, var2, param;
- tie(eq, var1, var2, param) = it.first;
+ tie(eq, var1, var2, param) = vectorToTuple<4>(it.first);
expr_t d2 = it.second;
int var1_col = getDynJacobianCol(var1) + 1;
@@ -5579,13 +5574,13 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
<< "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
<< "if nargout >= 3" << endl
- << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
+ << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
<< hessian1_output.str()
- << "gpp = zeros(" << jacobian_params_second_derivatives.size() << ",5);" << endl
+ << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
<< third_derivs_output.str()
<< "end" << endl
<< "if nargout >= 5" << endl
- << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl
+ << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
<< third_derivs1_output.str()
<< "end" << endl
<< "end" << endl;
@@ -5606,11 +5601,11 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
<< jacobian_output.str()
<< "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
- << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
+ << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
<< hessian1_output.str()
- << "gpp = zeros(" << jacobian_params_second_derivatives.size() << ",5);" << endl
+ << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
<< third_derivs_output.str()
- << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl
+ << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
<< third_derivs1_output.str()
<< "(rp, gp, rpp, gpp, hp)" << endl
<< "end" << endl
@@ -6371,7 +6366,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
deriv_node_temp_terms_t tef_terms;
temporary_terms_t temp_term_empty;
- temporary_terms_t temp_term_union = temporary_terms_res;
+ temporary_terms_t temp_term_union = temporary_terms_derivatives[0];
temporary_terms_t temp_term_union_m_1;
string concat = "";
@@ -6379,27 +6374,27 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
writeJsonModelLocalVariables(model_local_vars_output, tef_terms);
- writeJsonTemporaryTerms(temporary_terms_res, temp_term_union_m_1, model_output, tef_terms, concat);
+ writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union_m_1, model_output, tef_terms, concat);
model_output << ", ";
writeJsonModelEquations(model_output, true);
// Writing Jacobian
temp_term_union_m_1 = temp_term_union;
- temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
+ temp_term_union.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
concat = "jacobian";
writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, tef_terms, concat);
jacobian_output << ", \"jacobian\": {"
<< " \"nrows\": " << equations.size()
<< ", \"ncols\": " << dynJacobianColsNbr
<< ", \"entries\": [";
- for (auto it = first_derivatives.begin();
- it != first_derivatives.end(); it++)
+ for (auto it = derivatives[1].begin();
+ it != derivatives[1].end(); it++)
{
- if (it != first_derivatives.begin())
+ if (it != derivatives[1].begin())
jacobian_output << ", ";
int eq, var;
- tie(eq, var) = it->first;
+ tie(eq, var) = vectorToTuple<2>(it->first);
int col = getDynJacobianCol(var);
expr_t d1 = it->second;
@@ -6422,21 +6417,21 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
// Writing Hessian
temp_term_union_m_1 = temp_term_union;
- temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
+ temp_term_union.insert(temporary_terms_derivatives[2].begin(), temporary_terms_derivatives[2].end());
concat = "hessian";
writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, hessian_output, tef_terms, concat);
hessian_output << ", \"hessian\": {"
<< " \"nrows\": " << equations.size()
<< ", \"ncols\": " << hessianColsNbr
<< ", \"entries\": [";
- for (auto it = second_derivatives.begin();
- it != second_derivatives.end(); it++)
+ for (auto it = derivatives[2].begin();
+ it != derivatives[2].end(); it++)
{
- if (it != second_derivatives.begin())
+ if (it != derivatives[2].begin())
hessian_output << ", ";
int eq, var1, var2;
- tie(eq, var1, var2) = it->first;
+ tie(eq, var1, var2) = vectorToTuple<3>(it->first);
expr_t d2 = it->second;
int id1 = getDynJacobianCol(var1);
int id2 = getDynJacobianCol(var2);
@@ -6467,21 +6462,21 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
// Writing third derivatives
temp_term_union_m_1 = temp_term_union;
- temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
+ temp_term_union.insert(temporary_terms_derivatives[3].begin(), temporary_terms_derivatives[3].end());
concat = "third_derivatives";
writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, third_derivatives_output, tef_terms, concat);
third_derivatives_output << ", \"third_derivative\": {"
<< " \"nrows\": " << equations.size()
<< ", \"ncols\": " << hessianColsNbr * dynJacobianColsNbr
<< ", \"entries\": [";
- for (auto it = third_derivatives.begin();
- it != third_derivatives.end(); it++)
+ for (auto it = derivatives[3].begin();
+ it != derivatives[3].end(); it++)
{
- if (it != third_derivatives.begin())
+ if (it != derivatives[3].begin())
third_derivatives_output << ", ";
int eq, var1, var2, var3;
- tie(eq, var1, var2, var3) = it->first;
+ tie(eq, var1, var2, var3) = vectorToTuple<4>(it->first);
expr_t d3 = it->second;
if (writeDetails)
@@ -6538,11 +6533,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
void
DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const
{
- if (!residuals_params_derivatives.size()
- && !residuals_params_second_derivatives.size()
- && !jacobian_params_derivatives.size()
- && !jacobian_params_second_derivatives.size()
- && !hessian_params_derivatives.size())
+ if (!params_derivatives.size())
return;
ostringstream model_local_vars_output; // Used for storing model local vars
@@ -6563,14 +6554,14 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
<< " \"neqs\": " << equations.size()
<< ", \"nparamcols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = residuals_params_derivatives.begin();
- it != residuals_params_derivatives.end(); it++)
+ auto &rp = params_derivatives.find({ 0, 1 })->second;
+ for (auto it = rp.begin(); it != rp.end(); it++)
{
- if (it != residuals_params_derivatives.begin())
+ if (it != rp.begin())
jacobian_output << ", ";
int eq, param;
- tie(eq, param) = it->first;
+ tie(eq, param) = vectorToTuple<2>(it->first);
expr_t d1 = it->second;
int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
@@ -6590,19 +6581,20 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
jacobian_output << "\"}" << endl;
}
jacobian_output << "]}";
+
hessian_output << "\"deriv_jacobian_wrt_params\": {"
<< " \"neqs\": " << equations.size()
<< ", \"nvarcols\": " << dynJacobianColsNbr
<< ", \"nparamcols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = jacobian_params_derivatives.begin();
- it != jacobian_params_derivatives.end(); it++)
+ auto &gp = params_derivatives.find({ 1, 1 })->second;
+ for (auto it = gp.begin(); it != gp.end(); it++)
{
- if (it != jacobian_params_derivatives.begin())
+ if (it != gp.begin())
hessian_output << ", ";
int eq, var, param;
- tie(eq, var, param) = it->first;
+ tie(eq, var, param) = vectorToTuple<3>(it->first);
expr_t d2 = it->second;
int var_col = getDynJacobianCol(var) + 1;
@@ -6632,14 +6624,14 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
<< ", \"nparam1cols\": " << symbol_table.param_nbr()
<< ", \"nparam2cols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = residuals_params_second_derivatives.begin();
- it != residuals_params_second_derivatives.end(); ++it)
+ auto &rpp = params_derivatives.find({ 0, 2 })->second;
+ for (auto it = rpp.begin(); it != rpp.end(); ++it)
{
- if (it != residuals_params_second_derivatives.begin())
+ if (it != rpp.begin())
hessian1_output << ", ";
int eq, param1, param2;
- tie(eq, param1, param2) = it->first;
+ tie(eq, param1, param2) = vectorToTuple<3>(it->first);
expr_t d2 = it->second;
int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
@@ -6661,20 +6653,21 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
hessian1_output << "\"}" << endl;
}
hessian1_output << "]}";
+
third_derivs_output << "\"second_deriv_jacobian_wrt_params\": {"
<< " \"neqs\": " << equations.size()
<< ", \"nvarcols\": " << dynJacobianColsNbr
<< ", \"nparam1cols\": " << symbol_table.param_nbr()
<< ", \"nparam2cols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = jacobian_params_second_derivatives.begin();
- it != jacobian_params_second_derivatives.end(); ++it)
+ auto &gpp = params_derivatives.find({ 1, 2 })->second;
+ for (auto it = gpp.begin(); it != gpp.end(); ++it)
{
- if (it != jacobian_params_second_derivatives.begin())
+ if (it != gpp.begin())
third_derivs_output << ", ";
int eq, var, param1, param2;
- tie(eq, var, param1, param2) = it->first;
+ tie(eq, var, param1, param2) = vectorToTuple<4>(it->first);
expr_t d2 = it->second;
int var_col = getDynJacobianCol(var) + 1;
@@ -6708,14 +6701,14 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
<< ", \"nvar2cols\": " << dynJacobianColsNbr
<< ", \"nparamcols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = hessian_params_derivatives.begin();
- it != hessian_params_derivatives.end(); ++it)
+ auto &hp = params_derivatives.find({ 2, 1 })->second;
+ for (auto it = hp.begin(); it != hp.end(); ++it)
{
- if (it != hessian_params_derivatives.begin())
+ if (it != hp.begin())
third_derivs1_output << ", ";
int eq, var1, var2, param;
- tie(eq, var1, var2, param) = it->first;
+ tie(eq, var1, var2, param) = vectorToTuple<4>(it->first);
expr_t d2 = it->second;
int var1_col = getDynJacobianCol(var1) + 1;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 4e3d2959..76d878eb 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -43,51 +43,41 @@ ModelTree::copyHelper(const ModelTree &m)
for (const auto & it : m.aux_equations)
aux_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
+ auto convert_deriv_map = [f](map<vector<int>, expr_t> dm)
+ {
+ map<vector<int>, expr_t> dm2;
+ for (const auto &it : dm)
+ dm2.emplace(it.first, f(it.second));
+ return dm2;
+ };
+
// Derivatives
- for (const auto & it : m.first_derivatives)
- first_derivatives[it.first] = f(it.second);
- for (const auto & it : m.second_derivatives)
- second_derivatives[it.first] = f(it.second);
- for (const auto & it : m.third_derivatives)
- third_derivatives[it.first] = f(it.second);
- for (const auto & it : m.residuals_params_derivatives)
- residuals_params_derivatives[it.first] = f(it.second);
- for (const auto & it : m.residuals_params_second_derivatives)
- residuals_params_second_derivatives[it.first] = f(it.second);
- for (const auto & it : m.jacobian_params_derivatives)
- jacobian_params_derivatives[it.first] = f(it.second);
- for (const auto & it : m.jacobian_params_second_derivatives)
- jacobian_params_second_derivatives[it.first] = f(it.second);
- for (const auto & it : m.hessian_params_derivatives)
- hessian_params_derivatives[it.first] = f(it.second);
+ for (const auto &it : m.derivatives)
+ derivatives.push_back(convert_deriv_map(it));
+ for (const auto &it : m.params_derivatives)
+ params_derivatives[it.first] = convert_deriv_map(it.second);
+
+ auto convert_temporary_terms_t = [f](temporary_terms_t tt)
+ {
+ temporary_terms_t tt2;
+ for (const auto &it : tt)
+ tt2.insert(f(it));
+ return tt2;
+ };
// Temporary terms
for (const auto & it : m.temporary_terms)
temporary_terms.insert(f(it));
for (const auto & it : m.temporary_terms_mlv)
temporary_terms_mlv[f(it.first)] = f(it.second);
- for (const auto & it : m.temporary_terms_res)
- temporary_terms_res.insert(f(it));
- for (const auto & it : m.temporary_terms_g1)
- temporary_terms_g1.insert(f(it));
- for (const auto & it : m.temporary_terms_g2)
- temporary_terms_g2.insert(f(it));
- for (const auto & it : m.temporary_terms_g3)
- temporary_terms_g3.insert(f(it));
+ for (const auto &it : m.temporary_terms_derivatives)
+ temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
for (const auto & it : m.temporary_terms_idxs)
temporary_terms_idxs[f(it.first)] = it.second;
for (const auto & it : m.params_derivs_temporary_terms)
params_derivs_temporary_terms.insert(f(it));
- for (const auto & it : m.params_derivs_temporary_terms_res)
- params_derivs_temporary_terms_res.insert(f(it));
- for (const auto & it : m.params_derivs_temporary_terms_g1)
- params_derivs_temporary_terms_g1.insert(f(it));
- for (const auto & it : m.params_derivs_temporary_terms_res2)
- params_derivs_temporary_terms_res2.insert(f(it));
- for (const auto & it : m.params_derivs_temporary_terms_g12)
- params_derivs_temporary_terms_g12.insert(f(it));
- for (const auto & it : m.params_derivs_temporary_terms_g2)
- params_derivs_temporary_terms_g2.insert(f(it));
+ for (const auto & it : m.params_derivs_temporary_terms_split)
+ params_derivs_temporary_terms_split[it.first] = convert_temporary_terms_t(it.second);
for (const auto & it : m.params_derivs_temporary_terms_idxs)
params_derivs_temporary_terms_idxs[f(it.first)] = it.second;
@@ -129,28 +119,14 @@ ModelTree::operator=(const ModelTree &m)
equation_tags = m.equation_tags;
NNZDerivatives = m.NNZDerivatives;
- first_derivatives.clear();
- second_derivatives.clear();
- third_derivatives.clear();
- residuals_params_derivatives.clear();
- residuals_params_second_derivatives.clear();
- jacobian_params_derivatives.clear();
- jacobian_params_second_derivatives.clear();
- hessian_params_derivatives.clear();
+ derivatives.clear();
+ params_derivatives.clear();
temporary_terms.clear();
temporary_terms_mlv.clear();
- temporary_terms_res.clear();
- temporary_terms_g1.clear();
- temporary_terms_g2.clear();
- temporary_terms_g3.clear();
- temporary_terms_idxs.clear();
+ temporary_terms_derivatives.clear();
params_derivs_temporary_terms.clear();
- params_derivs_temporary_terms_res.clear();
- params_derivs_temporary_terms_g1.clear();
- params_derivs_temporary_terms_res2.clear();
- params_derivs_temporary_terms_g12.clear();
- params_derivs_temporary_terms_g2.clear();
+ params_derivs_temporary_terms_split.clear();
params_derivs_temporary_terms_idxs.clear();
trend_symbols_map.clear();
@@ -350,8 +326,8 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
contemporaneous_jacobian[{ it->first.first, it->first.second }] = 0;
try
{
- if (first_derivatives.find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == first_derivatives.end())
- first_derivatives[{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
+ if (derivatives[1].find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == derivatives[1].end())
+ derivatives[1][{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
}
catch (DataTree::UnknownDerivIDException &e)
{
@@ -397,15 +373,15 @@ void
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose)
{
int nb_elements_contemparenous_Jacobian = 0;
- set<pair<int, int>> jacobian_elements_to_delete;
- for (first_derivatives_t::const_iterator it = first_derivatives.begin();
- it != first_derivatives.end(); it++)
+ set<vector<int>> jacobian_elements_to_delete;
+ for (auto it = derivatives[1].begin();
+ it != derivatives[1].end(); it++)
{
- int deriv_id = it->first.second;
+ int deriv_id = it->first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
expr_t Id = it->second;
- int eq = it->first.first;
+ int eq = it->first[0];
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
@@ -429,7 +405,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
{
if (verbose)
cout << "the coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")" << endl;
- jacobian_elements_to_delete.emplace(eq, deriv_id);
+ jacobian_elements_to_delete.insert({ eq, deriv_id });
}
else
{
@@ -449,11 +425,11 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
// Get rid of the elements of the Jacobian matrix below the cutoff
for (const auto & it : jacobian_elements_to_delete)
- first_derivatives.erase(it);
+ derivatives[1].erase(it);
if (jacobian_elements_to_delete.size() > 0)
{
- cout << jacobian_elements_to_delete.size() << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
+ cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
<< "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
}
}
@@ -1277,10 +1253,11 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
ExternalFunctionsTable &external_functions_table_arg,
bool is_dynamic_arg) :
- DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg}
+ DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
+ derivatives(4),
+ NNZDerivatives(4, 0),
+ temporary_terms_derivatives(4)
{
- for (int & NNZDerivative : NNZDerivatives)
- NNZDerivative = 0;
}
int
@@ -1294,8 +1271,8 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
ExprNodeOutputType output_type,
const temporary_terms_t &temporary_terms) const
{
- auto it = first_derivatives.find({ eq, getDerivID(symb_id, lag) });
- if (it != first_derivatives.end())
+ auto it = derivatives[1].find({ eq, getDerivID(symb_id, lag) });
+ if (it != derivatives[1].end())
(it->second)->writeOutput(output, output_type, temporary_terms, {});
else
output << 0;
@@ -1311,8 +1288,8 @@ ModelTree::computeJacobian(const set<int> &vars)
expr_t d1 = equations[eq]->getDerivative(var);
if (d1 == Zero)
continue;
- first_derivatives[{ eq, var }] = d1;
- ++NNZDerivatives[0];
+ derivatives[1][{ eq, var }] = d1;
+ ++NNZDerivatives[1];
}
}
}
@@ -1320,10 +1297,10 @@ ModelTree::computeJacobian(const set<int> &vars)
void
ModelTree::computeHessian(const set<int> &vars)
{
- for (const auto &it : first_derivatives)
+ for (const auto &it : derivatives[1])
{
int eq, var1;
- tie(eq, var1) = it.first;
+ tie(eq, var1) = vectorToTuple<2>(it.first);
expr_t d1 = it.second;
// Store only second derivatives with var2 <= var1
@@ -1335,11 +1312,11 @@ ModelTree::computeHessian(const set<int> &vars)
expr_t d2 = d1->getDerivative(var2);
if (d2 == Zero)
continue;
- second_derivatives[{ eq, var1, var2 }] = d2;
+ derivatives[2][{ eq, var1, var2 }] = d2;
if (var2 == var1)
- ++NNZDerivatives[1];
+ ++NNZDerivatives[2];
else
- NNZDerivatives[1] += 2;
+ NNZDerivatives[2] += 2;
}
}
}
@@ -1347,10 +1324,10 @@ ModelTree::computeHessian(const set<int> &vars)
void
ModelTree::computeThirdDerivatives(const set<int> &vars)
{
- for (const auto &it : second_derivatives)
+ for (const auto &it : derivatives[2])
{
int eq, var1, var2;
- tie(eq, var1, var2) = it.first;
+ tie(eq, var1, var2) = vectorToTuple<3>(it.first);
// By construction, var2 <= var1
expr_t d2 = it.second;
@@ -1364,13 +1341,13 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
expr_t d3 = d2->getDerivative(var3);
if (d3 == Zero)
continue;
- third_derivatives[{ eq, var1, var2, var3 }] = d3;
+ derivatives[3][{ eq, var1, var2, var3 }] = d3;
if (var3 == var2 && var2 == var1)
- ++NNZDerivatives[2];
+ ++NNZDerivatives[3];
else if (var3 == var2 || var2 == var1)
- NNZDerivatives[2] += 3;
+ NNZDerivatives[3] += 3;
else
- NNZDerivatives[2] += 6;
+ NNZDerivatives[3] += 6;
}
}
}
@@ -1381,10 +1358,8 @@ 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_res.clear();
- temporary_terms_g1.clear();
- temporary_terms_g2.clear();
- temporary_terms_g3.clear();
+ temporary_terms_derivatives.clear();
+ temporary_terms_derivatives.resize(4);
// Collect all model local variables appearing in equations. See #101
// All used model local variables are automatically set as temporary variables
@@ -1400,10 +1375,10 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
}
map<NodeTreeReference, temporary_terms_t> temp_terms_map;
- temp_terms_map[NodeTreeReference::residuals] = temporary_terms_res;
- temp_terms_map[NodeTreeReference::firstDeriv] = temporary_terms_g1;
- temp_terms_map[NodeTreeReference::secondDeriv] = temporary_terms_g2;
- temp_terms_map[NodeTreeReference::thirdDeriv] = temporary_terms_g3;
+ 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];
if (!no_tmp_terms)
{
@@ -1412,17 +1387,17 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
temp_terms_map,
is_matlab, NodeTreeReference::residuals);
- for (auto & first_derivative : first_derivatives)
+ for (auto & first_derivative : derivatives[1])
first_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
is_matlab, NodeTreeReference::firstDeriv);
- for (auto & second_derivative : second_derivatives)
+ for (auto & second_derivative : derivatives[2])
second_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
is_matlab, NodeTreeReference::secondDeriv);
- for (auto & third_derivative : third_derivatives)
+ for (auto & third_derivative : derivatives[3])
third_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
is_matlab, NodeTreeReference::thirdDeriv);
@@ -1432,26 +1407,26 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
it != temp_terms_map.end(); it++)
temporary_terms.insert(it->second.begin(), it->second.end());
- temporary_terms_res = temp_terms_map[NodeTreeReference::residuals];
- temporary_terms_g1 = temp_terms_map[NodeTreeReference::firstDeriv];
- temporary_terms_g2 = temp_terms_map[NodeTreeReference::secondDeriv];
- temporary_terms_g3 = temp_terms_map[NodeTreeReference::thirdDeriv];
+ 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];
int idx = 0;
for (map<expr_t, expr_t, ExprNodeLess>::const_iterator it = temporary_terms_mlv.begin();
it != temporary_terms_mlv.end(); it++)
temporary_terms_idxs[it->first] = idx++;
- for (auto it : temporary_terms_res)
+ for (auto it : temporary_terms_derivatives[0])
temporary_terms_idxs[it] = idx++;
- for (auto it : temporary_terms_g1)
+ for (auto it : temporary_terms_derivatives[1])
temporary_terms_idxs[it] = idx++;
- for (auto it : temporary_terms_g2)
+ for (auto it : temporary_terms_derivatives[2])
temporary_terms_idxs[it] = idx++;
- for (auto it : temporary_terms_g3)
+ for (auto it : temporary_terms_derivatives[3])
temporary_terms_idxs[it] = idx++;
}
@@ -1887,12 +1862,12 @@ ModelTree::Write_Inf_To_Bin_File(const string &filename,
exit(EXIT_FAILURE);
}
u_count_int = 0;
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
+ int deriv_id = first_derivative.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
- int eq = first_derivative.first.first;
+ int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
@@ -2084,7 +2059,7 @@ ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, Expr
if (isMatlabOutput(output_type) || isJuliaOutput(output_type))
output << row_nb + 1 << "," << col_nb + 1;
else
- output << row_nb + col_nb * NNZDerivatives[order-1];
+ output << row_nb + col_nb * NNZDerivatives[order];
output << RIGHT_ARRAY_SUBSCRIPT(output_type);
}
@@ -2103,58 +2078,58 @@ ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
expr_t d1 = equations[eq]->getDerivative(param);
if (d1 == Zero)
continue;
- residuals_params_derivatives[{ eq, param }] = d1;
+ params_derivatives[{ 0, 1 }][{ eq, param }] = d1;
}
if (paramsDerivsOrder == 2)
- for (const auto &it : residuals_params_derivatives)
+ for (const auto &it : params_derivatives[{ 0, 1 }])
{
int eq, param1;
- tie(eq, param1) = it.first;
+ tie(eq, param1) = vectorToTuple<2>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
- residuals_params_second_derivatives[{ eq, param1, param }] = d2;
+ params_derivatives[{ 0, 2 }][{ eq, param1, param }] = d2;
}
- for (const auto &it : first_derivatives)
+ for (const auto &it : derivatives[1])
{
int eq, var;
- tie(eq, var) = it.first;
+ tie(eq, var) = vectorToTuple<2>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
- jacobian_params_derivatives[{ eq, var, param }] = d2;
+ params_derivatives[{ 1, 1 }][{ eq, var, param }] = d2;
}
if (paramsDerivsOrder == 2)
{
- for (const auto &it : jacobian_params_derivatives)
+ for (const auto &it : params_derivatives[{ 1, 1 }])
{
int eq, var, param1;
- tie(eq, var, param1) = it.first;
+ tie(eq, var, param1) = vectorToTuple<3>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
- jacobian_params_second_derivatives[{ eq, var, param1, param }] = d2;
+ params_derivatives[{ 1, 2 }][{ eq, var, param1, param }] = d2;
}
- for (const auto &it : second_derivatives)
+ for (const auto &it : derivatives[2])
{
int eq, var1, var2;
- tie(eq, var1, var2) = it.first;
+ tie(eq, var1, var2) = vectorToTuple<3>(it.first);
expr_t d1 = it.second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
- hessian_params_derivatives[{ eq, var1, var2, param }] = d2;
+ params_derivatives[{ 2, 1 }][{ eq, var1, var2, param }] = d2;
}
}
}
@@ -2166,64 +2141,61 @@ ModelTree::computeParamsDerivativesTemporaryTerms()
map<expr_t, pair<int, NodeTreeReference >> reference_count;
params_derivs_temporary_terms.clear();
map<NodeTreeReference, temporary_terms_t> temp_terms_map;
- temp_terms_map[NodeTreeReference::residualsParamsDeriv] = params_derivs_temporary_terms_res;
- temp_terms_map[NodeTreeReference::jacobianParamsDeriv] = params_derivs_temporary_terms_g1;
- temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv] = params_derivs_temporary_terms_res2;
- temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv] = params_derivs_temporary_terms_g12;
- temp_terms_map[NodeTreeReference::hessianParamsDeriv] = params_derivs_temporary_terms_g2;
+ temp_terms_map[NodeTreeReference::residualsParamsDeriv] = params_derivs_temporary_terms_split[{ 0, 1 }];
+ temp_terms_map[NodeTreeReference::jacobianParamsDeriv] = params_derivs_temporary_terms_split[{ 1, 1 }];
+ temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv] = params_derivs_temporary_terms_split[{ 0, 2 }];
+ temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv] = params_derivs_temporary_terms_split[{ 1, 2 }];
+ temp_terms_map[NodeTreeReference::hessianParamsDeriv] = params_derivs_temporary_terms_split[{ 2, 1}];
- for (auto & residuals_params_derivative : residuals_params_derivatives)
+ for (const auto &residuals_params_derivative : params_derivatives[{ 0, 1 }])
residuals_params_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::residualsParamsDeriv);
- for (auto & jacobian_params_derivative : jacobian_params_derivatives)
+ for (const auto &jacobian_params_derivative : params_derivatives[{ 1, 1 }])
jacobian_params_derivative.second->computeTemporaryTerms(reference_count,
temp_terms_map,
true, NodeTreeReference::jacobianParamsDeriv);
- for (second_derivatives_t::const_iterator it = residuals_params_second_derivatives.begin();
- it != residuals_params_second_derivatives.end(); ++it)
- it->second->computeTemporaryTerms(reference_count,
- temp_terms_map,
- true, NodeTreeReference::residualsParamsSecondDeriv);
+ for (const auto &it : params_derivatives[{ 0, 2 }])
+ it.second->computeTemporaryTerms(reference_count,
+ temp_terms_map,
+ true, NodeTreeReference::residualsParamsSecondDeriv);
- for (third_derivatives_t::const_iterator it = jacobian_params_second_derivatives.begin();
- it != jacobian_params_second_derivatives.end(); ++it)
- it->second->computeTemporaryTerms(reference_count,
- temp_terms_map,
- true, NodeTreeReference::jacobianParamsSecondDeriv);
+ for (const auto &it : params_derivatives[{ 1, 2 }])
+ it.second->computeTemporaryTerms(reference_count,
+ temp_terms_map,
+ true, NodeTreeReference::jacobianParamsSecondDeriv);
- for (third_derivatives_t::const_iterator it = hessian_params_derivatives.begin();
- it != hessian_params_derivatives.end(); ++it)
- it->second->computeTemporaryTerms(reference_count,
- temp_terms_map,
- true, NodeTreeReference::hessianParamsDeriv);
+ for (const auto &it : params_derivatives[{ 2, 1 }])
+ it.second->computeTemporaryTerms(reference_count,
+ temp_terms_map,
+ true, NodeTreeReference::hessianParamsDeriv);
for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
it != temp_terms_map.end(); it++)
params_derivs_temporary_terms.insert(it->second.begin(), it->second.end());
- params_derivs_temporary_terms_res = temp_terms_map[NodeTreeReference::residualsParamsDeriv];
- params_derivs_temporary_terms_g1 = temp_terms_map[NodeTreeReference::jacobianParamsDeriv];
- params_derivs_temporary_terms_res2 = temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv];
- params_derivs_temporary_terms_g12 = temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv];
- params_derivs_temporary_terms_g2 = temp_terms_map[NodeTreeReference::hessianParamsDeriv];
+ params_derivs_temporary_terms_split[{ 0, 1 }] = temp_terms_map[NodeTreeReference::residualsParamsDeriv];
+ params_derivs_temporary_terms_split[{ 1, 1 }] = temp_terms_map[NodeTreeReference::jacobianParamsDeriv];
+ params_derivs_temporary_terms_split[{ 0, 2 }] = temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv];
+ params_derivs_temporary_terms_split[{ 1, 2 }] = temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv];
+ params_derivs_temporary_terms_split[{ 2, 1 }] = temp_terms_map[NodeTreeReference::hessianParamsDeriv];
int idx = 0;
- for (auto tt : params_derivs_temporary_terms_res)
+ for (auto tt : params_derivs_temporary_terms_split[{ 0, 1 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
- for (auto tt : params_derivs_temporary_terms_g1)
+ for (auto tt : params_derivs_temporary_terms_split[{ 1, 1 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
- for (auto tt : params_derivs_temporary_terms_res2)
+ for (auto tt : params_derivs_temporary_terms_split[{ 0, 2 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
- for (auto tt : params_derivs_temporary_terms_g12)
+ for (auto tt : params_derivs_temporary_terms_split[{ 1, 2 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
- for (auto tt : params_derivs_temporary_terms_g2)
+ for (auto tt : params_derivs_temporary_terms_split[{ 2, 1 }])
params_derivs_temporary_terms_idxs[tt] = idx++;
}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index ccf27a40..1b21f300 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -34,6 +34,17 @@ using namespace std;
#include "DataTree.hh"
#include "ExtendedPreprocessorTypes.hh"
+// Helper to convert a vector into a tuple
+template <typename T, size_t... Indices>
+auto vectorToTupleHelper(const vector<T>& v, index_sequence<Indices...>) {
+ return make_tuple(v[Indices]...);
+}
+template <size_t N, typename T>
+auto vectorToTuple(const vector<T>& v) {
+ assert(v.size() >= N);
+ return vectorToTupleHelper(v, make_index_sequence<N>());
+}
+
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t >>;
@@ -76,87 +87,44 @@ protected:
//! Stores equation tags
vector<pair<int, pair<string, string>>> equation_tags;
- //! Number of non-zero derivatives
- array<int, 3> NNZDerivatives;
-
- using first_derivatives_t = map<pair<int, int>, expr_t>;
- //! First order derivatives
- /*! First index is equation number, second is variable w.r. to which is computed the derivative.
- Only non-null derivatives are stored in the map.
- Variable indices are those of the getDerivID() method.
- */
- first_derivatives_t first_derivatives;
-
- using second_derivatives_t = map<tuple<int, int, int>, expr_t>;
- //! Second order derivatives
- /*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
- Only non-null derivatives are stored in the map.
- Contains only second order derivatives where var1 >= var2 (for obvious symmetry reasons).
- Variable indices are those of the getDerivID() method.
- */
- second_derivatives_t second_derivatives;
-
- using third_derivatives_t = map<tuple<int, int, int, int>, expr_t>;
- //! Third order derivatives
- /*! First index is equation number, second, third and fourth are variables w.r. to which is computed the derivative.
- Only non-null derivatives are stored in the map.
- Contains only third order derivatives where var1 >= var2 >= var3 (for obvious symmetry reasons).
- Variable indices are those of the getDerivID() method.
- */
- third_derivatives_t third_derivatives;
-
- //! Derivatives of the residuals w.r. to parameters
- /*! First index is equation number, second is parameter.
- Only non-null derivatives are stored in the map.
- Parameter indices are those of the getDerivID() method.
- */
- first_derivatives_t residuals_params_derivatives;
-
- //! Second derivatives of the residuals w.r. to parameters
- /*! First index is equation number, second and third indeces are parameters.
- Only non-null derivatives are stored in the map.
- Parameter indices are those of the getDerivID() method.
- */
- second_derivatives_t residuals_params_second_derivatives;
+ //! Stores derivatives
+ /*! Index 0 is not used, index 1 contains first derivatives, ...
+ For each derivation order, stores a map whose key is a vector of integer: the
+ first integer is the equation index, the remaining ones are the derivation
+ IDs of variables */
+ vector<map<vector<int>, expr_t>> derivatives;
- //! Derivatives of the jacobian w.r. to parameters
- /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter.
- Only non-null derivatives are stored in the map.
- Variable and parameter indices are those of the getDerivID() method.
- */
- second_derivatives_t jacobian_params_derivatives;
-
- //! Second derivatives of the jacobian w.r. to parameters
- /*! First index is equation number, second is endo/exo/exo_det variable, and third and fourth are parameters.
- Only non-null derivatives are stored in the map.
- Variable and parameter indices are those of the getDerivID() method.
- */
- third_derivatives_t jacobian_params_second_derivatives;
-
- //! Derivatives of the hessian w.r. to parameters
- /*! First index is equation number, first and second are endo/exo/exo_det variable, and third is parameter.
- Only non-null derivatives are stored in the map.
- Variable and parameter indices are those of the getDerivID() method.
- */
- third_derivatives_t hessian_params_derivatives;
+ //! Number of non-zero derivatives
+ /*! Index 0 is not used, index 1 contains number of non-zero first
+ derivatives, ... */
+ vector<int> NNZDerivatives;
+
+ //! Derivatives with respect to parameters
+ /*! The key of the outer map is a pair (derivation order w.r.t. endogenous,
+ 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, then the IDs of parameters */
+ map<pair<int,int>, map<vector<int>, expr_t>> params_derivatives;
//! Temporary terms for the static/dynamic file (those which will be noted T[x])
temporary_terms_t temporary_terms;
+
+ //! Temporary terms for model local variables
map<expr_t, expr_t, ExprNodeLess> temporary_terms_mlv;
- temporary_terms_t temporary_terms_res;
- temporary_terms_t temporary_terms_g1;
- temporary_terms_t temporary_terms_g2;
- temporary_terms_t temporary_terms_g3;
+
+ //! Temporary terms for residuals and derivatives
+ /*! Index 0 is temp. terms of residuals, index 1 for first derivatives, ... */
+ vector<temporary_terms_t> temporary_terms_derivatives;
temporary_terms_idxs_t temporary_terms_idxs;
//! Temporary terms for the file containing parameters derivatives
temporary_terms_t params_derivs_temporary_terms;
- temporary_terms_t params_derivs_temporary_terms_res;
- temporary_terms_t params_derivs_temporary_terms_g1;
- temporary_terms_t params_derivs_temporary_terms_res2;
- temporary_terms_t params_derivs_temporary_terms_g12;
- temporary_terms_t params_derivs_temporary_terms_g2;
+
+ //! 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_split;
temporary_terms_idxs_t params_derivs_temporary_terms_idxs;
@@ -384,8 +352,8 @@ public:
/*! Writes either (i+1,j+1) or [i+j*no_eq] */
void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Helper for writing the sparse Hessian or third derivatives in MATLAB and C
- /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[1]]
- If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[2]] */
+ /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[2]]
+ If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[3]] */
void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
inline static std::string
c_Equation_Type(int type)
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 80851ee9..9d6d52e0 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -218,8 +218,8 @@ StaticModel::StaticModel(const DynamicModel &m) :
void
StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
{
- auto it = first_derivatives.find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
- if (it != first_derivatives.end())
+ auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
+ if (it != derivatives[1].end())
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
{
@@ -247,7 +247,6 @@ StaticModel::computeTemporaryTermsOrdered()
map<expr_t, pair<int, int>> first_occurence;
map<expr_t, int> reference_count;
BinaryOpNode *eq_node;
- first_derivatives_t::const_iterator it;
first_chain_rule_derivatives_t::const_iterator it_chr;
ostringstream tmp_s;
v_temporary_terms.clear();
@@ -627,23 +626,23 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
fjmp_if_eval.write(code_file, instruction_number);
int prev_instruction_number = instruction_number;
- vector<vector<pair<int, int>>> derivatives;
- derivatives.resize(symbol_table.endo_nbr());
+ vector<vector<pair<int, int>>> my_derivatives;
+ my_derivatives.resize(symbol_table.endo_nbr());
count_u = symbol_table.endo_nbr();
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
+ int deriv_id = first_derivative.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
expr_t d1 = first_derivative.second;
- unsigned int eq = first_derivative.first.first;
+ unsigned int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
fnumexpr.write(code_file, instruction_number);
- if (!derivatives[eq].size())
- derivatives[eq].clear();
- derivatives[eq].emplace_back(var, count_u);
+ if (!my_derivatives[eq].size())
+ my_derivatives[eq].clear();
+ my_derivatives[eq].emplace_back(var, count_u);
d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
@@ -656,10 +655,10 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
{
FLDR_ fldr(i);
fldr.write(code_file, instruction_number);
- if (derivatives[i].size())
+ if (my_derivatives[i].size())
{
- for (vector<pair<int, int>>::const_iterator it = derivatives[i].begin();
- it != derivatives[i].end(); it++)
+ for (vector<pair<int, int>>::const_iterator it = my_derivatives[i].begin();
+ it != my_derivatives[i].end(); it++)
{
FLDSU_ fldsu(it->second);
fldsu.write(code_file, instruction_number);
@@ -667,7 +666,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
fldsv.write(code_file, instruction_number);
FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
fbinary.write(code_file, instruction_number);
- if (it != derivatives[i].begin())
+ if (it != my_derivatives[i].begin())
{
FBINARY_ fbinary{static_cast<int>(BinaryOpcode::plus)};
fbinary.write(code_file, instruction_number);
@@ -697,20 +696,20 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
tt3.clear();
// The Jacobian if we have to solve the block determinsitic bloc
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
+ int deriv_id = first_derivative.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
expr_t d1 = first_derivative.second;
- unsigned int eq = first_derivative.first.first;
+ unsigned int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
fnumexpr.write(code_file, instruction_number);
- if (!derivatives[eq].size())
- derivatives[eq].clear();
- derivatives[eq].emplace_back(var, count_u);
+ if (!my_derivatives[eq].size())
+ my_derivatives[eq].clear();
+ my_derivatives[eq].emplace_back(var, count_u);
d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
FSTPG2_ fstpg2(eq, var);
@@ -1201,12 +1200,12 @@ map<pair<int, pair<int, int >>, expr_t>
StaticModel::collect_first_order_derivatives_endogenous()
{
map<pair<int, pair<int, int >>, expr_t> endo_derivatives;
- for (auto & first_derivative : first_derivatives)
+ for (auto & first_derivative : derivatives[1])
{
- if (getTypeByDerivID(first_derivative.first.second) == SymbolType::endogenous)
+ if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous)
{
- int eq = first_derivative.first.first;
- int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
+ int eq = first_derivative.first[0];
+ int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
int lag = 0;
endo_derivatives[{ eq, { var, lag } }] = first_derivative.second;
}
@@ -1248,7 +1247,6 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
if (!nopreprocessoroutput)
cout << "Computing static model derivatives:" << endl
<< " - order 1" << endl;
- first_derivatives.clear();
computeJacobian(vars);
@@ -1470,7 +1468,7 @@ StaticModel::writeStaticMatlabCompatLayer(const string &basename) const
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
- int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size();
+ int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
output << "function [residual, g1, g2, g3] = static(y, x, params)" << endl
<< " T = NaN(" << ntt << ", 1);" << endl
@@ -1526,11 +1524,11 @@ StaticModel::writeStaticModel(const string &basename,
writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
model_tt_output, output_type, tef_terms);
- writeTemporaryTerms(temporary_terms_res,
+ writeTemporaryTerms(temporary_terms_derivatives[0],
temp_term_union,
temporary_terms_idxs,
model_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_res.begin(), temporary_terms_res.end());
+ temp_term_union.insert(temporary_terms_derivatives[0].begin(), temporary_terms_derivatives[0].end());
writeModelEquations(model_output, output_type, temp_term_union);
@@ -1539,18 +1537,18 @@ StaticModel::writeStaticModel(const string &basename,
int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
// Write Jacobian w.r. to endogenous only
- if (!first_derivatives.empty())
+ if (!derivatives[1].empty())
{
- writeTemporaryTerms(temporary_terms_g1,
+ writeTemporaryTerms(temporary_terms_derivatives[1],
temp_term_union,
temporary_terms_idxs,
jacobian_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
+ temp_term_union.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
}
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
int eq, var;
- tie(eq, var) = first_derivative.first;
+ tie(eq, var) = vectorToTuple<2>(first_derivative.first);
expr_t d1 = first_derivative.second;
int symb_id = getSymbIDByDerivID(var);
@@ -1563,19 +1561,19 @@ StaticModel::writeStaticModel(const string &basename,
int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
// Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
- if (!second_derivatives.empty())
+ if (!derivatives[2].empty())
{
- writeTemporaryTerms(temporary_terms_g2,
+ writeTemporaryTerms(temporary_terms_derivatives[2],
temp_term_union,
temporary_terms_idxs,
hessian_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
+ temp_term_union.insert(temporary_terms_derivatives[2].begin(), temporary_terms_derivatives[2].end());
int k = 0; // Keep the line of a 2nd derivative in v2
- for (const auto & second_derivative : second_derivatives)
+ for (const auto & second_derivative : derivatives[2])
{
int eq, var1, var2;
- tie(eq, var1, var2) = second_derivative.first;
+ tie(eq, var1, var2) = vectorToTuple<3>(second_derivative.first);
expr_t d2 = second_derivative.second;
int symb_id1 = getSymbIDByDerivID(var1);
@@ -1634,19 +1632,19 @@ StaticModel::writeStaticModel(const string &basename,
}
// Writing third derivatives
- if (!third_derivatives.empty())
+ if (!derivatives[3].empty())
{
- writeTemporaryTerms(temporary_terms_g3,
+ writeTemporaryTerms(temporary_terms_derivatives[3],
temp_term_union,
temporary_terms_idxs,
third_derivatives_tt_output, output_type, tef_terms);
- temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
+ temp_term_union.insert(temporary_terms_derivatives[3].begin(), temporary_terms_derivatives[3].end());
int k = 0; // Keep the line of a 3rd derivative in v3
- for (const auto & third_derivative : third_derivatives)
+ for (const auto & third_derivative : derivatives[3])
{
int eq, var1, var2, var3;
- tie(eq, var1, var2, var3) = third_derivative.first;
+ tie(eq, var1, var2, var3) = vectorToTuple<4>(third_derivative.first);
expr_t d3 = third_derivative.second;
int id1 = getSymbIDByDerivID(var1);
@@ -1731,7 +1729,7 @@ StaticModel::writeStaticModel(const string &basename,
<< " residual = real(residual)+imag(residual).^2;" << endl
<< "end";
writeStaticModelHelper(basename, "static_resid", "residual", "static_resid_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(),
"", init_output, end_output,
model_output, model_tt_output);
@@ -1744,7 +1742,7 @@ StaticModel::writeStaticModel(const string &basename,
<< " g1 = real(g1)+2*imag(g1);" << endl
<< "end";
writeStaticModelHelper(basename, "static_g1", "g1", "static_g1_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
"static_resid_tt",
init_output, end_output,
jacobian_output, jacobian_tt_output);
@@ -1754,16 +1752,16 @@ StaticModel::writeStaticModel(const string &basename,
init_output.clear();
end_output.str(string());
end_output.clear();
- if (second_derivatives.size())
+ if (derivatives[2].size())
{
- init_output << "v2 = zeros(" << NNZDerivatives[1] << ",3);";
+ init_output << "v2 = zeros(" << NNZDerivatives[2] << ",3);";
end_output << "g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");";
}
else
init_output << "g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");";
writeStaticModelHelper(basename, "static_g2", "g2", "static_g2_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size()
- + temporary_terms_g2.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()
+ + temporary_terms_derivatives[2].size(),
"static_g1_tt",
init_output, end_output,
hessian_output, hessian_tt_output);
@@ -1774,16 +1772,16 @@ StaticModel::writeStaticModel(const string &basename,
end_output.str(string());
end_output.clear();
int ncols = hessianColsNbr * JacobianColsNbr;
- if (third_derivatives.size())
+ if (derivatives[3].size())
{
- init_output << "v3 = zeros(" << NNZDerivatives[2] << ",3);";
+ init_output << "v3 = zeros(" << NNZDerivatives[3] << ",3);";
end_output << "g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");";
}
else
init_output << "g3 = sparse([],[],[]," << nrows << "," << ncols << ");";
writeStaticModelHelper(basename, "static_g3", "g3", "static_g3_tt",
- temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size()
- + temporary_terms_g2.size() + temporary_terms_g3.size(),
+ temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()
+ + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size(),
"static_g2_tt",
init_output, end_output,
third_derivatives_output, third_derivatives_tt_output);
@@ -1888,15 +1886,15 @@ StaticModel::writeStaticModel(const string &basename,
// Write the number of temporary terms
output << "tmp_nbr = zeros(Int,4)" << endl
- << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_res.size() << "# Number of temporary terms for the residuals" << endl
- << "tmp_nbr[2] = " << temporary_terms_g1.size() << "# Number of temporary terms for g1 (jacobian)" << endl
- << "tmp_nbr[3] = " << temporary_terms_g2.size() << "# Number of temporary terms for g2 (hessian)" << endl
- << "tmp_nbr[4] = " << temporary_terms_g3.size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
+ << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl
+ << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl
+ << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl
+ << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
// staticResidTT!
output << "function staticResidTT!(T::Vector{Float64}," << endl
<< " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
- << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_res.size() << endl
+ << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl
<< model_tt_output.str()
<< " return nothing" << endl
<< "end" << endl << endl;
@@ -1932,7 +1930,7 @@ StaticModel::writeStaticModel(const string &basename,
output << "function staticG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
<< " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl
<< " @assert length(T) >= "
- << temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() << endl
+ << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
<< " @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << endl
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
@@ -1962,7 +1960,7 @@ StaticModel::writeStaticModel(const string &basename,
output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
<< " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
<< " @assert length(T) >= "
- << temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() << endl
+ << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl
<< " @assert size(g2) == (" << equations.size() << ", " << g2ncols << ")" << endl
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
@@ -1990,7 +1988,7 @@ StaticModel::writeStaticModel(const string &basename,
output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
<< " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
<< " @assert length(T) >= "
- << temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size() << endl
+ << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl
<< " @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
@@ -2043,7 +2041,7 @@ StaticModel::writeStaticCFile(const string &basename) const
string filename = basename + "/model/src/static.c";
string filename_mex = basename + "/model/src/static_mex.c";
- int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size();
+ int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
ofstream output;
output.open(filename, ios::out | ios::binary);
@@ -2146,7 +2144,7 @@ StaticModel::writeStaticCFile(const string &basename) const
<< " if (nlhs >= 3)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v2. */" << endl
- << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3
+ << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3
<< ", mxREAL);" << endl
<< " double *v2 = mxGetPr(plhs[2]);" << endl
<< " static_g2_tt(y, x, nb_row_x, params, T);" << endl
@@ -2250,10 +2248,10 @@ void
StaticModel::writeOutput(ostream &output, bool block) const
{
output << "M_.static_tmp_nbr = zeros(4,1); % Number of temporaries used for the static model" <<endl
- << "M_.static_tmp_nbr(1) = " << temporary_terms_res.size() << "; % Number of temporaries used for the evaluation of the residuals" << endl
- << "M_.static_tmp_nbr(2) = " << temporary_terms_g1.size() << "; % Number of temporaries used for the evaluation of g1 (jacobian)" << endl
- << "M_.static_tmp_nbr(3) = " << temporary_terms_g2.size() << "; % Number of temporaries used for the evaluation of g2 (hessian)" << endl
- << "M_.static_tmp_nbr(4) = " << temporary_terms_g3.size() << "; % Number of temporaries used for the evaluation of g3 (third order derivatives)" << endl;
+ << "M_.static_tmp_nbr(1) = " << temporary_terms_derivatives[0].size() << "; % Number of temporaries used for the evaluation of the residuals" << endl
+ << "M_.static_tmp_nbr(2) = " << temporary_terms_derivatives[1].size() << "; % Number of temporaries used for the evaluation of g1 (jacobian)" << endl
+ << "M_.static_tmp_nbr(3) = " << temporary_terms_derivatives[2].size() << "; % Number of temporaries used for the evaluation of g2 (hessian)" << endl
+ << "M_.static_tmp_nbr(4) = " << temporary_terms_derivatives[3].size() << "; % Number of temporaries used for the evaluation of g3 (third order derivatives)" << endl;
if (!block)
return;
@@ -2290,12 +2288,12 @@ StaticModel::writeOutput(ostream &output, bool block) const
output << "];\n";
map<pair<int, int>, int> row_incidence;
- for (const auto & first_derivative : first_derivatives)
+ for (const auto & first_derivative : derivatives[1])
{
- int deriv_id = first_derivative.first.second;
+ int deriv_id = first_derivative.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
- int eq = first_derivative.first.first;
+ int eq = first_derivative.first[0];
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
//int lag = getLagByDerivID(deriv_id);
@@ -2444,7 +2442,7 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
int eqr = it_l.second.first;
int varr = it_l.second.second;
if (Deriv_type == 0)
- first_chain_rule_derivatives[{ eqr, { varr, lag } }] = first_derivatives[{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
+ first_chain_rule_derivatives[{ eqr, { varr, lag } }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
else if (Deriv_type == 1)
first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
@@ -2505,10 +2503,10 @@ StaticModel::collect_block_first_order_derivatives()
derivative_endo = vector<derivative_t>(nb_blocks);
endo_max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
- for (auto & first_derivative : first_derivatives)
+ for (auto & first_derivative : derivatives[1])
{
- int eq = first_derivative.first.first;
- int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
+ int eq = first_derivative.first[0];
+ int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
int lag = 0;
int block_eq = equation_2_block[eq];
int block_var = variable_2_block[var];
@@ -2518,10 +2516,10 @@ StaticModel::collect_block_first_order_derivatives()
endo_max_leadlag_block[block_eq] = { 0, 0 };
derivative_t tmp_derivative;
lag_var_t lag_var;
- if (getTypeByDerivID(first_derivative.first.second) == SymbolType::endogenous && block_eq == block_var)
+ if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous && block_eq == block_var)
{
tmp_derivative = derivative_endo[block_eq];
- tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
+ tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
derivative_endo[block_eq] = tmp_derivative;
}
}
@@ -2642,11 +2640,7 @@ StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
void
StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) const
{
- if (!residuals_params_derivatives.size()
- && !residuals_params_second_derivatives.size()
- && !jacobian_params_derivatives.size()
- && !jacobian_params_second_derivatives.size()
- && !hessian_params_derivatives.size())
+ if (!params_derivatives.size())
return;
ExprNodeOutputType output_type = (julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel);
@@ -2663,10 +2657,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
writeTemporaryTerms(params_derivs_temporary_terms, {}, params_derivs_temporary_terms_idxs, model_output, output_type, tef_terms);
- for (const auto & residuals_params_derivative : residuals_params_derivatives)
+ for (const auto & residuals_params_derivative : params_derivatives.find({ 0, 1 })->second)
{
int eq, param;
- tie(eq, param) = residuals_params_derivative.first;
+ tie(eq, param) = vectorToTuple<2>(residuals_params_derivative.first);
expr_t d1 = residuals_params_derivative.second;
int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
@@ -2678,10 +2672,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
jacobian_output << ";" << endl;
}
- for (const auto & jacobian_params_derivative : jacobian_params_derivatives)
+ for (const auto & jacobian_params_derivative : params_derivatives.find({ 1, 1 })->second)
{
int eq, var, param;
- tie(eq, var, param) = jacobian_params_derivative.first;
+ tie(eq, var, param) = vectorToTuple<3>(jacobian_params_derivative.first);
expr_t d2 = jacobian_params_derivative.second;
int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
@@ -2695,10 +2689,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
}
int i = 1;
- for (const auto &it : residuals_params_second_derivatives)
+ for (const auto &it : params_derivatives.find({ 0, 2 })->second)
{
int eq, param1, param2;
- tie(eq, param1, param2) = it.first;
+ tie(eq, param1, param2) = vectorToTuple<3>(it.first);
expr_t d2 = it.second;
int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
@@ -2719,10 +2713,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
}
i = 1;
- for (const auto &it : jacobian_params_second_derivatives)
+ for (const auto &it : params_derivatives.find({ 1, 2 })->second)
{
int eq, var, param1, param2;
- tie(eq, var, param1, param2) = it.first;
+ tie(eq, var, param1, param2) = vectorToTuple<4>(it.first);
expr_t d2 = it.second;
int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
@@ -2746,10 +2740,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
}
i = 1;
- for (const auto &it : hessian_params_derivatives)
+ for (const auto &it : params_derivatives.find({ 2, 1 })->second)
{
int eq, var1, var2, param;
- tie(eq, var1, var2, param) = it.first;
+ tie(eq, var1, var2, param) = vectorToTuple<4>(it.first);
expr_t d2 = it.second;
int var1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var1)) + 1;
@@ -2836,13 +2830,13 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
<< symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
<< "if nargout >= 3" << endl
- << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
+ << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
<< hessian1_output.str()
- << "gpp = zeros(" << jacobian_params_second_derivatives.size() << ",5);" << endl
+ << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
<< third_derivs_output.str()
<< "end" << endl
<< "if nargout >= 5" << endl
- << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl
+ << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
<< third_derivs1_output.str()
<< "end" << endl
<< "end" << endl;
@@ -2863,11 +2857,11 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
<< "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
- << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
+ << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
<< hessian1_output.str()
- << "gpp = zeros(" << jacobian_params_second_derivatives.size() << ",5);" << endl
+ << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
<< third_derivs_output.str()
- << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl
+ << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
<< third_derivs1_output.str()
<< "(rp, gp, rpp, gpp, hp)" << endl
<< "end" << endl
@@ -2892,14 +2886,14 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
ostringstream third_derivatives_output; // Used for storing third order derivatives equations
deriv_node_temp_terms_t tef_terms;
- temporary_terms_t temp_term_union = temporary_terms_res;
+ temporary_terms_t temp_term_union = temporary_terms_derivatives[0];
temporary_terms_t temp_term_union_m_1;
string concat = "";
writeJsonModelLocalVariables(model_local_vars_output, tef_terms);
- writeJsonTemporaryTerms(temporary_terms_res, temp_term_union_m_1, model_output, tef_terms, concat);
+ writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union_m_1, model_output, tef_terms, concat);
model_output << ", ";
writeJsonModelEquations(model_output, true);
@@ -2909,21 +2903,21 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
// Write Jacobian w.r. to endogenous only
temp_term_union_m_1 = temp_term_union;
- temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
+ temp_term_union.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
concat = "jacobian";
writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, tef_terms, concat);
jacobian_output << ", \"jacobian\": {"
<< " \"nrows\": " << nrows
<< ", \"ncols\": " << JacobianColsNbr
<< ", \"entries\": [";
- for (auto it = first_derivatives.begin();
- it != first_derivatives.end(); it++)
+ for (auto it = derivatives[1].begin();
+ it != derivatives[1].end(); it++)
{
- if (it != first_derivatives.begin())
+ if (it != derivatives[1].begin())
jacobian_output << ", ";
int eq, var;
- tie(eq, var) = it->first;
+ tie(eq, var) = vectorToTuple<2>(it->first);
int symb_id = getSymbIDByDerivID(var);
int col = symbol_table.getTypeSpecificID(symb_id);
expr_t d1 = it->second;
@@ -2947,21 +2941,21 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
// Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
temp_term_union_m_1 = temp_term_union;
- temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
+ temp_term_union.insert(temporary_terms_derivatives[2].begin(), temporary_terms_derivatives[2].end());
concat = "hessian";
writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, hessian_output, tef_terms, concat);
hessian_output << ", \"hessian\": {"
<< " \"nrows\": " << equations.size()
<< ", \"ncols\": " << g2ncols
<< ", \"entries\": [";
- for (auto it = second_derivatives.begin();
- it != second_derivatives.end(); it++)
+ for (auto it = derivatives[2].begin();
+ it != derivatives[2].end(); it++)
{
- if (it != second_derivatives.begin())
+ if (it != derivatives[2].begin())
hessian_output << ", ";
int eq, var1, var2;
- tie(eq, var1, var2) = it->first;
+ tie(eq, var1, var2) = vectorToTuple<3>(it->first);
int symb_id1 = getSymbIDByDerivID(var1);
int symb_id2 = getSymbIDByDerivID(var2);
expr_t d2 = it->second;
@@ -2994,21 +2988,21 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
// Writing third derivatives
temp_term_union_m_1 = temp_term_union;
- temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
+ temp_term_union.insert(temporary_terms_derivatives[3].begin(), temporary_terms_derivatives[3].end());
concat = "third_derivatives";
writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, third_derivatives_output, tef_terms, concat);
third_derivatives_output << ", \"third_derivative\": {"
<< " \"nrows\": " << equations.size()
<< ", \"ncols\": " << hessianColsNbr * JacobianColsNbr
<< ", \"entries\": [";
- for (auto it = third_derivatives.begin();
- it != third_derivatives.end(); it++)
+ for (auto it = derivatives[3].begin();
+ it != derivatives[3].end(); it++)
{
- if (it != third_derivatives.begin())
+ if (it != derivatives[3].begin())
third_derivatives_output << ", ";
int eq, var1, var2, var3;
- tie(eq, var1, var2, var3) = it->first;
+ tie(eq, var1, var2, var3) = vectorToTuple<4>(it->first);
expr_t d3 = it->second;
if (writeDetails)
@@ -3062,11 +3056,7 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
void
StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const
{
- if (!residuals_params_derivatives.size()
- && !residuals_params_second_derivatives.size()
- && !jacobian_params_derivatives.size()
- && !jacobian_params_second_derivatives.size()
- && !hessian_params_derivatives.size())
+ if (!params_derivatives.size())
return;
ostringstream model_local_vars_output; // Used for storing model local vars
@@ -3087,14 +3077,14 @@ StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
<< " \"neqs\": " << equations.size()
<< ", \"nparamcols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = residuals_params_derivatives.begin();
- it != residuals_params_derivatives.end(); it++)
+ auto &rp = params_derivatives.find({ 0, 1 })->second;
+ for (auto it = rp.begin(); it != rp.end(); it++)
{
- if (it != residuals_params_derivatives.begin())
+ if (it != rp.begin())
jacobian_output << ", ";
int eq, param;
- tie(eq, param) = it->first;
+ tie(eq, param) = vectorToTuple<2>(it->first);
expr_t d1 = it->second;
int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
@@ -3114,19 +3104,20 @@ StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
jacobian_output << "\"}" << endl;
}
jacobian_output << "]}";
+
hessian_output << "\"deriv_jacobian_wrt_params\": {"
<< " \"neqs\": " << equations.size()
<< ", \"nvarcols\": " << symbol_table.endo_nbr()
<< ", \"nparamcols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = jacobian_params_derivatives.begin();
- it != jacobian_params_derivatives.end(); it++)
+ auto &gp = params_derivatives.find({ 1, 1 })->second;
+ for (auto it = gp.begin(); it != gp.end(); it++)
{
- if (it != jacobian_params_derivatives.begin())
+ if (it != gp.begin())
hessian_output << ", ";
int eq, var, param;
- tie(eq, var, param) = it->first;
+ tie(eq, var, param) = vectorToTuple<3>(it->first);
expr_t d2 = it->second;
int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
@@ -3154,14 +3145,14 @@ StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
<< ", \"nparam1cols\": " << symbol_table.param_nbr()
<< ", \"nparam2cols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = residuals_params_second_derivatives.begin();
- it != residuals_params_second_derivatives.end(); ++it)
+ auto &rpp = params_derivatives.find({ 0, 2 })->second;
+ for (auto it = rpp.begin(); it != rpp.end(); ++it)
{
- if (it != residuals_params_second_derivatives.begin())
+ if (it != rpp.begin())
hessian1_output << ", ";
int eq, param1, param2;
- tie(eq, param1, param2) = it->first;
+ tie(eq, param1, param2) = vectorToTuple<3>(it->first);
expr_t d2 = it->second;
int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
@@ -3184,20 +3175,21 @@ StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
hessian1_output << "\"}" << endl;
}
hessian1_output << "]}";
+
third_derivs_output << "\"second_deriv_jacobian_wrt_params\": {"
<< " \"neqs\": " << equations.size()
<< ", \"nvarcols\": " << symbol_table.endo_nbr()
<< ", \"nparam1cols\": " << symbol_table.param_nbr()
<< ", \"nparam2cols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = jacobian_params_second_derivatives.begin();
- it != jacobian_params_second_derivatives.end(); ++it)
+ auto &gpp = params_derivatives.find({ 1, 2 })->second;
+ for (auto it = gpp.begin(); it != gpp.end(); ++it)
{
- if (it != jacobian_params_second_derivatives.begin())
+ if (it != gpp.begin())
third_derivs_output << ", ";
int eq, var, param1, param2;
- tie(eq, var, param1, param2) = it->first;
+ tie(eq, var, param1, param2) = vectorToTuple<4>(it->first);
expr_t d2 = it->second;
int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
@@ -3229,14 +3221,14 @@ StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
<< ", \"nvar2cols\": " << symbol_table.endo_nbr()
<< ", \"nparamcols\": " << symbol_table.param_nbr()
<< ", \"entries\": [";
- for (auto it = hessian_params_derivatives.begin();
- it != hessian_params_derivatives.end(); ++it)
+ auto &hp = params_derivatives.find({ 2, 1 })->second;
+ for (auto it = hp.begin(); it != hp.end(); ++it)
{
- if (it != hessian_params_derivatives.begin())
+ if (it != hp.begin())
third_derivs1_output << ", ";
int eq, var1, var2, param;
- tie(eq, var1, var2, param) = it->first;
+ tie(eq, var1, var2, param) = vectorToTuple<4>(it->first);
expr_t d2 = it->second;
int var1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var1)) + 1;
--
GitLab