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

Update - not yet working 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2260 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
george 2008-11-17 17:00:54 +00:00
parent 5fe4ebb670
commit ca0e886204
4 changed files with 569 additions and 445 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

173
mex/sources/korderpert/src/k_ord_dynare.cpp
View File

@ -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<int> DynareNameList::selectIndices(const vector<const char*>& ns) const
{
vector<int> 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<FSSparseTensor>::TensorContainer(dynare.md));
dnl = new DynareNameList(*this);
denl = new DynareExogNameList(*this);
// md= *(new TensorContainer<FSSparseTensor>(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<FSSparseTensor> (1);
// model derivatives FSSparseTensor instance
// model derivatives FSSparseTensor instance
FSSparseTensor mdTi=*(new FSSparseTensor (1, g1->ncols(),g1->nrows()));
for (int i = 0; i<g1->ncols(); i++){
for (int j = 0; j<g1->nrows(); j++){
@ -204,17 +170,15 @@ void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx)
}
}
// md container
md=*(new TensorContainer<FSSparseTensor>::TensorContainer(1));
// md=*(new TensorContainer<FSSparseTensor>(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<int> DynareNameList::selectIndices(const vector<const char*>& ns) const
{
vector<int> 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));
}

98
mex/sources/korderpert/src/k_ord_dynare.h
View File

@ -18,24 +18,40 @@
#ifndef K_ORD_DYNARE3_H
#define K_ORD_DYNARE3_H
#include <vector>
#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<const char*> names;
vector<const char*> 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<int> selectIndices(const vector<const char*>& ns) const;
};
class DynareExogNameList : public NameList {
vector<const char*> 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<const char*> 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<const char*> 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<const char*> 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<FSSparseTensor> 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<FSSparseTensor>& 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);

682
mex/sources/korderpert/src/k_order_perturbation.cpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
* 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 <http://www.gnu.org/licenses/>.
*/
// 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 <mxArray* >(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<int> 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<int> 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 <model>_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<DynamicFn*()>(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);
}
}

61
mex/sources/korderpert/src/k_order_perturbation.h
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
* 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 <http://www.gnu.org/licenses/>.
*/
// 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 <dlfcn.h> // 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 @@
// <model>_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();
};