diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 912f261e..2a1094dd 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -245,8 +245,8 @@ DynamicModel::operator=(const DynamicModel &m) void DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const { - auto it = first_derivatives.find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) }); - if (it != first_derivatives.end()) + auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) }); + if (it != derivatives[1].end()) (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false); else { @@ -274,7 +274,6 @@ DynamicModel::computeTemporaryTermsOrdered() map> first_occurence; map 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>, expr_t> first_derivatives_reordered_endo; map, pair>, 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, int >>> derivatives; - derivatives.resize(symbol_table.endo_nbr()); + vector, 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, int>>::const_iterator it = derivatives[i].begin(); - it != derivatives[i].end(); it++) + for (vector, 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(BinaryOpcode::times)}; fbinary.write(code_file, instruction_number); - if (it != derivatives[i].begin()) + if (it != my_derivatives[i].begin()) { FBINARY_ fbinary{static_cast(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 &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>, 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>, expr_t> DynamicModel::collect_first_order_derivatives_endogenous() { map>, 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>(nb_blocks, { 0, 0 }); exo_det_max_leadlag_block = vector>(nb_blocks, { 0, 0 }); max_leadlag_block = vector>(nb_blocks, { 0, 0 }); - for (auto & first_derivative : first_derivatives) + for (auto & first_derivative : derivatives[1]) { - int eq = first_derivative.first.first; - int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second)); - int lag = getLagByDerivID(first_derivative.first.second); + int eq = first_derivative.first[0]; + int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1])); + int lag = getLagByDerivID(first_derivative.first[1]); int block_eq = equation_2_block[eq]; int block_var = 0; derivative_t tmp_derivative; lag_var_t lag_var; - switch (getTypeByDerivID(first_derivative.first.second)) + switch (getTypeByDerivID(first_derivative.first[1])) { case SymbolType::endogenous: block_var = variable_2_block[var]; @@ -4783,7 +4782,7 @@ DynamicModel::collect_block_first_order_derivatives() if (lag > 0 && lag > endo_max_leadlag_block[block_eq].second) endo_max_leadlag_block[block_eq] = { endo_max_leadlag_block[block_eq].first, lag }; tmp_derivative = derivative_endo[block_eq]; - tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }]; + tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }]; derivative_endo[block_eq] = tmp_derivative; } else @@ -4807,7 +4806,7 @@ DynamicModel::collect_block_first_order_derivatives() } } } - tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }]; + tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }]; derivative_other_endo[block_eq] = tmp_derivative; lag_var = other_endo_block[block_eq]; if (lag_var.find(lag) == lag_var.end()) @@ -4836,7 +4835,7 @@ DynamicModel::collect_block_first_order_derivatives() } } } - tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }]; + tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }]; derivative_exo[block_eq] = tmp_derivative; lag_var = exo_block[block_eq]; if (lag_var.find(lag) == lag_var.end()) @@ -4864,7 +4863,7 @@ DynamicModel::collect_block_first_order_derivatives() } } } - tmp_derivative[{ lag, { eq, var } }] = first_derivatives[{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }]; + tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }]; derivative_exo_det[block_eq] = tmp_derivative; lag_var = exo_det_block[block_eq]; if (lag_var.find(lag) == lag_var.end()) @@ -5378,11 +5377,7 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context void DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) const { - if (!residuals_params_derivatives.size() - && !residuals_params_second_derivatives.size() - && !jacobian_params_derivatives.size() - && !jacobian_params_second_derivatives.size() - && !hessian_params_derivatives.size()) + if (!params_derivatives.size()) return; ExprNodeOutputType output_type = (julia ? ExprNodeOutputType::juliaDynamicModel : ExprNodeOutputType::matlabDynamicModel); @@ -5398,10 +5393,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con writeTemporaryTerms(params_derivs_temporary_terms, {}, params_derivs_temporary_terms_idxs, model_output, output_type, tef_terms); - for (const auto & residuals_params_derivative : residuals_params_derivatives) + for (const auto & residuals_params_derivative : params_derivatives.find({ 0, 1 })->second) { int eq, param; - tie(eq, param) = residuals_params_derivative.first; + tie(eq, param) = vectorToTuple<2>(residuals_params_derivative.first); expr_t d1 = residuals_params_derivative.second; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; @@ -5412,10 +5407,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con jacobian_output << ";" << endl; } - for (const auto & jacobian_params_derivative : jacobian_params_derivatives) + for (const auto & jacobian_params_derivative : params_derivatives.find({ 1, 1 })->second) { int eq, var, param; - tie(eq, var, param) = jacobian_params_derivative.first; + tie(eq, var, param) = vectorToTuple<3>(jacobian_params_derivative.first); expr_t d2 = jacobian_params_derivative.second; int var_col = getDynJacobianCol(var) + 1; @@ -5428,10 +5423,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con } int i = 1; - for (const auto &it : residuals_params_second_derivatives) + for (const auto &it : params_derivatives.find({ 0, 2 })->second) { int eq, param1, param2; - tie(eq, param1, param2) = it.first; + tie(eq, param1, param2) = vectorToTuple<3>(it.first); expr_t d2 = it.second; int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1; @@ -5452,10 +5447,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con } i = 1; - for (const auto &it : jacobian_params_second_derivatives) + for (const auto &it : params_derivatives.find({ 1, 2 })->second) { int eq, var, param1, param2; - tie(eq, var, param1, param2) = it.first; + tie(eq, var, param1, param2) = vectorToTuple<4>(it.first); expr_t d2 = it.second; int var_col = getDynJacobianCol(var) + 1; @@ -5479,10 +5474,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con } i = 1; - for (const auto &it : hessian_params_derivatives) + for (const auto &it : params_derivatives.find({ 2, 1 })->second) { int eq, var1, var2, param; - tie(eq, var1, var2, param) = it.first; + tie(eq, var1, var2, param) = vectorToTuple<4>(it.first); expr_t d2 = it.second; int var1_col = getDynJacobianCol(var1) + 1; @@ -5579,13 +5574,13 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "if nargout >= 3" << endl - << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl + << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl << hessian1_output.str() - << "gpp = zeros(" << jacobian_params_second_derivatives.size() << ",5);" << endl + << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl << third_derivs_output.str() << "end" << endl << "if nargout >= 5" << endl - << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl + << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl << third_derivs1_output.str() << "end" << endl << "end" << endl; @@ -5606,11 +5601,11 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << jacobian_output.str() << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() - << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl + << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl << hessian1_output.str() - << "gpp = zeros(" << jacobian_params_second_derivatives.size() << ",5);" << endl + << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl << third_derivs_output.str() - << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl + << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl << third_derivs1_output.str() << "(rp, gp, rpp, gpp, hp)" << endl << "end" << endl @@ -6371,7 +6366,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c deriv_node_temp_terms_t tef_terms; temporary_terms_t temp_term_empty; - temporary_terms_t temp_term_union = temporary_terms_res; + temporary_terms_t temp_term_union = temporary_terms_derivatives[0]; temporary_terms_t temp_term_union_m_1; string concat = ""; @@ -6379,27 +6374,27 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c writeJsonModelLocalVariables(model_local_vars_output, tef_terms); - writeJsonTemporaryTerms(temporary_terms_res, temp_term_union_m_1, model_output, tef_terms, concat); + writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union_m_1, model_output, tef_terms, concat); model_output << ", "; writeJsonModelEquations(model_output, true); // Writing Jacobian temp_term_union_m_1 = temp_term_union; - temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end()); + temp_term_union.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end()); concat = "jacobian"; writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, tef_terms, concat); jacobian_output << ", \"jacobian\": {" << " \"nrows\": " << equations.size() << ", \"ncols\": " << dynJacobianColsNbr << ", \"entries\": ["; - for (auto it = first_derivatives.begin(); - it != first_derivatives.end(); it++) + for (auto it = derivatives[1].begin(); + it != derivatives[1].end(); it++) { - if (it != first_derivatives.begin()) + if (it != derivatives[1].begin()) jacobian_output << ", "; int eq, var; - tie(eq, var) = it->first; + tie(eq, var) = vectorToTuple<2>(it->first); int col = getDynJacobianCol(var); expr_t d1 = it->second; @@ -6422,21 +6417,21 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c // Writing Hessian temp_term_union_m_1 = temp_term_union; - temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end()); + temp_term_union.insert(temporary_terms_derivatives[2].begin(), temporary_terms_derivatives[2].end()); concat = "hessian"; writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, hessian_output, tef_terms, concat); hessian_output << ", \"hessian\": {" << " \"nrows\": " << equations.size() << ", \"ncols\": " << hessianColsNbr << ", \"entries\": ["; - for (auto it = second_derivatives.begin(); - it != second_derivatives.end(); it++) + for (auto it = derivatives[2].begin(); + it != derivatives[2].end(); it++) { - if (it != second_derivatives.begin()) + if (it != derivatives[2].begin()) hessian_output << ", "; int eq, var1, var2; - tie(eq, var1, var2) = it->first; + tie(eq, var1, var2) = vectorToTuple<3>(it->first); expr_t d2 = it->second; int id1 = getDynJacobianCol(var1); int id2 = getDynJacobianCol(var2); @@ -6467,21 +6462,21 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c // Writing third derivatives temp_term_union_m_1 = temp_term_union; - temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end()); + temp_term_union.insert(temporary_terms_derivatives[3].begin(), temporary_terms_derivatives[3].end()); concat = "third_derivatives"; writeJsonTemporaryTerms(temp_term_union, temp_term_union_m_1, third_derivatives_output, tef_terms, concat); third_derivatives_output << ", \"third_derivative\": {" << " \"nrows\": " << equations.size() << ", \"ncols\": " << hessianColsNbr * dynJacobianColsNbr << ", \"entries\": ["; - for (auto it = third_derivatives.begin(); - it != third_derivatives.end(); it++) + for (auto it = derivatives[3].begin(); + it != derivatives[3].end(); it++) { - if (it != third_derivatives.begin()) + if (it != derivatives[3].begin()) third_derivatives_output << ", "; int eq, var1, var2, var3; - tie(eq, var1, var2, var3) = it->first; + tie(eq, var1, var2, var3) = vectorToTuple<4>(it->first); expr_t d3 = it->second; if (writeDetails) @@ -6538,11 +6533,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c void DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const { - if (!residuals_params_derivatives.size() - && !residuals_params_second_derivatives.size() - && !jacobian_params_derivatives.size() - && !jacobian_params_second_derivatives.size() - && !hessian_params_derivatives.size()) + if (!params_derivatives.size()) return; ostringstream model_local_vars_output; // Used for storing model local vars @@ -6563,14 +6554,14 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) << " \"neqs\": " << equations.size() << ", \"nparamcols\": " << symbol_table.param_nbr() << ", \"entries\": ["; - for (auto it = residuals_params_derivatives.begin(); - it != residuals_params_derivatives.end(); it++) + auto &rp = params_derivatives.find({ 0, 1 })->second; + for (auto it = rp.begin(); it != rp.end(); it++) { - if (it != residuals_params_derivatives.begin()) + if (it != rp.begin()) jacobian_output << ", "; int eq, param; - tie(eq, param) = it->first; + tie(eq, param) = vectorToTuple<2>(it->first); expr_t d1 = it->second; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; @@ -6590,19 +6581,20 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) jacobian_output << "\"}" << endl; } jacobian_output << "]}"; + hessian_output << "\"deriv_jacobian_wrt_params\": {" << " \"neqs\": " << equations.size() << ", \"nvarcols\": " << dynJacobianColsNbr << ", \"nparamcols\": " << symbol_table.param_nbr() << ", \"entries\": ["; - for (auto it = jacobian_params_derivatives.begin(); - it != jacobian_params_derivatives.end(); it++) + auto &gp = params_derivatives.find({ 1, 1 })->second; + for (auto it = gp.begin(); it != gp.end(); it++) { - if (it != jacobian_params_derivatives.begin()) + if (it != gp.begin()) hessian_output << ", "; int eq, var, param; - tie(eq, var, param) = it->first; + tie(eq, var, param) = vectorToTuple<3>(it->first); expr_t d2 = it->second; int var_col = getDynJacobianCol(var) + 1; @@ -6632,14 +6624,14 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) << ", \"nparam1cols\": " << symbol_table.param_nbr() << ", \"nparam2cols\": " << symbol_table.param_nbr() << ", \"entries\": ["; - for (auto it = residuals_params_second_derivatives.begin(); - it != residuals_params_second_derivatives.end(); ++it) + auto &rpp = params_derivatives.find({ 0, 2 })->second; + for (auto it = rpp.begin(); it != rpp.end(); ++it) { - if (it != residuals_params_second_derivatives.begin()) + if (it != rpp.begin()) hessian1_output << ", "; int eq, param1, param2; - tie(eq, param1, param2) = it->first; + tie(eq, param1, param2) = vectorToTuple<3>(it->first); expr_t d2 = it->second; int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1; @@ -6661,20 +6653,21 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) hessian1_output << "\"}" << endl; } hessian1_output << "]}"; + third_derivs_output << "\"second_deriv_jacobian_wrt_params\": {" << " \"neqs\": " << equations.size() << ", \"nvarcols\": " << dynJacobianColsNbr << ", \"nparam1cols\": " << symbol_table.param_nbr() << ", \"nparam2cols\": " << symbol_table.param_nbr() << ", \"entries\": ["; - for (auto it = jacobian_params_second_derivatives.begin(); - it != jacobian_params_second_derivatives.end(); ++it) + auto &gpp = params_derivatives.find({ 1, 2 })->second; + for (auto it = gpp.begin(); it != gpp.end(); ++it) { - if (it != jacobian_params_second_derivatives.begin()) + if (it != gpp.begin()) third_derivs_output << ", "; int eq, var, param1, param2; - tie(eq, var, param1, param2) = it->first; + tie(eq, var, param1, param2) = vectorToTuple<4>(it->first); expr_t d2 = it->second; int var_col = getDynJacobianCol(var) + 1; @@ -6708,14 +6701,14 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) << ", \"nvar2cols\": " << dynJacobianColsNbr << ", \"nparamcols\": " << symbol_table.param_nbr() << ", \"entries\": ["; - for (auto it = hessian_params_derivatives.begin(); - it != hessian_params_derivatives.end(); ++it) + auto &hp = params_derivatives.find({ 2, 1 })->second; + for (auto it = hp.begin(); it != hp.end(); ++it) { - if (it != hessian_params_derivatives.begin()) + if (it != hp.begin()) third_derivs1_output << ", "; int eq, var1, var2, param; - tie(eq, var1, var2, param) = it->first; + tie(eq, var1, var2, param) = vectorToTuple<4>(it->first); expr_t d2 = it->second; int var1_col = getDynJacobianCol(var1) + 1; diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 4e3d2959..76d878eb 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -43,51 +43,41 @@ ModelTree::copyHelper(const ModelTree &m) for (const auto & it : m.aux_equations) aux_equations.push_back(dynamic_cast(f(it))); + auto convert_deriv_map = [f](map, expr_t> dm) + { + map, 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> jacobian_elements_to_delete; - for (first_derivatives_t::const_iterator it = first_derivatives.begin(); - it != first_derivatives.end(); it++) + set> 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 &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 &vars) void ModelTree::computeHessian(const set &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 &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 &vars) void ModelTree::computeThirdDerivatives(const set &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 &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> 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 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::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> reference_count; params_derivs_temporary_terms.clear(); map 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::const_iterator it = temp_terms_map.begin(); it != temp_terms_map.end(); it++) params_derivs_temporary_terms.insert(it->second.begin(), it->second.end()); - params_derivs_temporary_terms_res = temp_terms_map[NodeTreeReference::residualsParamsDeriv]; - params_derivs_temporary_terms_g1 = temp_terms_map[NodeTreeReference::jacobianParamsDeriv]; - params_derivs_temporary_terms_res2 = temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv]; - params_derivs_temporary_terms_g12 = temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv]; - params_derivs_temporary_terms_g2 = temp_terms_map[NodeTreeReference::hessianParamsDeriv]; + params_derivs_temporary_terms_split[{ 0, 1 }] = temp_terms_map[NodeTreeReference::residualsParamsDeriv]; + params_derivs_temporary_terms_split[{ 1, 1 }] = temp_terms_map[NodeTreeReference::jacobianParamsDeriv]; + params_derivs_temporary_terms_split[{ 0, 2 }] = temp_terms_map[NodeTreeReference::residualsParamsSecondDeriv]; + params_derivs_temporary_terms_split[{ 1, 2 }] = temp_terms_map[NodeTreeReference::jacobianParamsSecondDeriv]; + params_derivs_temporary_terms_split[{ 2, 1 }] = temp_terms_map[NodeTreeReference::hessianParamsDeriv]; int idx = 0; - for (auto tt : params_derivs_temporary_terms_res) + for (auto tt : params_derivs_temporary_terms_split[{ 0, 1 }]) params_derivs_temporary_terms_idxs[tt] = idx++; - for (auto tt : params_derivs_temporary_terms_g1) + for (auto tt : params_derivs_temporary_terms_split[{ 1, 1 }]) params_derivs_temporary_terms_idxs[tt] = idx++; - for (auto tt : params_derivs_temporary_terms_res2) + for (auto tt : params_derivs_temporary_terms_split[{ 0, 2 }]) params_derivs_temporary_terms_idxs[tt] = idx++; - for (auto tt : params_derivs_temporary_terms_g12) + for (auto tt : params_derivs_temporary_terms_split[{ 1, 2 }]) params_derivs_temporary_terms_idxs[tt] = idx++; - for (auto tt : params_derivs_temporary_terms_g2) + for (auto tt : params_derivs_temporary_terms_split[{ 2, 1 }]) params_derivs_temporary_terms_idxs[tt] = idx++; } diff --git a/src/ModelTree.hh b/src/ModelTree.hh index ccf27a40..1b21f300 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -34,6 +34,17 @@ using namespace std; #include "DataTree.hh" #include "ExtendedPreprocessorTypes.hh" +// Helper to convert a vector into a tuple +template +auto vectorToTupleHelper(const vector& v, index_sequence) { + return make_tuple(v[Indices]...); +} +template +auto vectorToTuple(const vector& v) { + assert(v.size() >= N); + return vectorToTupleHelper(v, make_index_sequence()); +} + //! 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>; @@ -76,87 +87,44 @@ protected: //! Stores equation tags vector>> 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, expr_t>> derivatives; + //! Number of non-zero derivatives - array NNZDerivatives; + /*! Index 0 is not used, index 1 contains number of non-zero first + derivatives, ... */ + vector NNZDerivatives; - using first_derivatives_t = map, 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, 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, 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, map, 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 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_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, temporary_terms_t> params_derivs_temporary_terms_split; temporary_terms_idxs_t params_derivs_temporary_terms_idxs; @@ -384,8 +352,8 @@ public: /*! Writes either (i+1,j+1) or [i+j*no_eq] */ void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const; //! Helper for writing the sparse Hessian or third derivatives in MATLAB and C - /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[1]] - If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[2]] */ + /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[2]] + If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[3]] */ void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const; inline static std::string c_Equation_Type(int type) diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 80851ee9..9d6d52e0 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -218,8 +218,8 @@ StaticModel::StaticModel(const DynamicModel &m) : void StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const { - auto it = first_derivatives.find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) }); - if (it != first_derivatives.end()) + auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) }); + if (it != derivatives[1].end()) (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false); else { @@ -247,7 +247,6 @@ StaticModel::computeTemporaryTermsOrdered() map> first_occurence; map 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>> derivatives; - derivatives.resize(symbol_table.endo_nbr()); + vector>> 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>::const_iterator it = derivatives[i].begin(); - it != derivatives[i].end(); it++) + for (vector>::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(BinaryOpcode::times)}; fbinary.write(code_file, instruction_number); - if (it != derivatives[i].begin()) + if (it != my_derivatives[i].begin()) { FBINARY_ fbinary{static_cast(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>, expr_t> StaticModel::collect_first_order_derivatives_endogenous() { map>, 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" <, 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(nb_blocks); endo_max_leadlag_block = vector>(nb_blocks, { 0, 0 }); max_leadlag_block = vector>(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;