diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc index 9735e5918..f02ba300d 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cc +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2008-2009 Dynare Team + * Copyright (C) 2008-2010 Dynare Team * * This file is part of Dynare. * @@ -35,54 +35,27 @@ KordpDynare::KordpDynare(const char **endo, int num_endo, const char **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, + 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, DynamicModelDLL &dynamicDLL, double sstol, - const vector *var_order, const TwoDMatrix *llincidence, double criterium) throw (TLException) : + const vector &var_order, const TwoDMatrix &llincidence, double criterium) throw (TLException) : 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), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params(inParams), - md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder(var_order), + md(1), dnl(*this, endo), denl(*this, exo), dsnl(*this, dnl, denl), ss_tol(sstol), varOrder(var_order), ll_Incidence(llincidence), qz_criterium(criterium) { - dnl = new DynareNameList(*this, endo); - denl = new DynareExogNameList(*this, exo); - dsnl = new DynareStateNameList(*this, *dnl, *denl); - - JacobianIndices = ReorderDynareJacobianIndices(varOrder); + ReorderDynareJacobianIndices(); // Initialise ModelDerivativeContainer(*this, this->md, nOrder); for (int iord = 1; iord <= nOrder; iord++) - { - FSSparseTensor *t = new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY); - md.insert(t); - } - + md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY)); } KordpDynare::~KordpDynare() { - if (ySteady) - delete ySteady; - if (params) - delete params; - if (vCov) - delete vCov; - if (dnl) - delete dnl; - if (dsnl) - delete dsnl; - if (denl) - delete denl; - if (JacobianIndices) - delete JacobianIndices; - if (varOrder) - delete varOrder; - if (ll_Incidence) - delete ll_Incidence; - if (NNZD) - delete NNZD; + // No need to manually delete tensors in "md", they are deleted by the TensorContainer destructor } void @@ -111,70 +84,56 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, * this is main derivative calculation functin that indirectly calls dynamic.dll * which performs actual calculation and reorders ***************************************************/ -void -KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException) -{ - TwoDMatrix *g2 = 0; // NULL; - TwoDMatrix *g3 = 0; // NULL; - TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian - g1->zeros(); - - if ((nJcols != g1->ncols()) || (nY != g1->nrows())) - throw DynareException(__FILE__, __LINE__, "Error in calcDerivatives: Created wrong jacobian"); - - if (nOrder > 1) - { - // generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1] - g2 = new TwoDMatrix((int)(*NNZD)[1], 3); - g2->zeros(); - } - if (nOrder > 2) - { - g3 = new TwoDMatrix((int)(*NNZD)[2], 3); - g3->zeros(); - } - Vector out(nY); - out.zeros(); - const Vector *llxYYp; // getting around the constantness - if ((nJcols - nExog) > yy.length()) - llxYYp = (LLxSteady(yy)); - else - llxYYp = &yy; - const Vector &llxYY = *(llxYYp); - - dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3); - - if ((nJcols != g1->ncols()) || (nY != g1->nrows())) - throw DynareException(__FILE__, __LINE__, "Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian"); - - populateDerivativesContainer(g1, 1, JacobianIndices); - delete g1; - if (nOrder > 1) - { - populateDerivativesContainer(g2, 2, JacobianIndices); - delete g2; - } - if (nOrder > 2) - { - populateDerivativesContainer(g3, 3, JacobianIndices); - delete g3; - } - delete llxYYp; -} - void KordpDynare::calcDerivativesAtSteady() throw (DynareException) { + TwoDMatrix g1(nY, nJcols); + g1.zeros(); + + TwoDMatrix *g2p = NULL, *g3p = NULL; + + if (nOrder > 1) + { + // allocate space for sparse Hessian + g2p = new TwoDMatrix((int) NNZD[1], 3); + g2p->zeros(); + } + + if (nOrder > 2) + { + g3p = new TwoDMatrix((int) NNZD[2], 3); + g3p->zeros(); + } + Vector xx(nexog()); xx.zeros(); - calcDerivatives(*ySteady, xx); + + Vector out(nY); + out.zeros(); + Vector llxSteady(nJcols-nExog); + LLxSteady(ySteady, llxSteady); + + dynamicDLL.eval(llxSteady, xx, ¶ms, out, &g1, g2p, g3p); + + populateDerivativesContainer(g1, 1, JacobianIndices); + + if (nOrder > 1) + { + populateDerivativesContainer(*g2p, 2, JacobianIndices); + delete g2p; + } + if (nOrder > 2) + { + populateDerivativesContainer(*g3p, 3, JacobianIndices); + delete g3p; + } } /******************************************************************************* * populateDerivatives to sparse Tensor and fit it in the Derivatives Container *******************************************************************************/ void -KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector *vOrder) +KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const vector &vOrder) { // model derivatives FSSparseTensor instance FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY)); @@ -183,15 +142,15 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vectorncols(); i++) + for (int i = 0; i < g.ncols(); i++) { - for (int j = 0; j < g->nrows(); j++) + for (int j = 0; j < g.nrows(); j++) { double x; if (s[0] < nJcols-nExog) - x = g->get(j, (*vOrder)[s[0]]); + x = g.get(j, vOrder[s[0]]); else - x = g->get(j, s[0]); + x = g.get(j, s[0]); if (x != 0.0) mdTi->insert(s, j, x); } @@ -203,11 +162,11 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector revOrder(nJcols1); for (int i = 0; i < nJcols1; i++) - revOrder[(*vOrder)[i]] = i; - for (int i = 0; i < g->nrows(); i++) + revOrder[vOrder[i]] = i; + for (int i = 0; i < g.nrows(); i++) { - int j = (int) g->get(i, 0)-1; // hessian indices start with 1 - int i1 = (int) g->get(i, 1) -1; + int j = (int) g.get(i, 0)-1; // hessian indices start with 1 + int i1 = (int) g.get(i, 1) -1; int s0 = (int) floor(((double) i1)/((double) nJcols)); int s1 = i1- (nJcols*s0); if (s0 < nJcols1) @@ -220,7 +179,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector= s[0]) { - double x = g->get(i, 2); + double x = g.get(i, 2); mdTi->insert(s, j, x); } } @@ -231,11 +190,11 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector revOrder(nJcols1); for (int i = 0; i < nJcols1; i++) - revOrder[(*vOrder)[i]] = i; - for (int i = 0; i < g->nrows(); i++) + revOrder[vOrder[i]] = i; + for (int i = 0; i < g.nrows(); i++) { - int j = (int) g->get(i, 0)-1; - int i1 = (int) g->get(i, 1) -1; + int j = (int) g.get(i, 0)-1; + int i1 = (int) g.get(i, 1) -1; int s0 = (int) floor(((double) i1)/((double) nJcols2)); int i2 = i1 - nJcols2*s0; int s1 = (int) floor(((double) i2)/((double) nJcols)); @@ -254,7 +213,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector= s[1]) && (s[1] >= s[0])) { - double x = g->get(i, 2); + double x = g.get(i, 2); mdTi->insert(s, j, x); } } @@ -263,6 +222,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector_dynamic DLL *************************************************************/ -Vector * -KordpDynare::LLxSteady(const Vector &yS) throw (DynareException, TLException) +void +KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) throw (DynareException, TLException) { if ((nJcols-nExog) == yS.length()) 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 - Vector *llxSteady = new Vector(nJcols-nExog); - for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) + 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)[((int) ll_Incidence->get(ll_row, i))-1] = yS[i]; + if (ll_Incidence.get(ll_row, i)) + llxSteady[((int) ll_Incidence.get(ll_row, i))-1] = yS[i]; } } - - return llxSteady; } /************************************ @@ -317,27 +277,27 @@ KordpDynare::LLxSteady(const Vector &yS) throw (DynareException, TLException) exogen ************************************/ -vector * -KordpDynare::ReorderDynareJacobianIndices(const vector *varOrder) throw (TLException) +void +KordpDynare::ReorderDynareJacobianIndices() throw (TLException) { // create temporary square 2D matrix size nEndo x nEndo (sparse) // for the lag, current and lead blocks of the jacobian - vector *JacobianIndices = new vector(nJcols); + JacobianIndices.resize(nJcols); vector tmp(nY); int i, j, rjoff = nJcols-nExog-1; - for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) + 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] = ((int) ll_Incidence->get(ll_row, (*varOrder)[i]-1)); + tmp[i] = ((int) ll_Incidence.get(ll_row, varOrder[i]-1)); // write the reordered blocks back to the jacobian // in reverse order for (j = nY-1; j >= 0; j--) if (tmp[j]) { - (*JacobianIndices)[rjoff] = tmp[j] -1; + JacobianIndices[rjoff] = tmp[j] -1; rjoff--; if (rjoff < 0) break; @@ -346,9 +306,7 @@ KordpDynare::ReorderDynareJacobianIndices(const vector *varOrder) throw (TL //add the indices for the nExog exogenous jacobians for (j = nJcols-nExog; j < nJcols; j++) - (*JacobianIndices)[j] = j; - - return JacobianIndices; + JacobianIndices[j] = j; } /**************************************************************************************/ @@ -375,7 +333,7 @@ DynareNameList::selectIndices(const vector &ns) const throw (Dynar DynareNameList::DynareNameList(const KordpDynare &dynare) { for (int i = 0; i < dynare.ny(); i++) - names.push_back(dynare.dnl->getName(i)); + names.push_back(dynare.dnl.getName(i)); } DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp) { @@ -386,7 +344,7 @@ DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp) DynareExogNameList::DynareExogNameList(const KordpDynare &dynare) { for (int i = 0; i < dynare.nexog(); i++) - names.push_back(dynare.denl->getName(i)); + names.push_back(dynare.denl.getName(i)); } DynareExogNameList::DynareExogNameList(const KordpDynare &dynare, const char **namesp) diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh index 1f1ee9139..10ad4f261 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.hh +++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh @@ -122,31 +122,31 @@ class KordpDynare : public DynamicModel 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 + const Vector &NNZD; /* the total number of non-zero derivative elements where hessian is 2nd : NZZD(order=2) */ const int nSteps; const int nOrder; Journal &journal; - Vector *ySteady; - Vector *params; - TwoDMatrix *vCov; + Vector &ySteady; + Vector ¶ms; + TwoDMatrix &vCov; TensorContainer md; // ModelDerivatives - DynareNameList *dnl; - DynareExogNameList *denl; - DynareStateNameList *dsnl; + DynareNameList dnl; + DynareExogNameList denl; + DynareStateNameList dsnl; const double ss_tol; - const vector *varOrder; - const TwoDMatrix *ll_Incidence; + const vector &varOrder; + const TwoDMatrix &ll_Incidence; double qz_criterium; - vector *JacobianIndices; + vector JacobianIndices; public: KordpDynare(const char **endo, int num_endo, const char **exo, int num_exo, int num_par, - Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred, - int nforw, int nboth, const int nJcols, const Vector *NNZD, + 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, DynamicModelDLL &dynamicDLL, double sstol, - const vector *varOrder, const TwoDMatrix *ll_Incidence, + const vector &varOrder, const TwoDMatrix &ll_Incidence, double qz_criterium) throw (TLException); virtual ~KordpDynare(); @@ -203,27 +203,27 @@ public: const NameList & getAllEndoNames() const { - return *dnl; + return dnl; } const NameList & getStateNames() const { - return *dsnl; + return dsnl; } const NameList & getExogNames() const { - return *denl; + return denl; } const TwoDMatrix & getVcov() const { - return *vCov; + return vCov; } Vector & getParams() { - return *params; + return params; } const TensorContainer & @@ -234,31 +234,31 @@ public: const Vector & getSteady() const { - return *ySteady; + return ySteady; } Vector & getSteady() { - return *ySteady; + return ySteady; } void solveDeterministicSteady(); void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) throw (DynareException); void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, const Vector &yyp, const Vector &xx) throw (DynareException); - void calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException); void calcDerivativesAtSteady() throw (DynareException); DynamicModelDLL &dynamicDLL; DynamicModel * clone() const { - return new KordpDynare(*this); + std::cerr << "KordpDynare::clone() not implemented" << std::endl; + exit(EXIT_FAILURE); } - Vector *LLxSteady(const Vector &yS) throw (DynareException, TLException); // returns ySteady extended with leads and lags + void LLxSteady(const Vector &yS, Vector &llxSteady) throw (DynareException, TLException); // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags private: - vector *ReorderDynareJacobianIndices(const vector *varOrder) throw (TLException); - void populateDerivativesContainer(TwoDMatrix *g, int ord, const vector *vOrder); + void ReorderDynareJacobianIndices() throw (TLException); + void populateDerivativesContainer(const TwoDMatrix &g, int ord, const vector &vOrder); }; #endif diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc index 5047dc2a2..93a8e31e0 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cc +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc @@ -130,17 +130,17 @@ extern "C" { mxFldp = mxGetField(M_, 0, "params"); double *dparams = (double *) mxGetData(mxFldp); int npar = (int) mxGetM(mxFldp); - Vector *modParams = new Vector(dparams, npar); + Vector modParams(dparams, npar); mxFldp = mxGetField(M_, 0, "Sigma_e"); dparams = (double *) mxGetData(mxFldp); npar = (int) mxGetN(mxFldp); - TwoDMatrix *vCov = new TwoDMatrix(npar, npar, dparams); + TwoDMatrix vCov(npar, npar, dparams); mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var dparams = (double *) mxGetData(mxFldp); const int nSteady = (int) mxGetM(mxFldp); - Vector *ySteady = new Vector(dparams, nSteady); + Vector ySteady(dparams, nSteady); mxFldp = mxGetField(dr, 0, "nstatic"); const int nStat = (int) mxGetScalar(mxFldp); @@ -172,9 +172,9 @@ extern "C" { npar = (int) mxGetM(mxFldp); if (npar != nEndo) //(nPar != npar) mexErrMsgTxt("Incorrect number of input var_order vars."); - vector *var_order_vp = (new vector(nEndo)); + vector var_order_vp(nEndo); for (int v = 0; v < nEndo; v++) - (*var_order_vp)[v] = (int)(*(dparams++)); + var_order_vp[v] = (int)(*(dparams++)); // the lag, current and lead blocks of the jacobian respectively mxFldp = mxGetField(M_, 0, "lead_lag_incidence"); @@ -182,14 +182,14 @@ extern "C" { npar = (int) mxGetN(mxFldp); int nrows = (int) mxGetM(mxFldp); - TwoDMatrix *llincidence = new TwoDMatrix(nrows, npar, dparams); + TwoDMatrix llincidence(nrows, npar, dparams); if (npar != nEndo) mexErrMsgIdAndTxt("dynare:k_order_perturbation", "Incorrect length of lead lag incidences: ncol=%d != nEndo=%d.", npar, nEndo); //get NNZH =NNZD(2) = the total number of non-zero Hessian elements mxFldp = mxGetField(M_, 0, "NNZDerivatives"); dparams = (double *) mxGetData(mxFldp); - Vector *NNZD = new Vector(dparams, (int) mxGetM(mxFldp)); + Vector NNZD(dparams, (int) mxGetM(mxFldp)); const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns @@ -258,7 +258,7 @@ extern "C" { app.getFoldDecisionRule().writeMMap(mm, string()); // get latest ysteady - ySteady = (Vector *)(&dynare.getSteady()); + ySteady = dynare.getSteady(); if (kOrder == 1) {