k_order_perturbation MEX: improve input handling and sanitization
parent
f8af21819e
commit
aaaa897a2a
|
@ -28,10 +28,10 @@
|
||||||
KordpDynare::KordpDynare(const std::vector<std::string> &endo,
|
KordpDynare::KordpDynare(const std::vector<std::string> &endo,
|
||||||
const std::vector<std::string> &exo, int nexog, int npar,
|
const std::vector<std::string> &exo, int nexog, int npar,
|
||||||
Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
|
Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
|
||||||
int npred, int nforw, int nboth, const Vector &nnzd,
|
int npred, int nforw, int nboth, const ConstVector &nnzd,
|
||||||
int nsteps, int norder,
|
int nsteps, int norder,
|
||||||
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
|
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
|
||||||
const std::vector<int> &dr_order, const TwoDMatrix &llincidence) :
|
const std::vector<int> &dr_order, const ConstTwoDMatrix &llincidence) :
|
||||||
nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
|
nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
|
||||||
nYs{npred + nboth}, nYss{nboth + nforw}, nY{nstat + npred + nboth + nforw},
|
nYs{npred + nboth}, nYss{nboth + nforw}, nY{nstat + npred + nboth + nforw},
|
||||||
nJcols{nExog+nY+nYs+nYss}, NNZD{nnzd}, nSteps{nsteps},
|
nJcols{nExog+nY+nYs+nYss}, NNZD{nnzd}, nSteps{nsteps},
|
||||||
|
|
|
@ -86,7 +86,7 @@ public:
|
||||||
const int nYss; // = nboth + nforw
|
const int nYss; // = nboth + nforw
|
||||||
const int nY; // = nstat + npred + nboth + nforw
|
const int nY; // = nstat + npred + nboth + nforw
|
||||||
const int nJcols; // nb of jacobian columns = nExog+nY+nYs+nYss
|
const int nJcols; // nb of jacobian columns = nExog+nY+nYs+nYss
|
||||||
const Vector &NNZD; /* the total number of non-zero derivative elements
|
const ConstVector &NNZD; /* the total number of non-zero derivative elements
|
||||||
where hessian is 2nd : NNZD(order=2) */
|
where hessian is 2nd : NNZD(order=2) */
|
||||||
const int nSteps;
|
const int nSteps;
|
||||||
const int nOrder;
|
const int nOrder;
|
||||||
|
@ -99,7 +99,7 @@ private:
|
||||||
TensorContainer<FSSparseTensor> md; // Model derivatives, in Dynare++ form
|
TensorContainer<FSSparseTensor> md; // Model derivatives, in Dynare++ form
|
||||||
DynareNameList dnl, denl;
|
DynareNameList dnl, denl;
|
||||||
DynareStateNameList dsnl;
|
DynareStateNameList dsnl;
|
||||||
const TwoDMatrix &ll_Incidence;
|
const ConstTwoDMatrix &ll_Incidence;
|
||||||
std::vector<int> dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones
|
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::vector<int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
|
||||||
|
|
||||||
|
@ -108,10 +108,10 @@ public:
|
||||||
KordpDynare(const std::vector<std::string> &endo,
|
KordpDynare(const std::vector<std::string> &endo,
|
||||||
const std::vector<std::string> &exo, int num_exo, int num_par,
|
const std::vector<std::string> &exo, int num_exo, int num_par,
|
||||||
Vector &ySteady, TwoDMatrix &vCov, Vector ¶ms, int nstat, int nPred,
|
Vector &ySteady, TwoDMatrix &vCov, Vector ¶ms, int nstat, int nPred,
|
||||||
int nforw, int nboth, const Vector &NNZD,
|
int nforw, int nboth, const ConstVector &NNZD,
|
||||||
int nSteps, int ord,
|
int nSteps, int ord,
|
||||||
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
|
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
|
||||||
const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence);
|
const std::vector<int> &varOrder, const ConstTwoDMatrix &ll_Incidence);
|
||||||
|
|
||||||
int
|
int
|
||||||
nstat() const override
|
nstat() const override
|
||||||
|
|
|
@ -54,7 +54,12 @@ DynareMxArrayToString(const mxArray *mxFldp)
|
||||||
assert(mxIsCell(mxFldp));
|
assert(mxIsCell(mxFldp));
|
||||||
std::vector<std::string> r;
|
std::vector<std::string> r;
|
||||||
for (size_t i = 0; i < mxGetNumberOfElements(mxFldp); i++)
|
for (size_t i = 0; i < mxGetNumberOfElements(mxFldp); i++)
|
||||||
r.emplace_back(mxArrayToString(mxGetCell(mxFldp, i)));
|
{
|
||||||
|
const mxArray *cell_mx = mxGetCell(mxFldp, i);
|
||||||
|
if (!(cell_mx && mxIsChar(cell_mx)))
|
||||||
|
mexErrMsgTxt("Cell is not a character array");
|
||||||
|
r.emplace_back(mxArrayToString(cell_mx));
|
||||||
|
}
|
||||||
|
|
||||||
return r;
|
return r;
|
||||||
}
|
}
|
||||||
|
@ -77,108 +82,120 @@ extern "C" {
|
||||||
mexFunction(int nlhs, mxArray *plhs[],
|
mexFunction(int nlhs, mxArray *plhs[],
|
||||||
int nrhs, const mxArray *prhs[])
|
int nrhs, const mxArray *prhs[])
|
||||||
{
|
{
|
||||||
if (nrhs < 3 || nlhs < 2)
|
if (nrhs < 3 || nlhs < 2 || nlhs > 3)
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("Must have at least 3 input parameters and takes at least 2 output parameters.");
|
DYN_MEX_FUNC_ERR_MSG_TXT("Must have at least 3 input parameters and takes 2 or 3 output parameters.");
|
||||||
|
|
||||||
const mxArray *dr = prhs[0];
|
// Give explicit names to input arguments
|
||||||
const mxArray *M_ = prhs[1];
|
const mxArray *dr_mx = prhs[0];
|
||||||
const mxArray *options_ = prhs[2];
|
const mxArray *M_mx = prhs[1];
|
||||||
bool use_dll = mxGetScalar(mxGetField(options_, 0, "use_dll")) != 0;
|
const mxArray *options_mx = prhs[2];
|
||||||
|
|
||||||
mxArray *mFname = mxGetField(M_, 0, "fname");
|
auto get_int_field = [](const mxArray *struct_mx, const std::string &fieldname)
|
||||||
if (!mxIsChar(mFname))
|
{
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("Input must be of type char.");
|
mxArray *field_mx = mxGetField(struct_mx, 0, fieldname.c_str());
|
||||||
|
if (!(field_mx && mxIsScalar(field_mx) && mxIsNumeric(field_mx)))
|
||||||
|
mexErrMsgTxt(("Field `" + fieldname + "' should be a numeric scalar").c_str());
|
||||||
|
return static_cast<int>(mxGetScalar(field_mx));
|
||||||
|
};
|
||||||
|
|
||||||
std::string fName = mxArrayToString(mFname);
|
// Extract various fields from options_
|
||||||
|
const int kOrder = get_int_field(options_mx, "order");
|
||||||
int kOrder;
|
|
||||||
mxArray *mxFldp = mxGetField(options_, 0, "order");
|
|
||||||
if (!mxIsNumeric(mxFldp))
|
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("options_.order must be a numeric value");
|
|
||||||
kOrder = static_cast<int>(mxGetScalar(mxFldp));
|
|
||||||
if (kOrder < 1)
|
if (kOrder < 1)
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("options_.order must be at least 1");
|
DYN_MEX_FUNC_ERR_MSG_TXT("options_.order must be at least 1");
|
||||||
|
|
||||||
|
const mxArray *use_dll_mx = mxGetField(options_mx, 0, "use_dll");
|
||||||
|
if (!(use_dll_mx && mxIsLogicalScalar(use_dll_mx)))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("options_.use_dll should be a logical scalar");
|
||||||
|
bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
|
||||||
|
|
||||||
double qz_criterium = 1+1e-6;
|
double qz_criterium = 1+1e-6;
|
||||||
mxFldp = mxGetField(options_, 0, "qz_criterium");
|
const mxArray *qz_criterium_mx = mxGetField(options_mx, 0, "qz_criterium");
|
||||||
if (mxGetNumberOfElements(mxFldp) > 0 && mxIsNumeric(mxFldp))
|
if (qz_criterium_mx && mxIsScalar(qz_criterium_mx) && mxIsNumeric(qz_criterium_mx))
|
||||||
qz_criterium = mxGetScalar(mxFldp);
|
qz_criterium = mxGetScalar(qz_criterium_mx);
|
||||||
|
|
||||||
mxFldp = mxGetField(M_, 0, "params");
|
// Extract various fields from M_
|
||||||
Vector modParams{mxFldp};
|
const mxArray *fname_mx = mxGetField(M_mx, 0, "fname");
|
||||||
|
if (!(fname_mx && mxIsChar(fname_mx) && mxGetM(fname_mx) == 1))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.fname should be a character string");
|
||||||
|
std::string fname{mxArrayToString(fname_mx)};
|
||||||
|
|
||||||
|
const mxArray *params_mx = mxGetField(M_mx, 0, "params");
|
||||||
|
if (!(params_mx && mxIsDouble(params_mx)))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.params should be a double precision array");
|
||||||
|
Vector modParams{ConstVector{params_mx}};
|
||||||
if (!modParams.isFinite())
|
if (!modParams.isFinite())
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("The parameters vector contains NaN or Inf");
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.params contains NaN or Inf");
|
||||||
|
|
||||||
mxFldp = mxGetField(M_, 0, "Sigma_e");
|
const mxArray *sigma_e_mx = mxGetField(M_mx, 0, "Sigma_e");
|
||||||
int dim = static_cast<int>(mxGetN(mxFldp));
|
if (!(sigma_e_mx && mxIsDouble(sigma_e_mx) && mxGetM(sigma_e_mx) == mxGetN(sigma_e_mx)))
|
||||||
TwoDMatrix vCov(dim, dim, Vector{mxFldp});
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.Sigma_e should be a double precision square matrix");
|
||||||
|
TwoDMatrix vCov{ConstTwoDMatrix{sigma_e_mx}};
|
||||||
if (!vCov.isFinite())
|
if (!vCov.isFinite())
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("The covariance matrix of shocks contains NaN or Inf");
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.Sigma_e contains NaN or Inf");
|
||||||
|
|
||||||
mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var
|
const int nStat = get_int_field(M_mx, "nstatic");
|
||||||
Vector ySteady{mxFldp};
|
const int nPred = get_int_field(M_mx, "npred");
|
||||||
if (!ySteady.isFinite())
|
const int nBoth = get_int_field(M_mx, "nboth");
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("The steady state vector contains NaN or Inf");
|
const int nForw = get_int_field(M_mx, "nfwrd");
|
||||||
|
|
||||||
mxFldp = mxGetField(M_, 0, "nstatic");
|
const int nExog = get_int_field(M_mx, "exo_nbr");
|
||||||
const int nStat = static_cast<int>(mxGetScalar(mxFldp));
|
const int nEndo = get_int_field(M_mx, "endo_nbr");
|
||||||
mxFldp = mxGetField(M_, 0, "npred");
|
const int nPar = get_int_field(M_mx, "param_nbr");
|
||||||
const int nPred = static_cast<int>(mxGetScalar(mxFldp));
|
|
||||||
mxFldp = mxGetField(M_, 0, "nboth");
|
|
||||||
const int nBoth = static_cast<int>(mxGetScalar(mxFldp));
|
|
||||||
mxFldp = mxGetField(M_, 0, "nfwrd");
|
|
||||||
const int nForw = static_cast<int>(mxGetScalar(mxFldp));
|
|
||||||
|
|
||||||
mxFldp = mxGetField(M_, 0, "exo_nbr");
|
const mxArray *lead_lag_incidence_mx = mxGetField(M_mx, 0, "lead_lag_incidence");
|
||||||
const int nExog = static_cast<int>(mxGetScalar(mxFldp));
|
if (!(lead_lag_incidence_mx && mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == 3
|
||||||
mxFldp = mxGetField(M_, 0, "endo_nbr");
|
&& mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(nEndo)))
|
||||||
const int nEndo = static_cast<int>(mxGetScalar(mxFldp));
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.lead_lag_incidence should be a double precision matrix with 3 rows and M_.endo_nbr columns");
|
||||||
mxFldp = mxGetField(M_, 0, "param_nbr");
|
ConstTwoDMatrix llincidence{lead_lag_incidence_mx};
|
||||||
const int nPar = static_cast<int>(mxGetScalar(mxFldp));
|
|
||||||
|
|
||||||
mxFldp = mxGetField(dr, 0, "order_var");
|
const mxArray *nnzderivatives_mx = mxGetField(M_mx, 0, "NNZDerivatives");
|
||||||
dim = static_cast<int>(mxGetM(mxFldp));
|
if (!(nnzderivatives_mx && mxIsDouble(nnzderivatives_mx)))
|
||||||
if (dim != nEndo)
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.NNZDerivatives should be a double precision array");
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of dr.order_var");
|
ConstVector NNZD{nnzderivatives_mx};
|
||||||
std::vector<int> dr_order(nEndo);
|
if (NNZD.length() < kOrder || NNZD[kOrder-1] == -1)
|
||||||
std::transform(mxGetPr(mxFldp), mxGetPr(mxFldp)+dim, dr_order.begin(),
|
|
||||||
[](double x) { return static_cast<int>(x)-1; });
|
|
||||||
|
|
||||||
// the lag, current and lead blocks of the jacobian respectively
|
|
||||||
TwoDMatrix llincidence(mxGetField(M_, 0, "lead_lag_incidence"));
|
|
||||||
if (llincidence.nrows() != 3 || llincidence.ncols() != nEndo)
|
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of M_.lead_lag_incidence");
|
|
||||||
|
|
||||||
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");
|
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, "endo_names");
|
const mxArray *endo_names_mx = mxGetField(M_mx, 0, "endo_names");
|
||||||
std::vector<std::string> endoNames = DynareMxArrayToString(mxFldp);
|
if (!(endo_names_mx && mxIsCell(endo_names_mx) && mxGetNumberOfElements(endo_names_mx) == static_cast<size_t>(nEndo)))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.endo_names should be a cell array of M_.endo_nbr elements");
|
||||||
|
std::vector<std::string> endoNames = DynareMxArrayToString(endo_names_mx);
|
||||||
|
|
||||||
mxFldp = mxGetField(M_, 0, "exo_names");
|
const mxArray *exo_names_mx = mxGetField(M_mx, 0, "exo_names");
|
||||||
std::vector<std::string> exoNames = DynareMxArrayToString(mxFldp);
|
if (!(exo_names_mx && mxIsCell(exo_names_mx) && mxGetNumberOfElements(exo_names_mx) == static_cast<size_t>(nExog)))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.exo_names should be a cell array of M_.exo_nbr elements");
|
||||||
|
std::vector<std::string> exoNames = DynareMxArrayToString(exo_names_mx);
|
||||||
|
|
||||||
if (nEndo != static_cast<int>(endoNames.size()) || nExog != static_cast<int>(exoNames.size()))
|
const mxArray *dynamic_tmp_nbr_mx = mxGetField(M_mx, 0, "dynamic_tmp_nbr");
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of M_.endo_names or M_.exo_names");
|
if (!(dynamic_tmp_nbr_mx && mxIsDouble(dynamic_tmp_nbr_mx) && mxGetNumberOfElements(dynamic_tmp_nbr_mx) >= static_cast<size_t>(kOrder+1)))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("M_.dynamic_tmp_nbr should be a double precision array with strictly more elements than the order of derivation");
|
||||||
|
int ntt = std::accumulate(mxGetPr(dynamic_tmp_nbr_mx), mxGetPr(dynamic_tmp_nbr_mx)+kOrder+1, 0);
|
||||||
|
|
||||||
mxFldp = mxGetField(M_, 0, "dynamic_tmp_nbr");
|
// Extract various fields from dr
|
||||||
if (static_cast<int>(mxGetM(mxFldp)) < kOrder+1 || mxGetN(mxFldp) != 1)
|
const mxArray *ys_mx = mxGetField(dr_mx, 0, "ys"); // and not in order of dr.order_var
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of M_.dynamic_tmp_nbr");
|
if (!(ys_mx && mxIsDouble(ys_mx)))
|
||||||
int ntt = std::accumulate(mxGetPr(mxFldp), mxGetPr(mxFldp)+kOrder+1, 0);
|
DYN_MEX_FUNC_ERR_MSG_TXT("dr.ys should be a double precision array");
|
||||||
|
Vector ySteady{ConstVector{ys_mx}};
|
||||||
|
if (!ySteady.isFinite())
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("dr.ys contains NaN or Inf");
|
||||||
|
|
||||||
|
const mxArray *order_var_mx = mxGetField(dr_mx, 0, "order_var");
|
||||||
|
if (!(order_var_mx && mxIsDouble(order_var_mx) && mxGetNumberOfElements(order_var_mx) == static_cast<size_t>(nEndo)))
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("dr.order_var should be a double precision array of M_.endo_nbr elements");
|
||||||
|
std::vector<int> dr_order(nEndo);
|
||||||
|
std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx)+nEndo, dr_order.begin(),
|
||||||
|
[](double x) { return static_cast<int>(x)-1; });
|
||||||
|
|
||||||
const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
|
const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
|
||||||
|
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
Journal journal(fName + ".jnl");
|
Journal journal(fname + ".jnl");
|
||||||
|
|
||||||
std::unique_ptr<DynamicModelAC> dynamicModelFile;
|
std::unique_ptr<DynamicModelAC> dynamicModelFile;
|
||||||
if (use_dll)
|
if (use_dll)
|
||||||
dynamicModelFile = std::make_unique<DynamicModelDLL>(fName, ntt, kOrder);
|
dynamicModelFile = std::make_unique<DynamicModelDLL>(fname, ntt, kOrder);
|
||||||
else
|
else
|
||||||
dynamicModelFile = std::make_unique<DynamicModelMFile>(fName, ntt);
|
dynamicModelFile = std::make_unique<DynamicModelMFile>(fname, ntt);
|
||||||
|
|
||||||
// intiate tensor library
|
// intiate tensor library
|
||||||
TLStatic::init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog);
|
TLStatic::init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog);
|
||||||
|
|
Loading…
Reference in New Issue