From 6f2af6943fae65b3f236d78de375b2f1be995503 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 7 Feb 2024 18:50:34 +0100 Subject: [PATCH] k_order_welfare MEX: use the sparse representation of the static model --- .../k_order_welfare/k_ord_objective.cc | 69 +------ .../k_order_welfare/k_ord_objective.hh | 5 +- .../k_order_welfare/k_order_welfare.cc | 51 +++-- mex/sources/k_order_welfare/objective_m.cc | 191 ++++++++++-------- mex/sources/k_order_welfare/objective_m.hh | 25 ++- 5 files changed, 159 insertions(+), 182 deletions(-) diff --git a/mex/sources/k_order_welfare/k_ord_objective.cc b/mex/sources/k_order_welfare/k_ord_objective.cc index 943ce8dd5..480999345 100644 --- a/mex/sources/k_order_welfare/k_ord_objective.cc +++ b/mex/sources/k_order_welfare/k_ord_objective.cc @@ -22,21 +22,19 @@ #include #include -KordwDynare::KordwDynare(KordpDynare& m, ConstVector& NNZD_arg, Journal& jr, Vector& inParams, +KordwDynare::KordwDynare(KordpDynare& m, Journal& jr, Vector& inParams, std::unique_ptr objectiveFile_arg, const std::vector& dr_order) : model {m}, - NNZD {NNZD_arg}, journal {jr}, params {inParams}, resid(1), ud {1}, objectiveFile {std::move(objectiveFile_arg)} { - dynppToDyn = dr_order; dynToDynpp.resize(model.ny()); for (int i = 0; i < model.ny(); i++) - dynToDynpp[dynppToDyn[i]] = i; + dynToDynpp[dr_order[i]] = i; } void @@ -45,71 +43,10 @@ KordwDynare::calcDerivativesAtSteady() assert(ud.begin() == ud.end()); - std::vector dyn_ud; // Planner's objective derivatives, in Dynare form - dyn_ud.emplace_back(1, model.ny()); // Allocate Jacobian - dyn_ud.back().zeros(); - - for (int i = 2; i <= model.order(); i++) - { - // Higher order derivatives, as sparse (3-column) matrices - dyn_ud.emplace_back(static_cast(NNZD[i - 1]), 3); - dyn_ud.back().zeros(); - } - Vector xx(model.nexog()); xx.zeros(); resid.zeros(); - objectiveFile->eval(model.getSteady(), xx, params, resid, dyn_ud); - - for (int i = 1; i <= model.order(); i++) - populateDerivativesContainer(dyn_ud, i); -} - -void -KordwDynare::populateDerivativesContainer(const std::vector& dyn_ud, int ord) -{ - const TwoDMatrix& u = dyn_ud[ord - 1]; - - // utility derivatives FSSparseTensor instance - auto udTi = std::make_unique(ord, model.ny(), 1); - - IntSequence s(ord, 0); - - if (ord == 1) - for (int i = 0; i < u.ncols(); i++) - { - for (int j = 0; j < u.nrows(); j++) - { - double x = u.get(j, dynppToDyn[s[0]]); - if (x != 0.0) - udTi->insert(s, j, x); - } - s[0]++; - } - else // ord ≥ 2 - for (int i = 0; i < u.nrows(); i++) - { - int j = static_cast(u.get(i, 0)) - 1; - int i1 = static_cast(u.get(i, 1)) - 1; - if (j < 0 || i1 < 0) - continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix()) - - for (int k = 0; k < ord; k++) - { - s[k] = dynToDynpp[i1 % model.ny()]; - i1 /= model.ny(); - } - - if (ord == 2 && !s.isSorted()) - continue; // Skip symmetric elements (only needed at order 2) - else if (ord > 2) - s.sort(); // For higher order, canonicalize the multi-index - - double x = u.get(i, 2); - udTi->insert(s, j, x); - } - - ud.insert(std::move(udTi)); + objectiveFile->eval(model.getSteady(), xx, params, resid, dynToDynpp, ud); } template<> diff --git a/mex/sources/k_order_welfare/k_ord_objective.hh b/mex/sources/k_order_welfare/k_ord_objective.hh index 0acc4932e..329ffe97b 100644 --- a/mex/sources/k_order_welfare/k_ord_objective.hh +++ b/mex/sources/k_order_welfare/k_ord_objective.hh @@ -29,22 +29,19 @@ class KordwDynare { public: KordpDynare& model; - const ConstVector& NNZD; private: Journal& journal; Vector& params; Vector resid; TensorContainer ud; // planner's objective derivatives, in Dynare++ form - std::vector dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones std::vector dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones std::unique_ptr objectiveFile; public: - KordwDynare(KordpDynare& m, ConstVector& NNZD_arg, Journal& jr, Vector& inParams, + KordwDynare(KordpDynare& m, Journal& jr, Vector& inParams, std::unique_ptr objectiveFile_arg, const std::vector& varOrder); void calcDerivativesAtSteady(); - void populateDerivativesContainer(const std::vector& dyn_ud, int ord); [[nodiscard]] const TensorContainer& getPlannerObjDerivatives() const { diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc index 04c2cbe7c..130afa58a 100644 --- a/mex/sources/k_order_welfare/k_order_welfare.cc +++ b/mex/sources/k_order_welfare/k_order_welfare.cc @@ -31,6 +31,7 @@ #include #include +#include #include "dynmex.h" @@ -214,15 +215,6 @@ extern "C" mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you " "used the right order option inside the `stoch_simul' command"); - const mxArray* nnzderivatives_obj_mx = mxGetField(M_mx, 0, "NNZDerivatives_objective"); - if (!(nnzderivatives_obj_mx && mxIsDouble(nnzderivatives_obj_mx) - && !mxIsComplex(nnzderivatives_obj_mx) && !mxIsSparse(nnzderivatives_obj_mx))) - mexErrMsgTxt("M_.NNZDerivatives should be a real dense array"); - ConstVector NNZD_obj {nnzderivatives_obj_mx}; - if (NNZD.length() < kOrder || NNZD_obj[kOrder - 1] == -1) - mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you " - "used the right order option inside the `stoch_simul' command"); - const mxArray* endo_names_mx = mxGetField(M_mx, 0, "endo_names"); if (!(endo_names_mx && mxIsCell(endo_names_mx) && mxGetNumberOfElements(endo_names_mx) == static_cast(nEndo))) @@ -261,6 +253,32 @@ extern "C" std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx) + nEndo, dr_order.begin(), [](double x) { return static_cast(x) - 1; }); + const mxArray* objective_g1_sparse_rowval_mx + = mxGetField(M_mx, 0, "objective_g1_sparse_rowval"); + if (!(objective_g1_sparse_rowval_mx && mxIsInt32(objective_g1_sparse_rowval_mx))) + mexErrMsgTxt("M_.objective_g1_sparse_rowval should be an int32 array"); + + const mxArray* objective_g1_sparse_colval_mx + = mxGetField(M_mx, 0, "objective_g1_sparse_colval"); + if (!(objective_g1_sparse_colval_mx && mxIsInt32(objective_g1_sparse_colval_mx))) + mexErrMsgTxt("M_.objective_g1_sparse_colval should be an int32 array"); + + const mxArray* objective_g1_sparse_colptr_mx + = mxGetField(M_mx, 0, "objective_g1_sparse_colptr"); + if (!(objective_g1_sparse_colptr_mx && mxIsInt32(objective_g1_sparse_colptr_mx))) + mexErrMsgTxt("M_.objective_g1_sparse_colptr should be an int32 array"); + + std::vector objective_gN_sparse_indices; + for (int o {2}; o <= kOrder; o++) + { + using namespace std::string_literals; + auto fieldname {"objective_g"s + std::to_string(o) + "_sparse_indices"}; + const mxArray* indices = mxGetField(M_mx, 0, fieldname.c_str()); + if (!(indices && mxIsInt32(indices))) + mexErrMsgTxt(("M_."s + fieldname + " should be an int32 array").c_str()); + objective_gN_sparse_indices.push_back(indices); + } + const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state @@ -291,21 +309,14 @@ extern "C" // run stochastic steady app.walkStochSteady(); - const mxArray* objective_tmp_nbr_mx = mxGetField(M_mx, 0, "objective_tmp_nbr"); - if (!(objective_tmp_nbr_mx && mxIsDouble(objective_tmp_nbr_mx) - && !mxIsComplex(objective_tmp_nbr_mx) && !mxIsSparse(objective_tmp_nbr_mx) - && mxGetNumberOfElements(objective_tmp_nbr_mx) >= static_cast(kOrder + 1))) - mexErrMsgTxt("M_.objective_tmp_nbr should be a real dense array with strictly more elements " - "than the order of derivation"); - int ntt_objective = std::accumulate(mxGetPr(objective_tmp_nbr_mx), - mxGetPr(objective_tmp_nbr_mx) + kOrder + 1, 0); - // Getting derivatives of the planner's objective function std::unique_ptr objectiveFile; - objectiveFile = std::make_unique(fname, ntt_objective); + objectiveFile = std::make_unique( + fname, kOrder, objective_g1_sparse_rowval_mx, objective_g1_sparse_colval_mx, + objective_g1_sparse_colptr_mx, objective_gN_sparse_indices); // make KordwDynare object - KordwDynare welfare(dynare, NNZD_obj, journal, modParams, std::move(objectiveFile), dr_order); + KordwDynare welfare(dynare, journal, modParams, std::move(objectiveFile), dr_order); // construct main K-order approximation class of welfare ApproximationWelfare appwel(welfare, discount_factor, app.get_rule_ders(), diff --git a/mex/sources/k_order_welfare/objective_m.cc b/mex/sources/k_order_welfare/objective_m.cc index 1996cbb39..ff523d90b 100644 --- a/mex/sources/k_order_welfare/objective_m.cc +++ b/mex/sources/k_order_welfare/objective_m.cc @@ -26,88 +26,41 @@ #include "objective_m.hh" -ObjectiveMFile::ObjectiveMFile(const std::string& modName, int ntt_arg) : - ntt {ntt_arg}, ObjectiveMFilename {modName + ".objective.static"} +ObjectiveMFile::ObjectiveMFile(const std::string& modName, int kOrder_arg, + const mxArray* objective_g1_sparse_rowval_mx_arg, + const mxArray* objective_g1_sparse_colval_mx_arg, + const mxArray* objective_g1_sparse_colptr_mx_arg, + const std::vector objective_gN_sparse_indices_arg) : + ObjectiveMFilename {modName + ".objective.sparse.static"}, + kOrder {kOrder_arg}, + objective_g1_sparse_rowval_mx {objective_g1_sparse_rowval_mx_arg}, + objective_g1_sparse_colval_mx {objective_g1_sparse_colval_mx_arg}, + objective_g1_sparse_colptr_mx {objective_g1_sparse_colptr_mx_arg}, + objective_gN_sparse_indices {objective_gN_sparse_indices_arg} { } -void -ObjectiveMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray* sparseMat, TwoDMatrix& tdm) -{ - int totalCols = mxGetN(sparseMat); - mwIndex* rowIdxVector = mxGetIr(sparseMat); - mwIndex* colIdxVector = mxGetJc(sparseMat); - - assert(tdm.ncols() == 3); - /* Under MATLAB, the following check always holds at equality; under Octave, - there may be an inequality, because Octave diminishes nzmax if one gives - zeros in the values vector when calling sparse(). */ - assert(tdm.nrows() >= static_cast(mxGetNzmax(sparseMat))); - - double* ptr = mxGetPr(sparseMat); - - int rind = 0; - int output_row = 0; - - for (int i = 0; i < totalCols; i++) - for (int j = 0; j < static_cast((colIdxVector[i + 1] - colIdxVector[i])); j++, rind++) - { - tdm.get(output_row, 0) = rowIdxVector[rind] + 1; - tdm.get(output_row, 1) = i + 1; - tdm.get(output_row, 2) = ptr[rind]; - output_row++; - } - - /* If there are less elements than expected (that might happen if some - derivative is symbolically not zero but numerically zero at the evaluation - point), then fill in the matrix with empty entries, that will be - recognized as such by KordpDynare::populateDerivativesContainer() */ - while (output_row < tdm.nrows()) - { - tdm.get(output_row, 0) = 0; - tdm.get(output_row, 1) = 0; - tdm.get(output_row, 2) = 0; - output_row++; - } -} - void ObjectiveMFile::eval(const Vector& y, const Vector& x, const Vector& modParams, Vector& residual, - std::vector& md) noexcept(false) + const std::vector& dynToDynpp, + TensorContainer& derivatives) const { - mxArray* T_m = mxCreateDoubleMatrix(ntt, 1, mxREAL); + mxArray* y_mx = mxCreateDoubleMatrix(y.length(), 1, mxREAL); + std::copy_n(y.base(), y.length(), mxGetPr(y_mx)); - mxArray* y_m = mxCreateDoubleMatrix(y.length(), 1, mxREAL); - std::copy_n(y.base(), y.length(), mxGetPr(y_m)); + mxArray* x_mx = mxCreateDoubleMatrix(1, x.length(), mxREAL); + std::copy_n(x.base(), x.length(), mxGetPr(x_mx)); - mxArray* x_m = mxCreateDoubleMatrix(1, x.length(), mxREAL); - std::copy_n(x.base(), x.length(), mxGetPr(x_m)); + mxArray* params_mx = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL); + std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_mx)); - mxArray* params_m = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL); - std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_m)); - - mxArray* T_flag_m = mxCreateLogicalScalar(false); - - { - // Compute temporary terms (for all orders) - std::string funcname = ObjectiveMFilename + "_g" + std::to_string(md.size()) + "_tt"; - std::array plhs; - std::array prhs {T_m, y_m, x_m, params_m}; - - int retVal - = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str()); - if (retVal != 0) - throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname); - - mxDestroyArray(T_m); - T_m = plhs[0]; - } + mxArray *T_order_mx, *T_mx; { // Compute residuals std::string funcname = ObjectiveMFilename + "_resid"; - std::array plhs; - std::array prhs {T_m, y_m, x_m, params_m, T_flag_m}; + std::array plhs; + std::array prhs {y_mx, x_mx, params_mx}; int retVal = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str()); @@ -116,35 +69,103 @@ ObjectiveMFile::eval(const Vector& y, const Vector& x, const Vector& modParams, residual = Vector {plhs[0]}; mxDestroyArray(plhs[0]); + + T_order_mx = plhs[1]; + T_mx = plhs[2]; } - for (size_t i = 1; i <= md.size(); i++) + { + // Compute Jacobian + std::string funcname = ObjectiveMFilename + "_g1"; + std::array plhs; + std::array prhs {y_mx, + x_mx, + params_mx, + const_cast(objective_g1_sparse_rowval_mx), + const_cast(objective_g1_sparse_colval_mx), + const_cast(objective_g1_sparse_colptr_mx), + T_order_mx, + T_mx}; + + int retVal + = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str()); + if (retVal != 0) + throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname); + + assert(static_cast(mxGetN(plhs[0])) == y.length()); + double* g1_v {mxGetPr(plhs[0])}; + mwIndex* g1_ir {mxGetIr(plhs[0])}; + mwIndex* g1_jc {mxGetJc(plhs[0])}; + + IntSequence s(1, 0); + auto tensor = std::make_unique(1, y.length(), 1); + + for (int j {0}; j < y.length(); j++) + for (mwIndex k {g1_jc[j]}; k < g1_jc[j + 1]; k++) + { + s[0] = dynToDynpp[j]; + tensor->insert(s, g1_ir[k], g1_v[k]); + } + + mxDestroyArray(plhs[0]); + + mxDestroyArray(T_order_mx); + T_order_mx = plhs[1]; + + mxDestroyArray(T_mx); + T_mx = plhs[2]; + + derivatives.insert(std::move(tensor)); + } + + for (int o {2}; o <= kOrder; o++) { - // Compute model derivatives - std::string funcname = ObjectiveMFilename + "_g" + std::to_string(i); - std::array plhs; - std::array prhs {T_m, y_m, x_m, params_m, T_flag_m}; + // Compute higher derivatives + std::string funcname = ObjectiveMFilename + "_g" + std::to_string(o); + std::array plhs; + std::array prhs {y_mx, x_mx, params_mx, T_order_mx, T_mx}; int retVal = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str()); if (retVal != 0) throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname); - if (i == 1) + const mxArray* sparse_indices_mx {objective_gN_sparse_indices[o - 2]}; + size_t nnz {mxGetM(sparse_indices_mx)}; +#if MX_HAS_INTERLEAVED_COMPLEX + const int32_T* sparse_indices {mxGetInt32s(sparse_indices_mx)}; +#else + const int32_T* sparse_indices {static_cast(mxGetData(sparse_indices_mx))}; +#endif + + assert(mxGetNumberOfElements(plhs[0]) == nnz); + double* gN_v {mxGetPr(plhs[0])}; + + IntSequence s(o, 0); + auto tensor = std::make_unique(o, y.length(), 1); + + for (size_t k {0}; k < nnz; k++) { - assert(static_cast(mxGetM(plhs[0])) == md[i - 1].nrows()); - assert(static_cast(mxGetN(plhs[0])) == md[i - 1].ncols()); - std::copy_n(mxGetPr(plhs[0]), mxGetM(plhs[0]) * mxGetN(plhs[0]), md[i - 1].base()); + for (int i {0}; i < o; i++) + s[i] = dynToDynpp[sparse_indices[k + (i + 1) * nnz] - 1]; + assert(s.isSorted()); + tensor->insert(s, sparse_indices[k] - 1, gN_v[k]); } - else - unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[0], md[i - 1]); mxDestroyArray(plhs[0]); + + mxDestroyArray(T_order_mx); + T_order_mx = plhs[1]; + + mxDestroyArray(T_mx); + T_mx = plhs[2]; + + derivatives.insert(std::move(tensor)); } - mxDestroyArray(T_m); - mxDestroyArray(y_m); - mxDestroyArray(x_m); - mxDestroyArray(params_m); - mxDestroyArray(T_flag_m); + mxDestroyArray(y_mx); + mxDestroyArray(x_mx); + mxDestroyArray(params_mx); + mxDestroyArray(T_order_mx); + mxDestroyArray(T_mx); } diff --git a/mex/sources/k_order_welfare/objective_m.hh b/mex/sources/k_order_welfare/objective_m.hh index f6c5f1ee8..d720f910b 100644 --- a/mex/sources/k_order_welfare/objective_m.hh +++ b/mex/sources/k_order_welfare/objective_m.hh @@ -22,21 +22,32 @@ #include -#include +#include "dynmex.h" -#include "twod_matrix.hh" +#include "Vector.hh" +#include "sparse_tensor.hh" +#include "t_container.hh" -// Handles calls to /+objective/static.m +// Handles calls to /+objective/+sparse/static*.m class ObjectiveMFile { private: - int ntt; // Size of vector of temporary terms const std::string ObjectiveMFilename; - static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray* sparseMat, TwoDMatrix& tdm); + const int kOrder; + const mxArray *const objective_g1_sparse_rowval_mx, *const objective_g1_sparse_colval_mx, + *const objective_g1_sparse_colptr_mx; + // Stores M_.objective_gN_sparse_indices, starting from N=2 + const std::vector objective_gN_sparse_indices; public: - explicit ObjectiveMFile(const std::string& modName, int ntt_arg); + ObjectiveMFile(const std::string& modName, int kOrder_arg, + const mxArray* objective_g1_sparse_rowval_mx_arg, + const mxArray* objective_g1_sparse_colval_mx_arg, + const mxArray* objective_g1_sparse_colptr_mx_arg, + const std::vector objective_gN_sparse_indices_arg); void eval(const Vector& y, const Vector& x, const Vector& params, Vector& residual, - std::vector& md); + const std::vector& dynToDynpp, TensorContainer& derivatives) const + noexcept(false); }; + #endif