diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc index 3368a963..05f02867 100644 --- a/src/ComputingTasks.cc +++ b/src/ComputingTasks.cc @@ -2167,7 +2167,7 @@ PlannerObjectiveStatement::getPlannerObjective() const void PlannerObjectiveStatement::computingPass() { - model_tree.computingPass({}, false, true, true, 0, false, false, false); + model_tree.computingPass(3, 0, {}, false, false, false, false); computing_pass_called = true; } diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 2f76f047..074e9af2 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -4367,11 +4367,11 @@ DynamicModel::substitutePacExpectation() } void -DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, +DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, - bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition) + bool bytecode, bool nopreprocessoroutput, bool linear_decomposition) { - assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder)); + assert(jacobianExo || (derivsOrder < 2 && paramsDerivsOrder == 0)); initializeVariablesAndEquations(); @@ -4393,38 +4393,18 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative // Launch computations if (!nopreprocessoroutput) - { - if (linear_decomposition) - cout << "Computing nonlinear dynamic model derivatives:" << endl - << " - order 1" << endl; - else - cout << "Computing dynamic model derivatives:" << endl - << " - order 1" << endl; - } + cout << "Computing " << (linear_decomposition ? "nonlinear " : "") + << "dynamic model derivatives (order " << derivsOrder << ")." << endl; - computeJacobian(vars); - - if (hessian) - { - if (!nopreprocessoroutput) - cout << " - order 2" << endl; - computeHessian(vars); - } + computeDerivatives(derivsOrder, vars); if (paramsDerivsOrder > 0) { if (!nopreprocessoroutput) - cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl; + cout << "Computing dynamic model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl; computeParamsDerivatives(paramsDerivsOrder); } - if (thirdDerivatives) - { - if (!nopreprocessoroutput) - cout << " - order 3" << endl; - computeThirdDerivatives(vars); - } - jacob_map_t contemporaneous_jacobian, static_jacobian; map >, expr_t> first_order_endo_derivatives; // for each block contains pair diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 341221c7..d84914b5 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -281,14 +281,13 @@ public: //! Execute computations (variable sorting + derivation) /*! \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed) - \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true) - \param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true) + \param derivsOrder order of derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true when order >= 2) \param paramsDerivsOrder order of derivatives w.r. to a pair (endo/exo/exo_det, parameter) to be computed (>0 implies jacobianExo = true) \param eval_context evaluation context for normalization \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */ - void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, - const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition); + void computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder, + const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool nopreprocessoroutput, bool linear_decomposition); //! Writes model initialization and lead/lag incidence matrix to output void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const; diff --git a/src/ModFile.cc b/src/ModFile.cc index ac14a8ee..018eca9d 100644 --- a/src/ModFile.cc +++ b/src/ModFile.cc @@ -701,7 +701,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri { non_linear_equations_dynamic_model = dynamic_model; non_linear_equations_dynamic_model.set_cutoff_to_zero(); - non_linear_equations_dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); + non_linear_equations_dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); } if (!no_static) { @@ -711,13 +711,14 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri || mod_file_struct.calib_smoother_present) static_model.set_cutoff_to_zero(); - const bool static_hessian = mod_file_struct.identification_present - || mod_file_struct.estimation_analytic_derivation; + int derivsOrder = 1; int paramsDerivsOrder = 0; if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation) - paramsDerivsOrder = params_derivs_order; - static_model.computingPass(global_eval_context, no_tmp_terms, static_hessian, - false, paramsDerivsOrder, block, byte_code, nopreprocessoroutput); + { + derivsOrder = 2; + paramsDerivsOrder = params_derivs_order; + } + static_model.computingPass(derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, byte_code, nopreprocessoroutput); } // Set things to compute for dynamic model if (mod_file_struct.perfect_foresight_solver_present || mod_file_struct.check_present @@ -727,7 +728,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri || mod_file_struct.calib_smoother_present) { if (mod_file_struct.perfect_foresight_solver_present) - dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); + dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); else { if (mod_file_struct.stoch_simul_present @@ -740,25 +741,21 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri cerr << "ERROR: Incorrect order option..." << endl; exit(EXIT_FAILURE); } - bool hessian = mod_file_struct.order_option >= 2 - || mod_file_struct.identification_present - || mod_file_struct.estimation_analytic_derivation - || linear - || output == FileOutputType::second - || output == FileOutputType::third; - bool thirdDerivatives = mod_file_struct.order_option == 3 - || mod_file_struct.estimation_analytic_derivation - || output == FileOutputType::third; + int derivsOrder = mod_file_struct.order_option; + if (mod_file_struct.identification_present || linear || output == FileOutputType::second) + derivsOrder = max(derivsOrder, 2); + if (mod_file_struct.estimation_analytic_derivation || output == FileOutputType::third) + derivsOrder = max(derivsOrder, 3); int paramsDerivsOrder = 0; if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation) paramsDerivsOrder = params_derivs_order; - dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); + dynamic_model.computingPass(true, derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); if (linear && mod_file_struct.ramsey_model_present) - orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); + orig_ramsey_dynamic_model.computingPass(true, 2, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); } } else // No computing task requested, compute derivatives up to 2nd order by default - dynamic_model.computingPass(true, true, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); + dynamic_model.computingPass(true, 2, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); map eqs; if (mod_file_struct.ramsey_model_present) @@ -784,7 +781,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri for (auto & statement : statements) statement->computingPass(); - epilogue.computingPass(true, true, false, 0, global_eval_context, true, false, false, false, true, false); + epilogue.computingPass(true, 2, 0, global_eval_context, true, false, false, false, true, false); } void diff --git a/src/ModelTree.cc b/src/ModelTree.cc index bcf4ddf8..db780b53 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -1276,77 +1276,45 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag, } void -ModelTree::computeJacobian(const set &vars) +ModelTree::computeDerivatives(int order, const set &vars) { + assert (order >= 1); + + // Do not shrink the vectors, since they have a minimal size of 4 (see constructor) + derivatives.resize(max(static_cast(order), derivatives.size())); + NNZDerivatives.resize(max(static_cast(order), NNZDerivatives.size()), 0); + + // First-order derivatives for (int var : vars) - { - for (int eq = 0; eq < (int) equations.size(); eq++) + for (int eq = 0; eq < (int) equations.size(); eq++) + { + expr_t d1 = equations[eq]->getDerivative(var); + if (d1 == Zero) + continue; + derivatives[1][{ eq, var }] = d1; + ++NNZDerivatives[1]; + } + + // Higher-order derivatives + for (int o = 2; o <= order; o++) + for (const auto &it : derivatives[o-1]) + for (int var : vars) { - expr_t d1 = equations[eq]->getDerivative(var); - if (d1 == Zero) + if (it.first.back() > var) continue; - derivatives[1][{ eq, var }] = d1; - ++NNZDerivatives[1]; + + expr_t d = it.second->getDerivative(var); + if (d == Zero) + continue; + + vector indices{it.first}; + indices.push_back(var); + // At this point, indices of endogenous variables are sorted in non-decreasing order + derivatives[o][indices] = d; + do + NNZDerivatives[o]++; + while (next_permutation(next(indices.begin()), indices.end())); } - } -} - -void -ModelTree::computeHessian(const set &vars) -{ - for (const auto &it : derivatives[1]) - { - int eq, var1; - tie(eq, var1) = vectorToTuple<2>(it.first); - expr_t d1 = it.second; - - // Store only second derivatives with var2 <= var1 - for (int var2 : vars) - { - if (var2 > var1) - continue; - - expr_t d2 = d1->getDerivative(var2); - if (d2 == Zero) - continue; - derivatives[2][{ eq, var1, var2 }] = d2; - if (var2 == var1) - ++NNZDerivatives[2]; - else - NNZDerivatives[2] += 2; - } - } -} - -void -ModelTree::computeThirdDerivatives(const set &vars) -{ - for (const auto &it : derivatives[2]) - { - int eq, var1, var2; - tie(eq, var1, var2) = vectorToTuple<3>(it.first); - // By construction, var2 <= var1 - - expr_t d2 = it.second; - - // Store only third derivatives such that var3 <= var2 <= var1 - for (int var3 : vars) - { - if (var3 > var2) - continue; - - expr_t d3 = d2->getDerivative(var3); - if (d3 == Zero) - continue; - derivatives[3][{ eq, var1, var2, var3 }] = d3; - if (var3 == var2 && var2 == var1) - ++NNZDerivatives[3]; - else if (var3 == var2 || var2 == var1) - NNZDerivatives[3] += 3; - else - NNZDerivatives[3] += 6; - } - } } void diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 2b2c076a..8612ded8 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -91,7 +91,8 @@ protected: /*! 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 */ + IDs of variables (in non-decreasing order, to avoid storing symmetric + elements several times) */ vector, expr_t>> derivatives; //! Number of non-zero derivatives @@ -104,7 +105,7 @@ protected: derivation order w.r.t. parameters). For e.g., { 1, 2 } corresponds to the jacobian differentiated twice w.r.t. to parameters. In inner maps, the vector of integers consists of: the equation index, then - the derivation IDs of endogenous, then the IDs of parameters */ + the derivation IDs of endogenous (in non-decreasing order), then the IDs of parameters */ map, map, expr_t>> params_derivatives; //! Storage for temporary terms in block/bytecode mode @@ -144,15 +145,10 @@ protected: //! Vector indicating if the equation is linear in endogenous variable (true) or not (false) vector is_equation_linear; - //! Computes 1st derivatives - /*! \param vars the derivation IDs w.r. to which compute the derivatives */ - void computeJacobian(const set &vars); - //! Computes 2nd derivatives - /*! \param vars the derivation IDs w.r. to which derive the 1st derivatives */ - void computeHessian(const set &vars); - //! Computes 3rd derivatives - /*! \param vars the derivation IDs w.r. to which derive the 2nd derivatives */ - void computeThirdDerivatives(const set &vars); + //! Computes derivatives + /*! \param order the derivation order + \param vars the derivation IDs w.r.t. which compute the derivatives */ + void computeDerivatives(int order, const set &vars); //! Computes derivatives of the Jacobian and Hessian w.r. to parameters void computeParamsDerivatives(int paramsDerivsOrder); //! Write derivative of an equation w.r. to a variable diff --git a/src/StaticModel.cc b/src/StaticModel.cc index e61a5b68..64229a76 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -1214,7 +1214,7 @@ StaticModel::collect_first_order_derivatives_endogenous() } void -StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, bool block, bool bytecode, const bool nopreprocessoroutput) +StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool bytecode, bool nopreprocessoroutput) { initializeVariablesAndEquations(); @@ -1233,9 +1233,9 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms equations.clear(); copy(neweqs.begin(), neweqs.end(), back_inserter(equations)); - // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic - set vars; + // Compute derivatives w.r. to all endogenous + set vars; for (int i = 0; i < symbol_table.endo_nbr(); i++) { int id = symbol_table.getID(SymbolType::endogenous, i); @@ -1245,29 +1245,14 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms // Launch computations if (!nopreprocessoroutput) - cout << "Computing static model derivatives:" << endl - << " - order 1" << endl; + cout << "Computing static model derivatives (order " << derivsOrder << ")." << endl; - computeJacobian(vars); - - if (hessian) - { - if (!nopreprocessoroutput) - cout << " - order 2" << endl; - computeHessian(vars); - } - - if (thirdDerivatives) - { - if (!nopreprocessoroutput) - cout << " - order 3" << endl; - computeThirdDerivatives(vars); - } + computeDerivatives(derivsOrder, vars); if (paramsDerivsOrder > 0) { if (!nopreprocessoroutput) - cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl; + cout << "Computing static model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl; computeParamsDerivatives(paramsDerivsOrder); } diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 28d4ad87..de40a312 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -193,10 +193,10 @@ public: /*! \param eval_context evaluation context for normalization \param no_tmp_terms if true, no temporary terms will be computed in the static files - \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed - \param paramsDerivsOrder order of derivatives w.r. to a pair (endo/exo/exo_det, parameter) to be computed + \param derivsOrder order of derivation with respect to endogenous + \param paramsDerivsOrder order of derivatives w.r. to a pair (endogenous, parameter) to be computed */ - void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatices, int paramsDerivsOrder, bool block, bool bytecode, const bool nopreprocessoroutput); + void computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool bytecode, bool nopreprocessoroutput); //! Adds informations for simulation in a binary file for a block decomposed model void Write_Inf_To_Bin_File_Block(const string &basename, const int &num,