From ca0e886204bee4fbc6cc569cb8fc684c311dea75 Mon Sep 17 00:00:00 2001 From: george Date: Mon, 17 Nov 2008 17:00:54 +0000 Subject: [PATCH] Update - not yet working git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2260 ac1d8469-bf42-47a9-8791-bf33cf982152 --- mex/sources/korderpert/src/k_ord_dynare.cpp | 173 +++-- mex/sources/korderpert/src/k_ord_dynare.h | 98 +-- .../korderpert/src/k_order_perturbation.cpp | 682 ++++++++++-------- .../korderpert/src/k_order_perturbation.h | 61 +- 4 files changed, 569 insertions(+), 445 deletions(-) diff --git a/mex/sources/korderpert/src/k_ord_dynare.cpp b/mex/sources/korderpert/src/k_ord_dynare.cpp index 63bd4adfb..66f2bd51a 100644 --- a/mex/sources/korderpert/src/k_ord_dynare.cpp +++ b/mex/sources/korderpert/src/k_ord_dynare.cpp @@ -15,19 +15,11 @@ #include "k_ord_dynare.h" -#include "dynare_exception.h" -#include "decision_rule.h" -#include "fs_tensor.h" -#include "SylvException.h" #include "mex.h" -///#include "forw_subst_builder.h" - #include "Dynare_pp/utils/cc/memory_file.h" #include "Dynare_pp/utils/cc/exception.h" -#include "Dynare_pp/tl/cc/tl_exception.h" -#include "Dynare_pp/kord/kord_exception.h" //#include "k_order_perturbation.h" @@ -36,51 +28,36 @@ #endif -/**************************************************************************************/ -/* DynareNameList class */ -/**************************************************************************************/ -vector DynareNameList::selectIndices(const vector& ns) const -{ - vector res; - for (unsigned int i = 0; i < ns.size(); i++) { - int j = 0; - while (j < getNum() && strcmp(getName(j), ns[i]) != 0) - j++; - if (j == getNum()) - throw DynareException(__FILE__, __LINE__, - string("Couldn't find name for ") + ns[i] + - " in DynareNameList::selectIndices"); - res.push_back(j); - } - return res; -} +FistOrderApproximation::GetRuleDers(double*dgy,double*dgu){ + +}; + ///////////////////////// /**************************************************************************************/ /* Dynare DynamicModel class */ /**************************************************************************************/ - +class KordpJacobian; KordpDynare::KordpDynare(const char** endo, int num_endo, - const char** exo, int nexog, - const char** par, int nPar, - Vector* ysteady, int nstat,int nPred, int nForw, int nboth, - const char* modName, int len, int nOrder, - double sstol, Journal& jr, DynamicModelDLL& dynamicDLL) - : journal(jr), md(1), ysteady(NULL), - dnl(NULL), denl(NULL), dsnl(NULL), nPar(nPar), nOrder(nOrder), - ss_tol(sstol), nStat(nstat), nBoth(nboth), - nForw(nForw), nExog(nexog), nPred(nPred), - nYs(nYs), nYss(nYss),nY(nY), dynamicDLL(dynamicDLL) - + const char** exo, int nexog, int nPar, //const char** par, + Vector* ySteady, TwoDMatrix* vCov, Vector* params, int nstat,int nPred, int nForw, int nboth, + const int nSteps, int nOrder, //const char* modName, + Journal& jr, DynamicModelDLL& dynamicDLL, double sstol) + : nStat(nstat), nBoth(nboth), nPred(nPred), nForw(nForw), nExog(nexog), nPar(nPar), + nYs(nYs), nYss(nYss),nY(nY), nSteps(nSteps), nOrder(nOrder), journal(jr), dynamicDLL(dynamicDLL), + ySteady(ySteady), vCov(vCov), params (params), md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol) { - try{ - - dnl = new DynareNameList(*this); - denl = new DynareExogNameList(*this); + dnl = new DynareNameList(*this, endo); + denl = new DynareExogNameList(*this, exo); dsnl = new DynareStateNameList(*this, *dnl, *denl); - throw DynareException(__FILE__, __LINE__, string("Could not open model file ")+modName); +/**** + ySteady = new Vector(*(ySteady)); + params = new Vector(*(params)); + vCov = new TwoDMatrix(*(vCov)); +******/ + // throw DynareException(__FILE__, __LINE__, string("Could not open model file ")+modName); } catch (...) {} @@ -89,50 +66,41 @@ KordpDynare::KordpDynare(const char** endo, int num_endo, /// writeModelInfo(journal); } -class KordpJacobian; - KordpDynare::KordpDynare(const KordpDynare& dynare) - : journal(dynare.journal),///, model(NULL), - ysteady(NULL), md(dynare.md), - dnl(NULL), denl(NULL), dsnl(NULL), - ss_tol(dynare.ss_tol), nStat(dynare.nStat), nBoth(dynare.nBoth), - nPred(dynare.nPred), nForw(dynare.nForw), nExog(dynare.nExog), - nYs(dynare.nYs), nYss(dynare.nYss),nY(dynare.nY), - dynamicDLL(dynare.dynamicDLL),nOrder(dynare.nOrder), nPar(dynare.nPar) + : nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred), + nForw(dynare.nForw), nExog(dynare.nExog), nPar(dynare.nPar), + nYs(dynare.nYs), nYss(dynare.nYss),nY(dynare.nY), + nSteps(dynare.nSteps), nOrder(dynare.nOrder), journal(dynare.journal), + dynamicDLL(dynare.dynamicDLL), //modName(dynare.modName), + ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md), + dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol) { /// model = dynare.model->clone(); - ysteady = new Vector(*(dynare.ysteady)); + ySteady = new Vector(*(dynare.ySteady)); params = new Vector(*(dynare.params)); - Vcov = new TwoDMatrix(*(dynare.Vcov)); + vCov = new TwoDMatrix(*(dynare.vCov)); // if (dynare.md)// !=NULL) // Inititalise ModelDerivatives md - md= *(new TensorContainer::TensorContainer(dynare.md)); - dnl = new DynareNameList(*this); - denl = new DynareExogNameList(*this); +// md= *(new TensorContainer(dynare.md)); + dnl = new DynareNameList(dynare);//(*this); + denl = new DynareExogNameList(dynare);//(*this); dsnl = new DynareStateNameList(*this, *dnl, *denl); } KordpDynare::~KordpDynare() { - if (ysteady) - delete ysteady; + if (ySteady) + delete ySteady; if (params) delete params; - if (Vcov) - delete Vcov; -/*************** May be needed + if (vCov) + delete vCov; if (dnl) delete dnl; if (dsnl) delete dsnl; if (denl) delete denl; - if (fe) - delete fe; - if (fde) - delete fde; -***********/// - } void KordpDynare::solveDeterministicSteady(Vector& steady) @@ -144,7 +112,7 @@ void KordpDynare::solveDeterministicSteady(Vector& steady) KordpJacobian dj(*this); ogu::NLSolver nls(dvf, dj, 500, ss_tol, journal); int iter; - if (! nls.solve(*ysteady, iter)) + if (! nls.solve(*ySteady, iter)) throw DynareException(__FILE__, __LINE__, "Could not obtain convergence in non-linear solver"); } @@ -193,9 +161,7 @@ void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx) params, //int it_, out, g1, NULL); -// if (!md){ -// md=new TensorContainer (1); - // model derivatives FSSparseTensor instance + // model derivatives FSSparseTensor instance FSSparseTensor mdTi=*(new FSSparseTensor (1, g1->ncols(),g1->nrows())); for (int i = 0; incols(); i++){ for (int j = 0; jnrows(); j++){ @@ -204,17 +170,15 @@ void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx) } } // md container - md=*(new TensorContainer::TensorContainer(1)); - +// md=*(new TensorContainer(1)); + md.clear(); md.insert(&mdTi); -// }else{} - } void KordpDynare::calcDerivativesAtSteady() { Vector xx(nexog()); xx.zeros(); - calcDerivatives(*ysteady, xx); + calcDerivatives(*ySteady, xx); } void KordpDynare::writeModelInfo(Journal& jr) const @@ -288,3 +252,58 @@ void KordpVectorFunction::eval(const ConstVector& in, Vector& out) d.evaluateSystem(out, in, xx); } +/**************************************************************************************/ +/* DynareNameList class */ +/**************************************************************************************/ +vector DynareNameList::selectIndices(const vector& ns) const +{ + vector res; + for (unsigned int i = 0; i < ns.size(); i++) { + int j = 0; + while (j < getNum() && strcmp(getName(j), ns[i]) != 0) + j++; + if (j == getNum()) + throw DynareException(__FILE__, __LINE__, + string("Couldn't find name for ") + ns[i] + + " in DynareNameList::selectIndices"); + res.push_back(j); + } + return res; +} + +DynareNameList::DynareNameList(const KordpDynare& dynare) +{ + for (int i = 0; i < dynare.ny(); i++) { + names.push_back(dynare.dnl->getName(i)); + } +} +DynareNameList::DynareNameList(const KordpDynare& dynare, const char ** namesp) +{ + for (int i = 0; i < dynare.ny(); i++) { + names.push_back(namesp[i]); + } +} + +DynareExogNameList::DynareExogNameList(const KordpDynare& dynare) +{ + for (int i = 0; i < dynare.nexog(); i++) { + names.push_back(dynare.denl->getName(i)); + } +} + +DynareExogNameList::DynareExogNameList(const KordpDynare& dynare, const char ** namesp) +{ + for (int i = 0; i < dynare.nexog(); i++) { + names.push_back(namesp[i]); + } +} + +DynareStateNameList::DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl, + const DynareExogNameList& denl) +{ + for (int i = 0; i < dynare.nys(); i++) + names.push_back(dnl.getName(i+dynare.nstat())); + for (int i = 0; i < dynare.nexog(); i++) + names.push_back(denl.getName(i)); +} + diff --git a/mex/sources/korderpert/src/k_ord_dynare.h b/mex/sources/korderpert/src/k_ord_dynare.h index 8a1f8c584..861f54b44 100644 --- a/mex/sources/korderpert/src/k_ord_dynare.h +++ b/mex/sources/korderpert/src/k_ord_dynare.h @@ -18,24 +18,40 @@ #ifndef K_ORD_DYNARE3_H #define K_ORD_DYNARE3_H #include -#include "Dynare_pp/tl/cc/t_container.h" -#include "Dynare_pp/tl/cc/sparse_tensor.h" -#include "Dynare_pp/kord/decision_rule.h" -#include "Dynare_pp/kord/dynamic_model.h" +#include "t_container.h" +#include "sparse_tensor.h" +#include "decision_rule.h" +#include "dynamic_model.h" +#include "dynare_exception.h" +#include "fs_tensor.h" +#include "SylvException.h" +#include "tl_exception.h" +#include "kord_exception.h" +#include "exception.h" #include "nlsolve.h" #include "k_order_perturbation.h" class KordpDynare; +// derive from Approximation to get protected derivatives +class FistOrderApproximation: public Approximation{ + public: + FGSContainer* GetRuleDers(){return rule_ders;}; + void GetRuleDers(double*dgy,double*dgu); + FGSContainer* GetRuleDersSS(){return rule_ders_ss;}; + void GetRuleDersSS(double*dgy, double*dgu); +} + // instantiations of pure abstract class NameList in dynamic_model.h: /*////////////////////////////////////////////*/ class DynareNameList : public NameList { -vector names; + vector names; public: DynareNameList(const KordpDynare& dynare); + DynareNameList(const KordpDynare& dynare, const char ** names); int getNum() const {return (int)names.size();} const char* getName(int i) const {return names[i];} /** This for each string of the input vector calculates its index @@ -44,26 +60,27 @@ public: vector selectIndices(const vector& ns) const; }; - class DynareExogNameList : public NameList { - vector names; - public: - DynareExogNameList(const KordpDynare& dynare); - int getNum() const - {return (int)names.size();} - const char* getName(int i) const - {return names[i];} - }; +class DynareExogNameList : public NameList { + vector names; +public: + DynareExogNameList(const KordpDynare& dynare); + DynareExogNameList(const KordpDynare& dynare, const char ** names); + int getNum() const + {return (int)names.size();} + const char* getName(int i) const + {return names[i];} +}; - class DynareStateNameList : public NameList { - vector names; - public: - DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl, - const DynareExogNameList& denl); - int getNum() const - {return (int)names.size();} - const char* getName(int i) const - {return names[i];} - }; +class DynareStateNameList : public NameList { + vector names; +public: + DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl, + const DynareExogNameList& denl); + int getNum() const + {return (int)names.size();} + const char* getName(int i) const + {return names[i];} +}; /*********************************************/ // The following only implements DynamicModel with help of ogdyn::DynareModel // instantiation of pure abstract DynamicModel decl. in dynamic_model.h @@ -74,6 +91,7 @@ friend class DynareNameList; friend class DynareExogNameList; friend class DynareStateNameList; friend class KordpDynareJacobian; +friend class DynamicModelDLL; ////////// const int nStat; const int nBoth; @@ -84,31 +102,25 @@ friend class KordpDynareJacobian; const int nYs; const int nYss; const int nY; + const int nSteps; const int nOrder; Journal& journal; - /// ogdyn::DynareModel* model; /// DynamicModel* model; - ///DynamicModelDLL& dynamicDLL; - const char* modName; - Vector* ysteady; + ///const char* modName; + Vector* ySteady; Vector* params; + TwoDMatrix* vCov; TensorContainer md; // ModelDerivatives DynareNameList* dnl; DynareExogNameList* denl; DynareStateNameList* dsnl; - - TwoDMatrix* Vcov; - - /// ogp::FormulaEvaluator* fe; - /// ogp::FormulaDerEvaluator* fde; const double ss_tol; public: KordpDynare(const char** endo, int num_endo, - const char** exo, int num_exo, - const char** par, int num_par, - Vector* ysteady, int nstat,int nPred, int nforw, int nboth, - const char* modName, int len, int ord, - double sstol, Journal& jr, DynamicModelDLL& dynamicDLL); + const char** exo, int num_exo, int num_par, //const char** par, + Vector* ySteady, TwoDMatrix* vCov, Vector* params, int nstat,int nPred, + int nforw, int nboth, const int nSteps, const int ord, //const char* modName, + Journal& jr, DynamicModelDLL& dynamicDLL, double sstol); /** Makes a deep copy of the object. */ KordpDynare(const KordpDynare& dyn); @@ -129,6 +141,8 @@ public: {return nYss;} int ny() const {return nY;} + int steps() const + {return nSteps;} int order() const {return nOrder;} const NameList& getAllEndoNames() const @@ -138,23 +152,23 @@ public: const NameList& getExogNames() const {return *denl;} TwoDMatrix& getVcov() const - {return *Vcov;} + {return *vCov;} Vector& getParams() {return *params;} const TensorContainer& getModelDerivatives() const {return md;} const Vector& getSteady() const - {return *ysteady;} + {return *ySteady;} Vector& getSteady() - {return *ysteady;} + {return *ySteady;} /// const ogdyn::DynareModel& getModel() const /// {return *model;} // here is true public interface void solveDeterministicSteady(Vector& steady); void solveDeterministicSteady() - {solveDeterministicSteady(*ysteady);} + {solveDeterministicSteady(*ySteady);} void evaluateSystem(Vector& out, const Vector& yy, const Vector& xx); void evaluateSystem(Vector& out, const Vector& yym, const Vector& yy, const Vector& yyp, const Vector& xx); diff --git a/mex/sources/korderpert/src/k_order_perturbation.cpp b/mex/sources/korderpert/src/k_order_perturbation.cpp index 8365f1596..c6bd86f99 100644 --- a/mex/sources/korderpert/src/k_order_perturbation.cpp +++ b/mex/sources/korderpert/src/k_order_perturbation.cpp @@ -1,42 +1,48 @@ /* 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 . - */ - +* 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 . +*/ + // k_order_perturbation.cpp : Defines the entry point for the DLL application. // +// Receives call from Dynare resol and runs instead of dr1.m +// [dr,info,M_,options_,oo_] = k_ord_perturbation(M_,options_,oo_); +// receives M_,options_,oo_ and returns dr1, in addition to M_,options_,oo_ +// -#include "stdafx.h" #include "k_ord_dynare.h" #include "math.h" -#include "mex.h" -#include "k_order_perturbation.h" +#include "string.h" +//#include "mex.h" +//#include "k_order_perturbation.h" #include "approximation.h" -/* +#ifdef _MSC_VER //&&WINDOWS + + BOOL APIENTRY DllMain( HANDLE hModule, - DWORD ul_reason_for_call, - LPVOID lpReserved - ) + DWORD ul_reason_for_call, + LPVOID lpReserved + ) { switch (ul_reason_for_call) { - case DLL_PROCESS_ATTACH: - case DLL_THREAD_ATTACH: - case DLL_THREAD_DETACH: - case DLL_PROCESS_DETACH: - break; + case DLL_PROCESS_ATTACH: + case DLL_THREAD_ATTACH: + case DLL_THREAD_DETACH: + case DLL_PROCESS_DETACH: + break; } return TRUE; } @@ -58,299 +64,385 @@ CK_order_perturbation::CK_order_perturbation() return; } -*/ +#endif // _MSC_VER && WINDOWS extern "C" { -void mexFunction(int nlhs, mxArray* plhs[], - int nrhs, const mxArray* prhs[]) -{ - if (nrhs < 3) - mexErrMsgTxt("Must have at least 3 input parameters.\n"); - if (nlhs == 0) - mexErrMsgTxt("Must have at least 1 output parameter.\n"); - -/* int order = (int)mxGetScalar(prhs[0]); - if (nrhs != 12 + order) { - mexErrMsgTxt("Must have exactly 11+order input parameters.\n"); -// return; - } -*/ -// const mxArray* const aa = prhs[0]; -/* const mxArray* rhs1in = (prhs[0]); - mxArray* rhs1[1];// = rhs1in; - rhs1[0] = const_cast (rhs1in); -*/ -// const mxArray* const mFname = prhs[1]; -// const char* mFname = mxArrayToString(prhs[nrhs-1]); - - char* mFname = mxArrayToString(prhs[nrhs-1]); - - - // SIM: - int nstat = (int)mxGetScalar(prhs[1]); - int npred = (int)mxGetScalar(prhs[2]); - int nboth = (int)mxGetScalar(prhs[3]); - int nforw = (int)mxGetScalar(prhs[4]); - int nexog = (int)mxGetScalar(prhs[5]); - - const mxArray* const ystart = prhs[6]; - const mxArray* const shocks = prhs[7]; - const mxArray* const vcov = prhs[8]; - int seed = (int)mxGetScalar(prhs[9]); - const mxArray* const ysteady = prhs[10]; - const int* const ystart_dim = mxGetDimensions(ystart); - const int* const shocks_dim = mxGetDimensions(shocks); - const int* const vcov_dim = mxGetDimensions(vcov); - const int* const ysteady_dim = mxGetDimensions(ysteady); - - - - int ny = nstat + npred + nboth + nforw; - if (ny != ystart_dim[0]) - mexErrMsgTxt("ystart has wrong number of rows.\n"); - if (1 != ystart_dim[1]) - mexErrMsgTxt("ystart has wrong number of cols.\n"); - int nper = shocks_dim[1]; - if (nexog != shocks_dim[0]) - mexErrMsgTxt("shocks has a wrong number of rows.\n"); - if (nexog != vcov_dim[0]) - mexErrMsgTxt("vcov has a wrong number of rows.\n"); - if (nexog != vcov_dim[1]) - mexErrMsgTxt("vcov has a wrong number of cols.\n"); - if (ny != ysteady_dim[0]) - mexErrMsgTxt("ysteady has wrong number of rows.\n"); - if (1 != ysteady_dim[1]) - mexErrMsgTxt("ysteady has wrong number of cols.\n"); - - mxArray* res = mxCreateDoubleMatrix(ny, nper, mxREAL); - THREAD_GROUP::max_parallel_threads = params.num_threads; - - try { - // make journal name and journal - std::string jname(params.basename); - jname += ".jnl"; - Journal journal(jname.c_str()); - - // make dynare object - Dynare dynare(params.modname, params.order, params.ss_tol, journal); - // make list of shocks for which we will do IRFs - vector irf_list_ind; - if (params.do_irfs_all) - for (int i = 0; i < dynare.nexog(); i++) - irf_list_ind.push_back(i); - else - irf_list_ind = ((const DynareNameList&)dynare.getExogNames()).selectIndices(params.irf_list); - - - try { - - system_random_generator.initSeed(params.seed); - - tls.init(dynare.order(), - dynare.nstat()+2*dynare.npred()+3*dynare.nboth()+ - 2*dynare.nforw()+dynare.nexog()); - - Approximation app(dynare, journal, params.num_steps); - - } catch (const KordException& e) { - // tell about the exception and continue - printf("Caught (not yet fatal) Kord exception: "); - e.print(); - JournalRecord rec(journal); - rec << "Solution routine not finished (" << e.get_message() - << "), see what happens" << endrec; - } catch (const TLException& e) { - mexErrMsgTxt("Caugth TL exception."); - } catch (SylvException& e) { - mexErrMsgTxt("Caught Sylv exception."); - } - - - try { - - app.walkStochSteady(); + + // Receives call from Dynare resol and runs instead of dr1.m + // but you need to run set_state_space first which s not in resol!!! + // [dr,info,M_,options_,oo_] = + // k_ord_perturbation(dr,check_flag,M_,options_,oo_); + // receives dr, check_flag,, M_,options_,oo_ and + // returns dr1, in addition to M_,options_,oo_ + + void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) + { + if (nrhs < 5) + mexErrMsgTxt("Must have at least 5 input parameters.\n"); + if (nlhs == 0) + mexErrMsgTxt("Must have at least 1 output parameter.\n"); - } catch (const KordException& e) { + const mxArray* dr = prhs[0]; + const int check_flag = (int)mxGetScalar(prhs[1]); + const mxArray* M_ = prhs[2]; + const mxArray* options_= (prhs[3]); + const mxArray* oo_ = (prhs[4]); + + mxArray* mFname =mxGetField(M_, 0, "fname"); + if(!mxIsChar(mFname)){ + mexErrMsgTxt("Input must be of type char."); + } + const char* fName = mxArrayToString(mFname); + +#ifdef DEBUG + mexPrintf("MexPrintf 2: check_flag = %d , fName = %s .\n", check_flag,fName); +#endif + int kOrder; + mxArray* mxFldp = mxGetField(options_, 0,"order" ); + if (mxIsNumeric(mxFldp)) + kOrder = (int)mxGetScalar(mxFldp); + else + kOrder = 1; + if (nrhs != 12 + kOrder) { + mexErrMsgTxt("Must have exactly 11+order input parameters.\n"); + return; + } + mxFldp = mxGetField(M_, 0,"params" ); + double * dparams = (double *) mxGetData(mxFldp); + int npar = (int)mxGetN(mxFldp); + Vector * modParams = new Vector(dparams, npar); + + mxFldp = mxGetField(M_, 0,"Sigma_e" ); + dparams = (double *) mxGetData(mxFldp); + npar = (int)mxGetN(mxFldp); + TwoDMatrix * vCov = new TwoDMatrix(npar, npar, dparams); + + + mxFldp = mxGetField(oo_, 0,"steady_state" ); // use in order of declaration +// mxFldp = mxGetField(dr, 0,"ys" ); // and not in order of dr.order_var + dparams = (double *) mxGetData(mxFldp); + npar = (int)mxGetN(mxFldp); + Vector * ySteady = new Vector(dparams, npar); + + + mxFldp = mxGetField(dr, 0,"nstatic" ); + const int nStat = (int)mxGetScalar(mxFldp); + mxFldp = mxGetField(dr, 0,"npred" ); + const int nPred = (int)mxGetScalar(mxFldp); + mxFldp = mxGetField(dr, 0,"nboth" ); + const int nBoth = (int)mxGetScalar(mxFldp); + mxFldp = mxGetField(dr, 0,"nfwrd" ); + const int nForw = (int)mxGetScalar(mxFldp); + + mxFldp = mxGetField(M_, 0,"exo_nbr" ); + const int nExog = (int)mxGetScalar(mxFldp); + mxFldp = mxGetField(M_, 0,"endo_nbr" ); + const int nEndo = (int)mxGetScalar(mxFldp); + mxFldp = mxGetField(M_, 0,"param_nbr" ); + const int nPar = (int)mxGetScalar(mxFldp); + + + mxFldp = mxGetField(M_, 0,"endo_names" ); + const char ** endoNames = (const char **) mxGetData(mxFldp); + const int nendo = (int)mxGetN(mxFldp); + + mxFldp = mxGetField(M_, 0,"exo_names" ); + const char ** exoNames = (const char **) mxGetData(mxFldp); + const int nexo = (int)mxGetN(mxFldp); +/****** + mxFldp = mxGetField(M_, 0,"param_names" ); + const char ** paramNames = (char **) mxGetData(mxFldp); + const int npar = (int)mxGetN(mxFldp); +************/ + if ((nEndo != nendo)||(nExog != nexo)) { //(nPar != npar) + mexErrMsgTxt("Incorrect number of input parameters.\n"); + //return; + } + +#ifdef DEBUG + mexPrintf("MexPrintf 2: check_flag = %d , fName = %s .\n", check_flag,fName); + mexPrintf("MexPrintf 2: nEndo = %d , nExo = %d .\n", nEndo,nExog); +#endif + /* Fetch time index */ + // int it_ = (int) mxGetScalar(prhs[3]) - 1; + + 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 + + THREAD_GROUP::max_parallel_threads = 2;//params.num_threads; + + try { + // make journal name and journal + std::string jName(fName); //params.basename); + jName += ".jnl"; + Journal journal(jName.c_str()); + + // DynamicFn * pDynamicFn = loadModelDynamicDLL (fname); + DynamicModelDLL dynamicDLL(fName); + + // make KordpDynare object + KordpDynare dynare(endoNames, nEndo, exoNames, nExog, nPar, // paramNames, + ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, + nSteps, kOrder, journal, dynamicDLL, sstol); + /************ + // make list of shocks for which we will do IRFs + vector irf_list_ind; + if (params.do_irfs_all){ + for (int i = 0; i < dynare.nexog(); i++) + irf_list_ind.push_back(i); + } + else + irf_list_ind = ((const DynareNameList&)dynare.getExogNames()).selectIndices(params.irf_list); + ****************/ + try { + // intiate tensor library + tls.init(dynare.order(), + dynare.nstat()+2*dynare.npred()+3*dynare.nboth()+ + 2*dynare.nforw()+dynare.nexog()); + + // construct main K-order approximation class + FistOrderApproximation app(dynare, journal, nSteps); + // run stochastic steady + app.walkStochSteady(); + + } catch (const KordException& e) { // tell about the exception and continue printf("Caught (not yet fatal) Kord exception: "); e.print(); JournalRecord rec(journal); rec << "Solution routine not finished (" << e.get_message() << "), see what happens" << endrec; - } catch (const TLException& e) { - mexErrMsgTxt("Caugth TL exception."); - } catch (SylvException& e) { - mexErrMsgTxt("Caught Sylv exception."); - } + } catch (const TLException& e) { + mexErrMsgTxt("Caugth TL exception."); + } catch (SylvException& e) { + mexErrMsgTxt("Caught Sylv exception."); + } - std::string ss_matrix_name(params.prefix); - ss_matrix_name += "_steady_states"; - ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str()); - - // check the approximation - if (params.check_along_path || params.check_along_shocks + /*********************** + std::string ss_matrix_name("K_ordp");//params.prefix); + ss_matrix_name += "_steady_states"; + + // ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str()); + + May be needed to + // check the approximation + if (params.check_along_path || params.check_along_shocks || params.check_on_ellipse) { GlobalChecker gcheck(app, THREAD_GROUP::max_parallel_threads, journal); if (params.check_along_shocks) - gcheck.checkAlongShocksAndSave(matfd, params.prefix, - params.getCheckShockPoints(), - params.getCheckShockScale(), - params.check_evals); + gcheck.checkAlongShocksAndSave(matfd, params.prefix, + params.getCheckShockPoints(), + params.getCheckShockScale(), + params.check_evals); if (params.check_on_ellipse) - gcheck.checkOnEllipseAndSave(matfd, params.prefix, + gcheck.checkOnEllipseAndSave(matfd, params.prefix, params.getCheckEllipsePoints(), params.getCheckEllipseScale(), params.check_evals); - if (params.check_along_path) - gcheck.checkAlongSimulationAndSave(matfd, params.prefix, - params.getCheckPathPoints(), - params.check_evals); - } - } catch (const KordException& e) { - printf("Caugth Kord exception: "); - e.print(); - return e.code(); - } catch (const TLException& e) { - printf("Caugth TL exception: "); - e.print(); - return 255; - } catch (SylvException& e) { - printf("Caught Sylv exception: "); - e.printMessage(); - return 255; - } catch (const DynareException& e) { - printf("Caught Dynare exception: %s\n", e.message()); - return 255; - } catch (const ogu::Exception& e) { - printf("Caught ogu::Exception: "); - e.print(); - return 255; + if (params.check_along_path) + gcheck.checkAlongSimulationAndSave(matfd, params.prefix, + params.getCheckPathPoints(), + params.check_evals); + } + *****************************/ + + // get protected derivatives from Approximation + FGSContainer* rule_ders = app.GetRuleDers(); + FGSContainer* rule_ders_ss = app.GetRuleDersSS(); + + // get latest ysteady + ysteady = dynare.getSteady()->base(); + + + + + + } catch (const KordException& e) { + printf("Caugth Kord exception: "); + e.print(); + return;// e.code(); + } catch (const TLException& e) { + printf("Caugth TL exception: "); + e.print(); + return;// 255; + } catch (SylvException& e) { + printf("Caught Sylv exception: "); + e.printMessage(); + return;// 255; + } catch (const DynareException& e) { + printf("Caught KordpDynare exception: %s\n", e.message()); + return;// 255; + } catch (const ogu::Exception& e) { + printf("Caught ogu::Exception: "); + e.print(); + return;// 255; + } + + double *gy, *gu; + int nb_row_x; + + ysteady = NULL; + if (nlhs >= 1) + { + /* Set the output pointer to the output matrix residual. */ + plhs[0] = mxCreateDoubleMatrix(nEndo,1, mxREAL); + /* Create a C pointer to a copy of the output matrix residual. */ + ysteady = mxGetPr(plhs[0]); + } + + gy = NULL; + if (nlhs >= 2) + { + /* Set the output pointer to the output matrix g1. */ + plhs[1] = mxCreateDoubleMatrix(nEndo, nEndo, mxREAL); + /* Create a C pointer to a copy of the output matrix g1. */ + gy = mxGetPr(plhs[1]); + } + + gu = NULL; + if (nlhs >= 3) + { + /* Set the output pointer to the output matrix g2. */ + plhs[2] = mxCreateDoubleMatrix(nExog, nExog, mxREAL); + /* Create a C pointer to a copy of the output matrix g1. */ + gu = mxGetPr(plhs[2]); + } + + + + /* Call the C subroutines. */ +// delete params; +// delete vCov; +// delete ySteady; + }; + +// Load _Dyamic DLL ); +DynamicModelDLL::DynamicModelDLL(const char * modName) +{ + char fName[MAX_MODEL_NAME]; + strcpy(fName,modName); + using namespace std; + // string sFname(mFname); + // string sFname(fname); + // string sExt("_.dll"); +// mexPrintf("MexPrintf: Call exp %d.\n", y[0]); + mexPrintf("MexPrintf: Call Load run DLL %s .\n", fName); + + try { + // typedef void * (__stdcall *DynamicFn)(); + +#ifdef WINDOWS + + HINSTANCE dynamicHinstance; +// dynamicHinstance=::LoadLibraryEx(strcat(fNname,"_.dll"),NULL,DONT_RESOLVE_DLL_REFERENCES);//sExt); //"_.dll"); + dynamicHinstance=::LoadLibrary(strcat(fName,"dynamic.dll"));//sExt); //"_.dll"); + if (dynamicHinstance==NULL) + throw 1; //alt: return; + // (DynamicFn*) typedef void * (__stdcall *DynamicFn)(); + mexPrintf("MexPrintf: Call GetProcAddress %s .\n", fName); + Dynamic = (DynamicFn *) ::GetProcAddress(dynamicHinstance,"Dynamic"); + +# else // linux + + void *dynamicHinstance = dlopen(strcat(fNname,"_dynamic.so"), RTLD_NOW); + if((dynamicHinstance == NULL) || dlerr()){ + cerr << dlerror() << endl; + mexPrintf("MexPrintf:Error loading DLL: %s", dlerror); + throw 1; + } + void *mkr = dlsym(dynamicHinstance, "Dynamic"); + if((mkr == NULL) || dlerr()){ + cerr << dlerror() << endl; + mexPrintf("MexPrintf:Error finding DLL function: %s", dlerror); + throw 2; + } + //The pointer to maker must be of type void *, since that is the type returned + //DynamicFn * + Dynamic = static_cast(mkr)(); +# endif + +// if (Dynamic == NULL) +// throw 3; //return; + + } catch (int i) { + mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i); + mexErrMsgTxt("Err: An error in Load and run DLL .\n"); + return; + + } catch (...) { + mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName); + mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n"); + return; + } } -void loadModelDynamicDLL () -{ - double *y, *x, *params; - double *residual, *g1, *g2; - int nb_row_x, it_; + /* + try { + + // Call int mexCallMATLAB(int nlhs, mxArray *plhs[], int nrhs, + // mxArray *prhs[], const char *name); + int success = mexCallMATLAB( nlhs, plhs, nrhs-1, rhs1 , mFname); + // int success = mexCallMATLAB( nlhs, plhs, nrhs-1, (struct mxArray *[]) &aa , mFname); + //plhs[0] = res; + + + + + } catch (...) { + mexErrMsgTxt("Err: Unknown error in Call MATLAB function .\n"); + mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", mFname); + return; + } + */ + // basic_string sFname(mFname); + + - /* Create a pointer to the input matrix y. */ - y = mxGetPr(prhs[0]); - - /* Create a pointer to the input matrix x. */ - x = mxGetPr(prhs[1]); - - /* Create a pointer to the input matrix params. */ - params = mxGetPr(prhs[2]); - - /* Fetch time index */ - it_ = (int) mxGetScalar(prhs[3]) - 1; - - /* Gets number of rows of matrix x. */ - nb_row_x = mxGetM(prhs[1]); - - residual = NULL; - if (nlhs >= 1) - { - /* Set the output pointer to the output matrix residual. */ - plhs[0] = mxCreateDoubleMatrix(16,1, mxREAL); - /* Create a C pointer to a copy of the output matrix residual. */ - residual = mxGetPr(plhs[0]); - } - - g1 = NULL; - if (nlhs >= 2) - { - /* Set the output pointer to the output matrix g1. */ - plhs[1] = mxCreateDoubleMatrix(16, 31, mxREAL); - /* Create a C pointer to a copy of the output matrix g1. */ - g1 = mxGetPr(plhs[1]); - } - - g2 = NULL; - if (nlhs >= 3) - { - /* Set the output pointer to the output matrix g2. */ - plhs[2] = mxCreateDoubleMatrix(16, 961, mxREAL); - /* Create a C pointer to a copy of the output matrix g1. */ - g2 = mxGetPr(plhs[2]); - } - - /* Call the C subroutines. */ +// close DLL: If the referenced object was successfully closed, +// close() returns 0, non 0 otherwise +int DynamicModelDLL::close(){ +#ifdef WINDOWS + // MS FreeLibrary returns non 0 if OK, 0 if fails. + bool rb=FreeLibrary(dynamicHinstance); + if (rb) + return 0; + else + return 1; +# else // linux + //If OK, dlclose() returns 0, non 0 otherwise + return dlclose(dynamicHinstance); +# endif +} +void DynamicModelDLL::eval(Vector&y, TwoDMatrix&x, Vector&modParams, + int it_, Vector&residual, TwoDMatrix&g1, TwoDMatrix&g2){ -/* - try { +// DynamicDLLfunc(double *y, double *x, int nb_row_x, double *params, +// int it_, double *residual, double *g1, double *g2) + double *dy, *dx, *dresidual, *dg1, *dg2, dbParams; + + dy=y.base(); + dx=x[it_]; + dbParams=modParams.base(); + Dynamic(dy, dx, int nb_row_x, dbParams, 1, dresidual, dg1, dg2) - // Call int mexCallMATLAB(int nlhs, mxArray *plhs[], int nrhs, - // mxArray *prhs[], const char *name); - int success = mexCallMATLAB( nlhs, plhs, nrhs-1, rhs1 , mFname); -// int success = mexCallMATLAB( nlhs, plhs, nrhs-1, (struct mxArray *[]) &aa , mFname); - //plhs[0] = res; - - - - - } catch (...) { - mexErrMsgTxt("Err: Unknown error in Call MATLAB function .\n"); - mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", mFname); - return; - } -*/ -// basic_string sFname(mFname); - using namespace std; - string sFname(mFname); - string sExt("_.dll"); - mexPrintf("MexPrintf: Call exp %d.\n", y[0]); - double dd = exp(y[0]); - mexPrintf("MexPrintf: exp (%d)= %d .\n", y[0],dd); - - try { -// typedef void * (__stdcall *DynamicFn)(); - typedef void * (DynamicFn) - (double *y, double *x, int nb_row_x, double *params, - int it_, double *residual, double *g1, double *g2); - HINSTANCE DynamicHinstance; - //DynamicFn * pDynamicFn; - mexPrintf("MexPrintf: Call Load run DLL %s .\n", mFname); -// DynamicHinstance=::LoadLibraryEx(strcat(mFname,"_.dll"),NULL,DONT_RESOLVE_DLL_REFERENCES);//sExt); //"_.dll"); - DynamicHinstance=::LoadLibrary(strcat(mFname,"_.dll"));//sExt); //"_.dll"); - if (DynamicHinstance==NULL) - //return; - throw 1; -// typedef void * (__stdcall *DynamicFn)(); - mexPrintf("MexPrintf: 2nd Call exp %d.\n", y[1]); - double dd = exp(y[1]); - mexPrintf("MexPrintf: exp (%d)= %d .\n", y[1],dd); - - mexPrintf("MexPrintf: Call GetProcAddress %s .\n", mFname); - DynamicFn * pDynamicFn = - (DynamicFn*) ::GetProcAddress(DynamicHinstance,"Dynamic"); - if (pDynamicFn == NULL) - //return; - throw 2; - mexPrintf("MexPrintf: Call Dynamic %s .\n", mFname); -// void * objptr = (*pDynamicFn)(y, x, nb_row_x, params, it_, residual, g1, g2); - try{ - (*pDynamicFn)(y, x, nb_row_x, params, it_, residual, g1, g2); -// if (objptr == NULL) - //return; - }catch (...){ - DWORD dw = GetLastError(); - mexPrintf("MexPrintf: error in Call Dynamic DLL %s, %d\n", mFname, dw); - throw 3; - } - - } catch (int i) { - mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", mFname, i); - mexErrMsgTxt("Err: An error in Load and run DLL .\n"); - return; - - } catch (...) { - mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", mFname); - mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n"); - return; - } +// g1= &(TwoDMatrix(double* d, int rows, int cols)); + g1= (TwoDMatrix(dg1, rows, cols)); + residual=Vector(dresidual); + } +void DynamicModelDLL::eval(Vector&y, Vector&x, Vector&modParams, + Vector&residual, TwoDMatrix&g1, TwoDMatrix&g2){ + + dy=y.base(); + dx=x.base(); + dbParams=modParams.base(); + Dynamic(dy, dx, int nb_row_x, dbParams, 1, dresidual, dg1, dg2) + +// g1= &(TwoDMatrix(double* d, int rows, int cols)); + g1= &(TwoDMatrix(dg1, rows, cols)); + + residual=Vector(dresidual); } -} \ No newline at end of file diff --git a/mex/sources/korderpert/src/k_order_perturbation.h b/mex/sources/korderpert/src/k_order_perturbation.h index d5413e97a..91aa1aafc 100644 --- a/mex/sources/korderpert/src/k_order_perturbation.h +++ b/mex/sources/korderpert/src/k_order_perturbation.h @@ -1,16 +1,16 @@ /* 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 . - */ +* 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 . +*/ // The following ifdef block is the standard way of creating macros which make exporting // from a DLL simpler. All files within this DLL are compiled with the K_ORDER_PERTURBATION_EXPORTS @@ -19,15 +19,15 @@ // K_ORDER_PERTURBATION_API functions as being imported from a DLL, wheras this DLL sees symbols // defined with this macro as being exported. - +#ifdef WINDOWS #ifdef K_ORDER_PERTURBATION_EXPORTS #define K_ORDER_PERTURBATION_API __declspec(dllexport) #else #define K_ORDER_PERTURBATION_API __declspec(dllimport) #endif -#ifdef WINDOWS #include "stdafx.h" + #else #include // unix/linux DLL (.so) handling routines #endif @@ -39,7 +39,7 @@ //class K_ORDER_PERTURBATION_API CK_order_perturbation { //public: // CK_order_perturbation(void); - // TODO: add your methods here. +// TODO: add your methods here. //extern void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *g2); @@ -53,33 +53,32 @@ // _Dynamic DLL pointer typedef void * (DynamicFn) - (double *y, double *x, int nb_row_x, double *params, - int it_, double *residual, double *g1, double *g2); +(double *y, double *x, int nb_row_x, double *params, + int it_, double *residual, double *g1, double *g2); typedef void *(mexFunctionPtr)(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]); -//typedef void * (DynamicFn()); +const int MAX_MODEL_NAME=100; class DynamicModelDLL { - -#ifdef WINDOWS - HINSTANCE dynamicHinstance; -# else // linux - void * dynamicHinstance ; -#endif private: - void Dynamic(double *y, double *x, int nb_row_x, double *params, - int it_, double *residual, double *g1, double *g2); -// const DynamicFn (*(DynamicDLLfunc())); + DynamicFn * Dynamic; +#ifdef WINDOWS + HINSTANCE dynamicHinstance; +# else // linux + void * dynamicHinstance ; +#endif +// void Dynamic(double *y, double *x, int nb_row_x, double *params, +// int it_, double *residual, double *g1, double *g2); public: // construct and load Dynamic model DLL DynamicModelDLL(const char* fname); ~DynamicModelDLL(){close();}; -// DynamicFn get(){return DynamicDLLfunc;}; -// void -// ((DynamicFn())*) get(){return Dynamic;}; + // DynamicFn get(){return DynamicDLLfunc;}; + // void + // ((DynamicFn())*) get(){return Dynamic;}; // evaluate Dynamic model DLL void eval(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *g2){ @@ -96,5 +95,5 @@ public: // close DLL: If the referenced object was successfully closed, // close() returns 0, non 0 otherwise int close(); - + };