diff --git a/mex/build/k_order_perturbation.am b/mex/build/k_order_perturbation.am index 0b3e8be84..e58d0ce43 100644 --- a/mex/build/k_order_perturbation.am +++ b/mex/build/k_order_perturbation.am @@ -15,7 +15,6 @@ nodist_k_order_perturbation_SOURCES = \ $(TOPDIR)/k_ord_dynare.hh \ $(TOPDIR)/dynamic_dll.cc \ $(TOPDIR)/dynamic_dll.hh \ - $(TOPDIR)/dynamic_abstract_class.cc \ $(TOPDIR)/dynamic_abstract_class.hh \ $(TOPDIR)/dynamic_m.cc \ $(TOPDIR)/dynamic_m.hh diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc deleted file mode 100644 index 76bcb6461..000000000 --- a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc +++ /dev/null @@ -1,72 +0,0 @@ -/* - * Copyright (C) 2010-2019 Dynare Team - * - * This file is part of Dynare. - * - * Dynare is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * Dynare is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with Dynare. If not, see . - */ - -#include - -#include "dynamic_abstract_class.hh" - -void -DynamicModelAC::copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols) -{ - assert(rows == tdm->nrows()); - assert(cols == tdm->ncols()); - - int dmIdx = 0; - for (int j = 0; j < cols; j++) - for (int i = 0; i < rows; i++) - tdm->get(i, j) = dm[dmIdx++]; -} - -void -DynamicModelAC::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm) -{ - int totalCols = mxGetN(sparseMat); - mwIndex *rowIdxVector = mxGetIr(sparseMat); - mwSize sizeRowIdxVector = mxGetNzmax(sparseMat); - mwIndex *colIdxVector = mxGetJc(sparseMat); - - assert(tdm->ncols() == 3); - assert(tdm->nrows() == sizeRowIdxVector); - - 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 Nzmax (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 < static_cast(sizeRowIdxVector)) - { - tdm->get(output_row, 0) = 0; - tdm->get(output_row, 1) = 0; - tdm->get(output_row, 2) = 0; - output_row++; - } -} diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh index 9cdd2621c..3a79df7db 100644 --- a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh +++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh @@ -20,15 +20,13 @@ #ifndef _DYNAMICMODELAC_HH #define _DYNAMICMODELAC_HH -#include "k_ord_dynare.hh" +#include "twod_matrix.hh" class DynamicModelAC { public: virtual ~DynamicModelAC() = default; - static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm); - static void copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols); virtual void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, - Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false) = 0; + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) = 0; }; #endif diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cc b/mex/sources/k_order_perturbation/dynamic_dll.cc index 719a93a46..ba8acd85c 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.cc +++ b/mex/sources/k_order_perturbation/dynamic_dll.cc @@ -19,6 +19,10 @@ #include "dynamic_dll.hh" +#include "dynare_exception.hh" + +#include + DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false) { std::string fName; @@ -27,71 +31,51 @@ DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false) #endif fName += "+" + modName + "/dynamic" + MEXEXT; - try - { #if defined(__CYGWIN32__) || defined(_WIN32) - dynamicHinstance = LoadLibrary(fName.c_str()); - if (!dynamicHinstance) - throw 1; - ntt = reinterpret_cast(GetProcAddress(dynamicHinstance, "ntt")); - dynamic_resid_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_resid_tt")); - dynamic_resid = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_resid")); - dynamic_g1_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g1_tt")); - dynamic_g1 = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g1")); - dynamic_g2_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g2_tt")); - dynamic_g2 = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g2")); - dynamic_g3_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g3_tt")); - dynamic_g3 = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g3")); - if (!ntt || !dynamic_resid_tt || !dynamic_resid - || !dynamic_g1_tt || !dynamic_g1 - || !dynamic_g2_tt || !dynamic_g2 - || !dynamic_g3_tt || !dynamic_g3) - { - FreeLibrary(dynamicHinstance); // Free the library - throw 2; - } + dynamicHinstance = LoadLibrary(fName.c_str()); + if (!dynamicHinstance) + throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": can't dynamically load the file"); + ntt = reinterpret_cast(GetProcAddress(dynamicHinstance, "ntt")); + dynamic_resid_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_resid_tt")); + dynamic_resid = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_resid")); + dynamic_g1_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g1_tt")); + dynamic_g1 = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g1")); + dynamic_g2_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g2_tt")); + dynamic_g2 = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g2")); + dynamic_g3_tt = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g3_tt")); + dynamic_g3 = reinterpret_cast(GetProcAddress(dynamicHinstance, "dynamic_g3")); + if (!ntt || !dynamic_resid_tt || !dynamic_resid + || !dynamic_g1_tt || !dynamic_g1 + || !dynamic_g2_tt || !dynamic_g2 + || !dynamic_g3_tt || !dynamic_g3) + { + FreeLibrary(dynamicHinstance); // Free the library + throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": can't locate the relevant dynamic symbols within the MEX file"); + } #else // Linux or Mac - dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW); - if (dynamicHinstance == nullptr) - { - std::cerr << dlerror() << std::endl; - throw 1; - } - ntt = reinterpret_cast(dlsym(dynamicHinstance, "ntt")); - dynamic_resid_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_resid_tt")); - dynamic_resid = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_resid")); - dynamic_g1_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g1_tt")); - dynamic_g1 = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g1")); - dynamic_g2_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g2_tt")); - dynamic_g2 = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g2")); - dynamic_g3_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g3_tt")); - dynamic_g3 = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g3")); - if (!ntt || !dynamic_resid_tt || !dynamic_resid - || !dynamic_g1_tt || !dynamic_g1 - || !dynamic_g2_tt || !dynamic_g2 - || !dynamic_g3_tt || !dynamic_g3) - { - dlclose(dynamicHinstance); // Free the library - std::cerr << dlerror() << std::endl; - throw 2; - } + dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW); + if (!dynamicHinstance) + throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": " + dlerror()); + ntt = reinterpret_cast(dlsym(dynamicHinstance, "ntt")); + dynamic_resid_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_resid_tt")); + dynamic_resid = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_resid")); + dynamic_g1_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g1_tt")); + dynamic_g1 = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g1")); + dynamic_g2_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g2_tt")); + dynamic_g2 = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g2")); + dynamic_g3_tt = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g3_tt")); + dynamic_g3 = reinterpret_cast(dlsym(dynamicHinstance, "dynamic_g3")); + if (!ntt || !dynamic_resid_tt || !dynamic_resid + || !dynamic_g1_tt || !dynamic_g1 + || !dynamic_g2_tt || !dynamic_g2 + || !dynamic_g3_tt || !dynamic_g3) + { + dlclose(dynamicHinstance); // Free the library + throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": " + dlerror()); + } #endif - } - catch (int i) - { - std::string msg{"Error when loading " + fName + " ("}; - if (i == 1) - msg += "can't dynamically load the file"; - if (i == 2) - msg += "can't locate the relevant dynamic symbols within the MEX file"; - msg += ")"; - throw DynareException(__FILE__, __LINE__, msg); - } - catch (...) - { - throw DynareException(__FILE__, __LINE__, "Can't find the relevant dynamic symbols in " + fName); - } + tt = std::make_unique(*ntt); } DynamicModelDLL::~DynamicModelDLL() @@ -112,21 +96,19 @@ void DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false) { - double *T = new double[*ntt]; - dynamic_resid_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T); - dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, residual.base()); + dynamic_resid_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get()); + dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), residual.base()); if (g1 || g2 || g3) - dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T); + dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get()); if (g1) - dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g1->base()); + dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g1->base()); if (g2 || g3) - dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T); + dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get()); if (g2) - dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g2->base()); + dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g2->base()); if (g3) { - dynamic_g3_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T); - dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g3->base()); + dynamic_g3_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get()); + dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g3->base()); } - delete[] T; } diff --git a/mex/sources/k_order_perturbation/dynamic_dll.hh b/mex/sources/k_order_perturbation/dynamic_dll.hh index 3e0568029..bfb1456dc 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.hh +++ b/mex/sources/k_order_perturbation/dynamic_dll.hh @@ -30,9 +30,9 @@ #endif #include +#include #include "dynamic_abstract_class.hh" -#include "dynare_exception.hh" using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T); using dynamic_resid_fct = void (*) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *residual); @@ -58,13 +58,14 @@ private: #else void *dynamicHinstance; // and in Linux or Mac #endif + std::unique_ptr tt; // Vector of temporary terms public: // construct and load Dynamic model DLL - explicit DynamicModelDLL(const std::string &fname) noexcept(false); + explicit DynamicModelDLL(const std::string &fname); virtual ~DynamicModelDLL(); void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, - Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false); + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override; }; #endif diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc index 1ef2dfa8d..0b05d53cf 100644 --- a/mex/sources/k_order_perturbation/dynamic_m.cc +++ b/mex/sources/k_order_perturbation/dynamic_m.cc @@ -18,6 +18,9 @@ */ #include +#include + +#include "dynare_exception.hh" #include "dynamic_m.hh" @@ -26,10 +29,49 @@ DynamicModelMFile::DynamicModelMFile(const std::string &modName) noexcept(false) { } +void +DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm) +{ + int totalCols = mxGetN(sparseMat); + mwIndex *rowIdxVector = mxGetIr(sparseMat); + mwSize sizeRowIdxVector = mxGetNzmax(sparseMat); + mwIndex *colIdxVector = mxGetJc(sparseMat); + + assert(tdm->ncols() == 3); + assert(tdm->nrows() == sizeRowIdxVector); + + 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 Nzmax (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 < static_cast(sizeRowIdxVector)) + { + tdm->get(output_row, 0) = 0; + tdm->get(output_row, 1) = 0; + tdm->get(output_row, 2) = 0; + output_row++; + } +} + void DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false) { + constexpr int nlhs_dynamic = 4, nrhs_dynamic = 5; mxArray *prhs[nrhs_dynamic], *plhs[nlhs_dynamic]; prhs[0] = mxCreateDoubleMatrix(y.length(), 1, mxREAL); @@ -48,8 +90,11 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam throw DynareException(__FILE__, __LINE__, "Trouble calling " + DynamicMFilename); residual = Vector{plhs[0]}; - copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, static_cast(mxGetM(plhs[1])), - static_cast(mxGetN(plhs[1]))); + + assert(static_cast(mxGetM(plhs[1])) == g1->nrows()); + assert(static_cast(mxGetN(plhs[1])) == g1->ncols()); + std::copy_n(mxGetPr(plhs[1]), mxGetM(plhs[1])*mxGetN(plhs[1]), g1->base()); + if (g2) unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2); if (g3) diff --git a/mex/sources/k_order_perturbation/dynamic_m.hh b/mex/sources/k_order_perturbation/dynamic_m.hh index 4632f3a33..95d6dc317 100644 --- a/mex/sources/k_order_perturbation/dynamic_m.hh +++ b/mex/sources/k_order_perturbation/dynamic_m.hh @@ -33,12 +33,11 @@ class DynamicModelMFile : public DynamicModelAC { private: const std::string DynamicMFilename; - constexpr static int nlhs_dynamic = 4; - constexpr static int nrhs_dynamic = 5; + static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm); public: - explicit DynamicModelMFile(const std::string &modName) noexcept(false); + explicit DynamicModelMFile(const std::string &modName); virtual ~DynamicModelMFile() = default; void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, - Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false); + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override; }; #endif diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc index be1dd34ff..744bc0519 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cc +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc @@ -19,74 +19,43 @@ // GP, based on work by O.Kamenik -#include -#include "first_order.hh" +#include "k_ord_dynare.hh" #include "dynamic_abstract_class.hh" +#include "dynare_exception.hh" -#include -#include +#include -#include -#include -/**************************************************************************************/ -/* Dynare DynamicModel class */ -/**************************************************************************************/ - -KordpDynare::KordpDynare(const std::vector &endo, int num_endo, +KordpDynare::KordpDynare(const std::vector &endo, const std::vector &exo, int nexog, int npar, Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat, - int npred, int nforw, int nboth, const int jcols, const Vector &nnzd, - const int nsteps, int norder, - Journal &jr, std::unique_ptr dynamicModelFile_arg, double sstol, - const std::vector &var_order, const TwoDMatrix &llincidence, double criterium) noexcept(false) : - nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar}, - nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps}, - nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov}, - md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order}, - ll_Incidence{llincidence}, qz_criterium{criterium}, - dynamicModelFile{move(dynamicModelFile_arg)} -{ - ReorderDynareJacobianIndices(); - - // Initialise ModelDerivativeContainer(*this, this->md, nOrder); - for (int iord = 1; iord <= nOrder; iord++) - md.insert(std::make_unique(iord, nY+nYs+nYss+nExog, nY)); -} - -KordpDynare::KordpDynare(const std::vector &endo, int num_endo, - const std::vector &exo, int nexog, int npar, - Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat, - int npred, int nforw, int nboth, const int jcols, const Vector &nnzd, - const int nsteps, int norder, - Journal &jr, std::unique_ptr dynamicModelFile_arg, double sstol, - const std::vector &var_order, const TwoDMatrix &llincidence, double criterium, + int npred, int nforw, int nboth, const Vector &nnzd, + int nsteps, int norder, + Journal &jr, std::unique_ptr dynamicModelFile_arg, + const std::vector &var_order, const TwoDMatrix &llincidence, std::unique_ptr g1_arg, std::unique_ptr g2_arg, - std::unique_ptr g3_arg) noexcept(false) : + std::unique_ptr g3_arg) : nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar}, - nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps}, + nYs{npred + nboth}, nYss{nboth + nforw}, nY{nstat + npred + nboth + nforw}, + nJcols{nExog+nY+nYs+nYss}, NNZD{nnzd}, nSteps{nsteps}, nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov}, - md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order}, - ll_Incidence{llincidence}, qz_criterium{criterium}, + md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, + ll_Incidence{llincidence}, g1p{std::move(g1_arg)}, g2p{std::move(g2_arg)}, g3p{std::move(g3_arg)}, - dynamicModelFile{move(dynamicModelFile_arg)} + dynamicModelFile{std::move(dynamicModelFile_arg)} { - ReorderDynareJacobianIndices(); - - // Initialise ModelDerivativeContainer(*this, this->md, nOrder); - for (int iord = 1; iord <= nOrder; iord++) - md.insert(std::make_unique(iord, nY+nYs+nYss+nExog, nY)); + computeJacobianPermutation(var_order); } void KordpDynare::solveDeterministicSteady() { JournalRecordPair pa(journal); - pa << "Non-linear solver for deterministic steady state By-passed " << endrec; + pa << "Non-linear solver for deterministic steady state skipped" << endrec; } void -KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false) +KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) { // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals out.zeros(); @@ -94,16 +63,12 @@ KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx void KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, - const ConstVector &yyp, const Vector &xx) noexcept(false) + const ConstVector &yyp, const Vector &xx) { // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals out.zeros(); } -/************************************************ - * this is main derivative calculation functin that indirectly calls dynamic.dll - * which performs actual calculation and reorders - ***************************************************/ void KordpDynare::calcDerivativesAtSteady() { @@ -136,20 +101,17 @@ KordpDynare::calcDerivativesAtSteady() dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p.get(), g2p.get(), g3p.get()); } - populateDerivativesContainer(*g1p, 1, JacobianIndices); + populateDerivativesContainer(*g1p, 1); if (nOrder > 1) - populateDerivativesContainer(*g2p, 2, JacobianIndices); + populateDerivativesContainer(*g2p, 2); if (nOrder > 2) - populateDerivativesContainer(*g3p, 3, JacobianIndices); + populateDerivativesContainer(*g3p, 3); } -/******************************************************************************* - * populateDerivatives to sparse Tensor and fit it in the Derivatives Container - *******************************************************************************/ void -KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const std::vector &vOrder) +KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord) { // model derivatives FSSparseTensor instance auto mdTi = std::make_unique(ord, nJcols, nY); @@ -161,11 +123,7 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st { for (int j = 0; j < g.nrows(); j++) { - double x; - if (s[0] < nJcols-nExog) - x = g.get(j, vOrder[s[0]]); - else - x = g.get(j, s[0]); + double x = g.get(j, dynppToDyn[s[0]]); if (x != 0.0) mdTi->insert(s, j, x); } @@ -173,26 +131,16 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st } else if (ord == 2) { - int nJcols1 = nJcols-nExog; - std::vector revOrder(nJcols1); - for (int i = 0; i < nJcols1; i++) - revOrder[vOrder[i]] = i; for (int i = 0; i < g.nrows(); i++) { int j = static_cast(g.get(i, 0))-1; // hessian indices start with 1 - int i1 = static_cast(g.get(i, 1)) -1; + int i1 = static_cast(g.get(i, 1))-1; if (j < 0 || i1 < 0) continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix()) int s0 = i1 / nJcols; int s1 = i1 % nJcols; - if (s0 < nJcols1) - s[0] = revOrder[s0]; - else - s[0] = s0; - if (s1 < nJcols1) - s[1] = revOrder[s1]; - else - s[1] = s1; + s[0] = dynToDynpp[s0]; + s[1] = dynToDynpp[s1]; if (s[1] >= s[0]) { double x = g.get(i, 2); @@ -202,67 +150,49 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st } else if (ord == 3) { - int nJcols1 = nJcols-nExog; int nJcols2 = nJcols*nJcols; - std::vector revOrder(nJcols1); - for (int i = 0; i < nJcols1; i++) - revOrder[vOrder[i]] = i; for (int i = 0; i < g.nrows(); i++) { int j = static_cast(g.get(i, 0))-1; - int i1 = static_cast(g.get(i, 1)) -1; + int i1 = static_cast(g.get(i, 1))-1; if (j < 0 || i1 < 0) continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix()) int s0 = i1 / nJcols2; int i2 = i1 % nJcols2; int s1 = i2 / nJcols; int s2 = i2 % nJcols; - if (s0 < nJcols1) - s[0] = revOrder[s0]; - else - s[0] = s0; - if (s1 < nJcols1) - s[1] = revOrder[s1]; - else - s[1] = s1; - if (s2 < nJcols1) - s[2] = revOrder[s2]; - else - s[2] = s2; - double x = g.get(i, 2); + s[0] = dynToDynpp[s0]; + s[1] = dynToDynpp[s1]; + s[2] = dynToDynpp[s2]; if (s.isSorted()) - mdTi->insert(s, j, x); + { + double x = g.get(i, 2); + mdTi->insert(s, j, x); + } } } - // md container - md.remove(Symmetry{ord}); md.insert(std::move(mdTi)); } -/********************************************************* - * LLxSteady() - * returns ySteady extended with leads and lags suitable for - * passing to _dynamic DLL - *************************************************************/ +/* Returns ySteady extended with leads and lags suitable for passing to + _dynamic */ void -KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false) +KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) { - if ((nJcols-nExog) == yS.length()) + if (yS.length() == nJcols-nExog) throw DynareException(__FILE__, __LINE__, "ySteady already of right size"); - // create temporary square 2D matrix size nEndo x nEndo (sparse) - // for the lag, current and lead blocks of the jacobian + /* Create temporary square 2D matrix size nEndo×nEndo (sparse) + for the lag, current and lead blocks of the jacobian */ if (llxSteady.length() != nJcols-nExog) throw DynareException(__FILE__, __LINE__, "llxSteady has wrong size"); for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++) // populate (non-sparse) vector with ysteady values for (int i = 0; i < nY; i++) - { - if (ll_Incidence.get(ll_row, i)) - llxSteady[static_cast(ll_Incidence.get(ll_row, i))-1] = yS[i]; - } + if (ll_Incidence.get(ll_row, i)) + llxSteady[static_cast(ll_Incidence.get(ll_row, i))-1] = yS[i]; } /************************************ @@ -291,26 +221,26 @@ KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false) ************************************/ void -KordpDynare::ReorderDynareJacobianIndices() noexcept(false) +KordpDynare::computeJacobianPermutation(const std::vector &var_order) { - // create temporary square 2D matrix size nEndo x nEndo (sparse) + dynppToDyn.resize(nJcols); + // create temporary square 2D matrix size nEndo×nEndo (sparse) // for the lag, current and lead blocks of the jacobian - JacobianIndices.resize(nJcols); std::vector tmp(nY); - int i, j, rjoff = nJcols-nExog-1; + int rjoff = nJcols-nExog-1; for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++) { // reorder in orde-var order & populate temporary nEndo (sparse) vector with // the lag, current and lead blocks of the jacobian respectively - for (i = 0; i < nY; i++) - tmp[i] = static_cast(ll_Incidence.get(ll_row, varOrder[i]-1)); + for (int i = 0; i < nY; i++) + tmp[i] = static_cast(ll_Incidence.get(ll_row, var_order[i]-1)); // write the reordered blocks back to the jacobian // in reverse order - for (j = nY-1; j >= 0; j--) + for (int j = nY-1; j >= 0; j--) if (tmp[j]) { - JacobianIndices[rjoff] = tmp[j] -1; + dynppToDyn[rjoff] = tmp[j]-1; rjoff--; if (rjoff < 0) break; @@ -318,13 +248,15 @@ KordpDynare::ReorderDynareJacobianIndices() noexcept(false) } //add the indices for the nExog exogenous jacobians - for (j = nJcols-nExog; j < nJcols; j++) - JacobianIndices[j] = j; + for (int j = nJcols-nExog; j < nJcols; j++) + dynppToDyn[j] = j; + + // Compute reverse ordering + dynToDynpp.resize(nJcols); + for (int i = 0; i < nJcols; i++) + dynToDynpp[dynppToDyn[i]] = i; } -/**************************************************************************************/ -/* DynareNameList class */ -/**************************************************************************************/ DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector names_arg) : names(std::move(names_arg)) @@ -334,7 +266,7 @@ DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector +#include #include +#include "journal.hh" +#include "Vector.hh" +#include "twod_matrix.hh" #include "t_container.hh" #include "sparse_tensor.hh" #include "decision_rule.hh" #include "dynamic_model.hh" - -#include "exception.hh" -#include "dynare_exception.hh" -#include "fs_tensor.hh" -#include "SylvException.hh" -#include "tl_exception.hh" -#include "kord_exception.hh" -#include "nlsolve.hh" -#include "approximation.hh" +#include "dynamic_abstract_class.hh" class KordpDynare; -/*////////////////////////////////////////////*/ -// instantiations of pure abstract class NameList in dynamic_model.h: -/*////////////////////////////////////////////*/ +// instantiations of pure abstract class NameList in dynamic_model.h class DynareNameList : public NameList { std::vector names; public: DynareNameList(const KordpDynare &dynare, std::vector names_arg); int - getNum() const + getNum() const override { return static_cast(names.size()); } const std::string & - getName(int i) const + getName(int i) const override { return names[i]; } @@ -66,44 +60,37 @@ public: DynareStateNameList(const KordpDynare &dynare, const DynareNameList &dnl, const DynareNameList &denl); int - getNum() const + getNum() const override { return static_cast(names.size()); } const std::string & - getName(int i) const + getName(int i) const override { return names[i]; } }; -/*********************************************/ + // The following only implements DynamicModel with help of ogdyn::DynareModel // instantiation of pure abstract DynamicModel decl. in dynamic_model.h -class DynamicModelAC; -class DynamicModelDLL; -class DynamicModelMFile; - class KordpDynare : public DynamicModel { - friend class DynareNameList; - friend class DynareStateNameList; - friend class DynamicModelDLL; - friend class DynamicModelMFile; - +public: const int nStat; const int nBoth; const int nPred; const int nForw; const int nExog; const int nPar; - const int nYs; // ={npred + nboth ; } - const int nYss; // nyss ={ nboth + nforw ; } - const int nY; // = num_endo={ nstat + npred + nboth + nforw ; } - const int nJcols; // no of jacobian columns= nExog+nEndo+nsPred+nsForw - const Vector &NNZD; /* the total number of non-zero derivative elements - where hessian is 2nd : NZZD(order=2) */ + const int nYs; // = npred + nboth + const int nYss; // = nboth + nforw + const int nY; // = nstat + npred + nboth + nforw + const int nJcols; // nb of jacobian columns = nExog+nY+nYs+nYss + const Vector &NNZD; /* the total number of non-zero derivative elements + where hessian is 2nd : NNZD(order=2) */ const int nSteps; const int nOrder; +private: Journal &journal; Vector &ySteady; Vector ¶ms; @@ -111,141 +98,110 @@ class KordpDynare : public DynamicModel TensorContainer md; // ModelDerivatives DynareNameList dnl, denl; DynareStateNameList dsnl; - const double ss_tol; - const std::vector &varOrder; const TwoDMatrix &ll_Incidence; - double qz_criterium; - std::vector JacobianIndices; + 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 g1p, g2p, g3p; + + std::unique_ptr dynamicModelFile; public: - KordpDynare(const std::vector &endo, int num_endo, + KordpDynare(const std::vector &endo, const std::vector &exo, int num_exo, int num_par, Vector &ySteady, TwoDMatrix &vCov, Vector ¶ms, int nstat, int nPred, - int nforw, int nboth, const int nJcols, const Vector &NNZD, - const int nSteps, const int ord, - Journal &jr, std::unique_ptr dynamicModelFile_arg, double sstol, + int nforw, int nboth, const Vector &NNZD, + int nSteps, int ord, + Journal &jr, std::unique_ptr dynamicModelFile_arg, const std::vector &varOrder, const TwoDMatrix &ll_Incidence, - double qz_criterium) noexcept(false); - KordpDynare(const std::vector &endo, int num_endo, - const std::vector &exo, int num_exo, int num_par, - Vector &ySteady, TwoDMatrix &vCov, Vector ¶ms, int nstat, int nPred, - int nforw, int nboth, const int nJcols, const Vector &NNZD, - const int nSteps, const int ord, - Journal &jr, std::unique_ptr dynamicModelFile_arg, double sstol, - const std::vector &varOrder, const TwoDMatrix &ll_Incidence, - double qz_criterium, std::unique_ptr g1_arg, - std::unique_ptr g2_arg, std::unique_ptr g3_arg) noexcept(false); + std::unique_ptr g1_arg, + std::unique_ptr g2_arg, std::unique_ptr g3_arg); int - nstat() const + nstat() const override { return nStat; } int - nboth() const + nboth() const override { return nBoth; } int - npred() const + npred() const override { return nPred; } int - nforw() const + nforw() const override { return nForw; } int - nexog() const + nexog() const override { return nExog; } int - nys() const - { - return nYs; - } - int - nyss() const - { - return nYss; - } - int - ny() const - { - return nY; - } - int - steps() const - { - return nSteps; - } - int - order() const + order() const override { return nOrder; } const NameList & - getAllEndoNames() const + getAllEndoNames() const override { return dnl; } const NameList & - getStateNames() const + getStateNames() const override { return dsnl; } const NameList & - getExogNames() const + getExogNames() const override { return denl; } const TwoDMatrix & - getVcov() const + getVcov() const override { return vCov; } - Vector & - getParams() - { - return params; - } const TensorContainer & - getModelDerivatives() const + getModelDerivatives() const override { return md; } const Vector & - getSteady() const + getSteady() const override { return ySteady; } Vector & - getSteady() + getSteady() override { return ySteady; } - void solveDeterministicSteady(); - void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false); + void solveDeterministicSteady() override; + void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) override; void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, - const ConstVector &yyp, const Vector &xx) noexcept(false); - void calcDerivativesAtSteady(); - std::unique_ptr dynamicModelFile; + const ConstVector &yyp, const Vector &xx) override; + void calcDerivativesAtSteady() override; std::unique_ptr - clone() const + clone() const override { std::cerr << "KordpDynare::clone() not implemented" << std::endl; exit(EXIT_FAILURE); } - void LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false); // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags - private: - void ReorderDynareJacobianIndices() noexcept(false); - void populateDerivativesContainer(const TwoDMatrix &g, int ord, const std::vector &vOrder); + // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags + void LLxSteady(const Vector &yS, Vector &llxSteady); + /* Computes the permutations mapping back and forth between Dynare and + Dynare++ orderings of variables */ + void computeJacobianPermutation(const std::vector &var_order); + // Fills tensor container (with reordering) at a given order + void populateDerivativesContainer(const TwoDMatrix &g, int ord); }; #endif diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc index dec222d11..19fc3bbfe 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cc +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc @@ -17,36 +17,33 @@ * along with Dynare. If not, see . */ -/* - Defines the entry point for the k-order perturbation application DLL. +/* Defines the entry point for the k-order perturbation application DLL. - Inputs: - 1) dr - 2) M_ - 3) options - - Outputs: - - if order == 1: only g_1 - - if order == 2: g_0, g_1, g_2 - - if order == 3: g_0, g_1, g_2, g_3 + See matlab/mex/k_order_perturbation.m for a description of inputs and + outputs. */ #include "dynamic_m.hh" #include "dynamic_dll.hh" +#include "k_ord_dynare.hh" + +#include "approximation.hh" +#include "exception.hh" +#include "dynare_exception.hh" +#include "kord_exception.hh" +#include "tl_exception.hh" +#include "SylvException.hh" #include -#include #include #include "dynmex.h" -////////////////////////////////////////////////////// -// Convert MATLAB Dynare endo and exo names array to a vector array of string pointers -// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows -// hence a rather low level approach is needed -/////////////////////////////////////////////////////// +/* Convert MATLAB Dynare endo and exo names array to a vector array of + string pointers. MATLAB “mx” function returns a long string concatenated by + columns rather than rows hence a rather low level approach is needed. */ void -DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, std::vector &out) +DynareMxArrayToString(const mxArray *mxFldp, int len, int width, std::vector &out) { char *cNamesCharStr = mxArrayToString(mxFldp); @@ -55,7 +52,7 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, std for (int i = 0; i < width; i++) for (int j = 0; j < len; j++) // Allow alphanumeric and underscores "_" only: - if (isalnum(cNamesCharStr[j+i*len]) || (cNamesCharStr[j+i*len] == '_')) + if (std::isalnum(cNamesCharStr[j+i*len]) || (cNamesCharStr[j+i*len] == '_')) out[j] += cNamesCharStr[j+i*len]; } @@ -122,14 +119,10 @@ extern "C" { const int nStat = static_cast(mxGetScalar(mxFldp)); mxFldp = mxGetField(M_, 0, "npred"); const int nPred = static_cast(mxGetScalar(mxFldp)); - mxFldp = mxGetField(M_, 0, "nspred"); - const int nsPred = static_cast(mxGetScalar(mxFldp)); mxFldp = mxGetField(M_, 0, "nboth"); const int nBoth = static_cast(mxGetScalar(mxFldp)); mxFldp = mxGetField(M_, 0, "nfwrd"); const int nForw = static_cast(mxGetScalar(mxFldp)); - mxFldp = mxGetField(M_, 0, "nsfwrd"); - const int nsForw = static_cast(mxGetScalar(mxFldp)); mxFldp = mxGetField(M_, 0, "exo_nbr"); const int nExog = static_cast(mxGetScalar(mxFldp)); @@ -158,13 +151,10 @@ extern "C" { DYN_MEX_FUNC_ERR_MSG_TXT(("dynare:k_order_perturbation: Incorrect length of lead lag incidences: ncol=" + std::to_string(npar) + " != nEndo=" + std::to_string(nEndo)).c_str()); - //get NNZH =NNZD(2) = the total number of non-zero Hessian elements mxFldp = mxGetField(M_, 0, "NNZDerivatives"); Vector NNZD{mxFldp}; if (NNZD[kOrder-1] == -1) - DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the stoch_simul command"); - - const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns + DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the 'stoch_simul' command"); mxFldp = mxGetField(M_, 0, "var_order_endo_names"); const int nendo = static_cast(mxGetM(mxFldp)); @@ -182,9 +172,9 @@ extern "C" { DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input parameters."); std::unique_ptr g1m, g2m, g3m; - // derivatives passed as arguments */ if (nrhs > 3) { + // Derivatives have been passed as arguments const mxArray *g1 = prhs[3]; int m = static_cast(mxGetM(g1)); int n = static_cast(mxGetN(g1)); @@ -206,7 +196,6 @@ extern "C" { } const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state - const double sstol = 1.e-13; //NL solver tolerance from try { @@ -222,10 +211,10 @@ extern "C" { TLStatic::init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog); // make KordpDynare object - KordpDynare dynare(endoNames, nEndo, exoNames, nExog, nPar, + KordpDynare dynare(endoNames, exoNames, nExog, nPar, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, - jcols, NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile), - sstol, var_order_vp, llincidence, qz_criterium, + NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile), + var_order_vp, llincidence, std::move(g1m), std::move(g2m), std::move(g3m)); // construct main K-order approximation class