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

k-order DLL: make model derivatives storage ready for an arbitrary order

time-shift
Sébastien Villemot 2019-04-10 18:56:14 +02:00
parent 85f472c74b
commit 2dc988b69d
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
8 changed files with 71 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
mex/sources/k_order_perturbation/dynamic_abstract_class.hh
View File

@ -20,6 +20,8 @@
#ifndef _DYNAMICMODELAC_HH
#define _DYNAMICMODELAC_HH
#include <vector>
#include "twod_matrix.hh"
class DynamicModelAC
@ -27,6 +29,6 @@ class DynamicModelAC
public:
virtual ~DynamicModelAC() = default;
virtual void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) = 0;
Vector &residual, std::vector<TwoDMatrix> &md) = 0;
};
#endif

21
mex/sources/k_order_perturbation/dynamic_dll.cc
View File

@ -94,21 +94,20 @@ DynamicModelDLL::~DynamicModelDLL()
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)
Vector &residual, std::vector<TwoDMatrix> &md) noexcept(false)
{
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, tt.get());
if (g1)
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, tt.get());
if (g2)
dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g2->base());
if (g3)
dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[0].base());
if (md.size() >= 2)
{
dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[1].base());
}
if (md.size() >= 3)
{
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());
dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[2].base());
}
}

2
mex/sources/k_order_perturbation/dynamic_dll.hh
View File

@ -66,6 +66,6 @@ public:
virtual ~DynamicModelDLL();
void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override;
Vector &residual, std::vector<TwoDMatrix> &md) override;
};
#endif

34
mex/sources/k_order_perturbation/dynamic_m.cc
View File

@ -30,15 +30,15 @@ DynamicModelMFile::DynamicModelMFile(const std::string &modName) noexcept(false)
}
void
DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm)
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);
assert(tdm.ncols() == 3);
assert(tdm.nrows() == sizeRowIdxVector);
double *ptr = mxGetPr(sparseMat);
@ -48,9 +48,9 @@ DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat,
for (int i = 0; i < totalCols; i++)
for (int j = 0; j < static_cast<int>((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];
tdm.get(output_row, 0) = rowIdxVector[rind] + 1;
tdm.get(output_row, 1) = i + 1;
tdm.get(output_row, 2) = ptr[rind];
output_row++;
}
@ -60,16 +60,16 @@ DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat,
recognized as such by KordpDynare::populateDerivativesContainer() */
while (output_row < static_cast<int>(sizeRowIdxVector))
{
tdm->get(output_row, 0) = 0;
tdm->get(output_row, 1) = 0;
tdm->get(output_row, 2) = 0;
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)
Vector &residual, std::vector<TwoDMatrix> &md) noexcept(false)
{
constexpr int nlhs_dynamic = 4, nrhs_dynamic = 5;
mxArray *prhs[nrhs_dynamic], *plhs[nlhs_dynamic];
@ -91,14 +91,14 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
residual = Vector{plhs[0]};
assert(static_cast<int>(mxGetM(plhs[1])) == g1->nrows());
assert(static_cast<int>(mxGetN(plhs[1])) == g1->ncols());
std::copy_n(mxGetPr(plhs[1]), mxGetM(plhs[1])*mxGetN(plhs[1]), g1->base());
assert(static_cast<int>(mxGetM(plhs[1])) == md[0].nrows());
assert(static_cast<int>(mxGetN(plhs[1])) == md[0].ncols());
std::copy_n(mxGetPr(plhs[1]), mxGetM(plhs[1])*mxGetN(plhs[1]), md[0].base());
if (g2)
unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2);
if (g3)
unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[3], g3);
if (md.size() >= 2)
unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], md[1]);
if (md.size() >= 3)
unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[3], md[2]);
for (int i = 0; i < nrhs_dynamic; i++)
mxDestroyArray(prhs[i]);

4
mex/sources/k_order_perturbation/dynamic_m.hh
View File

@ -33,11 +33,11 @@ class DynamicModelMFile : public DynamicModelAC
{
private:
const std::string DynamicMFilename;
static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm);
static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix &tdm);
public:
explicit DynamicModelMFile(const std::string &modName);
virtual ~DynamicModelMFile() = default;
void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override;
Vector &residual, std::vector<TwoDMatrix> &md) override;
};
#endif

46
mex/sources/k_order_perturbation/k_ord_dynare.cc
View File

@ -31,16 +31,13 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo,
int npred, int nforw, int nboth, const Vector &nnzd,
int nsteps, int norder,
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
const std::vector<int> &dr_order, const TwoDMatrix &llincidence,
std::unique_ptr<TwoDMatrix> g1_arg, std::unique_ptr<TwoDMatrix> g2_arg,
std::unique_ptr<TwoDMatrix> g3_arg) :
const std::vector<int> &dr_order, const TwoDMatrix &llincidence) :
nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
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},
ll_Incidence{llincidence},
g1p{std::move(g1_arg)}, g2p{std::move(g2_arg)}, g3p{std::move(g3_arg)},
dynamicModelFile{std::move(dynamicModelFile_arg)}
{
computeJacobianPermutation(dr_order);
@ -71,22 +68,16 @@ KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVect
void
KordpDynare::calcDerivativesAtSteady()
{
if (!g1p)
if (dyn_md.empty())
{
g1p = std::make_unique<TwoDMatrix>(nY, nJcols);
g1p->zeros();
dyn_md.emplace_back(nY, nJcols); // Allocate Jacobian
dyn_md.back().zeros();
if (nOrder > 1)
for (int i = 2; i <= nOrder; i++)
{
// allocate space for sparse Hessian
g2p = std::make_unique<TwoDMatrix>(static_cast<int>(NNZD[1]), 3);
g2p->zeros();
}
if (nOrder > 2)
{
g3p = std::make_unique<TwoDMatrix>(static_cast<int>(NNZD[2]), 3);
g3p->zeros();
// Higher order derivatives, as sparse (3-column) matrices
dyn_md.emplace_back(static_cast<int>(NNZD[i-1]), 3);
dyn_md.back().zeros();
}
Vector xx(nexog());
@ -97,21 +88,18 @@ KordpDynare::calcDerivativesAtSteady()
Vector llxSteady(nJcols-nExog);
LLxSteady(ySteady, llxSteady);
dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p.get(), g2p.get(), g3p.get());
dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, dyn_md);
}
populateDerivativesContainer(*g1p, 1);
if (nOrder > 1)
populateDerivativesContainer(*g2p, 2);
if (nOrder > 2)
populateDerivativesContainer(*g3p, 3);
for (int i = 1; i <= nOrder; i++)
populateDerivativesContainer(i);
}
void
KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord)
KordpDynare::populateDerivativesContainer(int ord)
{
const TwoDMatrix &g = dyn_md[ord-1];
// model derivatives FSSparseTensor instance
auto mdTi = std::make_unique<FSSparseTensor>(ord, nJcols, nY);
@ -249,6 +237,11 @@ KordpDynare::computeJacobianPermutation(const std::vector<int> &dr_order)
dynToDynpp[dynppToDyn[i]] = i;
}
void
KordpDynare::push_back_md(const mxArray *m)
{
dyn_md.emplace_back(ConstTwoDMatrix{m});
}
DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg)
: names(std::move(names_arg))
@ -263,3 +256,4 @@ DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const Dynare
for (int i = 0; i < dynare.nexog(); i++)
names.emplace_back(denl.getName(i));
}

15
mex/sources/k_order_perturbation/k_ord_dynare.hh
View File

@ -95,15 +95,14 @@ private:
Vector &ySteady;
Vector &params;
TwoDMatrix &vCov;
TensorContainer<FSSparseTensor> md; // ModelDerivatives
std::vector<TwoDMatrix> dyn_md; // Model derivatives, in Dynare form
TensorContainer<FSSparseTensor> md; // Model derivatives, in Dynare++ form
DynareNameList dnl, denl;
DynareStateNameList dsnl;
const TwoDMatrix &ll_Incidence;
std::vector<int> dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones
std::vector<int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
std::unique_ptr<TwoDMatrix> g1p, g2p, g3p;
std::unique_ptr<DynamicModelAC> dynamicModelFile;
public:
KordpDynare(const std::vector<std::string> &endo,
@ -112,9 +111,7 @@ public:
int nforw, int nboth, const Vector &NNZD,
int nSteps, int ord,
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence,
std::unique_ptr<TwoDMatrix> g1_arg,
std::unique_ptr<TwoDMatrix> g2_arg, std::unique_ptr<TwoDMatrix> g3_arg);
const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence);
int
nstat() const override
@ -194,14 +191,16 @@ public:
std::cerr << "KordpDynare::clone() not implemented" << std::endl;
exit(EXIT_FAILURE);
}
// Add model derivatives of a given order passed as argument to the MEX
void push_back_md(const mxArray *m);
private:
// 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<int> &var_order);
// Fills tensor container (with reordering) at a given order
void populateDerivativesContainer(const TwoDMatrix &g, int ord);
// Fills model derivatives in Dynare++ form (at a given order) given the Dynare form
void populateDerivativesContainer(int ord);
};
#endif

38
mex/sources/k_order_perturbation/k_order_perturbation.cc
View File

@ -163,30 +163,6 @@ extern "C" {
if (nEndo != static_cast<int>(endoNames.size()) || nExog != static_cast<int>(exoNames.size()))
DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of M_.endo_names or M_.exo_names");
std::unique_ptr<TwoDMatrix> g1m, g2m, g3m;
if (nrhs > 3)
{
// Derivatives have been passed as arguments
const mxArray *g1 = prhs[3];
int m = static_cast<int>(mxGetM(g1));
int n = static_cast<int>(mxGetN(g1));
g1m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g1}});
if (nrhs > 4)
{
const mxArray *g2 = prhs[4];
int m = static_cast<int>(mxGetM(g2));
int n = static_cast<int>(mxGetN(g2));
g2m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g2}});
if (nrhs > 5)
{
const mxArray *g3 = prhs[5];
int m = static_cast<int>(mxGetM(g3));
int n = static_cast<int>(mxGetN(g3));
g3m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g3}});
}
}
}
const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
try
@ -206,11 +182,19 @@ extern "C" {
KordpDynare dynare(endoNames, exoNames, nExog, nPar,
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
dr_order, llincidence,
std::move(g1m), std::move(g2m), std::move(g3m));
dr_order, llincidence);
// If model derivatives have been passed as arguments
if (nrhs > 3)
{
dynare.push_back_md(prhs[3]);
if (nrhs > 4)
dynare.push_back_md(prhs[4]);
if (nrhs > 5)
dynare.push_back_md(prhs[5]);
}
// construct main K-order approximation class
Approximation app(dynare, journal, nSteps, false, qz_criterium);
// run stochastic steady
app.walkStochSteady();