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

Make storage of derivatives ready for k-order 

This is the first step towards computing k-order derivatives.

Ref Dynare/dynare#217, #10
issue#70
Sébastien Villemot 2018-11-15 16:39:53 +01:00
parent 3eff2b91c2
commit e89dd9cb8b
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
4 changed files with 432 additions and 507 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

287
src/DynamicModel.cc
View File

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

268
src/ModelTree.cc
View File

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

112
src/ModelTree.hh
View File

@ -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;
//! 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;
//! Number of non-zero derivatives
array<int, 3> NNZDerivatives;
/*! Index 0 is not used, index 1 contains number of non-zero first
derivatives, ... */
vector<int> 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;
//! 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;
//! 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)

272
src/StaticModel.cc
View File

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