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

just formatting changes 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2659 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
george 2009-05-06 10:10:27 +00:00
parent aa3b7c78f3
commit 473b4d9a2f
5 changed files with 1982 additions and 2026 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

File diff suppressed because it is too large Load Diff

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

@ -1,283 +1,250 @@
/*
* Copyright (C) 2005 Ondra Kamenik
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* 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/>.
*/
// based on: $Id: dynare3.h 1764 2008-03-31 14:30:55Z kamenik $
// by 2005, Ondra Kamenik
#ifndef K_ORD_DYNARE3_H
#define K_ORD_DYNARE3_H
#include <vector>
#include "t_container.h"
#include "sparse_tensor.h"
#include "decision_rule.h"
#include "dynamic_model.h"
#include "exception.h"
#include "dynare_exception.h"
#include "fs_tensor.h"
#include "SylvException.h"
#include "tl_exception.h"
#include "kord_exception.h"
#include "nlsolve.h"
#include "approximation.h"
#include "k_order_perturbation.h"
class KordpDynare;
// derive from Approximation to get protected derivatives
/******
class FistOrderApproximation: public Approximation{
TwoDMatrix *gy;
TwoDMatrix *gu;
public:
FistOrderApproximation();
FistOrderApproximation(FistOrderApproximation& fo): ::Approximation(fo){
if (&(fo.GetGy())!=0){
gy= new TwoDMatrix(fo.GetGy());
gu= new TwoDMatrix(fo.GetGu());
}
};
virtual ~FistOrderApproximation(){
delete gy;
delete gu;
};
virtual void approxAtSteady(); // NOTE: change the Approximation parent to use virual too so that it can be overriden !!
//FGSContainer* GetRuleDers(){return rule_ders;};
//FGSContainer* GetRuleDersSS(){return rule_ders_ss;};
const TwoDMatrix& GetGy(){return (const TwoDMatrix&)(*gy);};
const TwoDMatrix& GetGu(){return (const TwoDMatrix&)(*gu);};
virtual void saveRuleDerivs(FirstOrder& fo);
};
***********/
/*////////////////////////////////////////////*/
// instantiations of pure abstract class NameList in dynamic_model.h:
/*////////////////////////////////////////////*/
class DynareNameList : public NameList {
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
* in the names. And returns the resulting vector of indices. If
* the name cannot be found, then an exception is raised. */
vector<int> selectIndices(const vector<const char*>& ns) const;
};
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];}
};
/*********************************************/
// The following only implements DynamicModel with help of ogdyn::DynareModel
// instantiation of pure abstract DynamicModel decl. in dynamic_model.h
//class DynamicModelDLL;
class KordpJacobian;
class KordpDynare : public DynamicModel {
friend class DynareNameList;
friend class DynareExogNameList;
friend class DynareStateNameList;
friend class KordpDynareJacobian;
friend class DynamicModelDLL;
//////////
const int nStat;
const int nBoth;
const int nPred;
const int nForw;
const int nExog;
const int nPar;
const int nYs; // ={npred + nboth ; }
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 int nSteps;
const int nOrder;
Journal& journal;
/// DynamicModel* model;
///const char* modName;
Vector* ySteady;
Vector* params;
TwoDMatrix* vCov;
TensorContainer<FSSparseTensor> md; // ModelDerivatives
DynareNameList* dnl;
DynareExogNameList* denl;
DynareStateNameList* dsnl;
const double ss_tol;
const vector<int>* varOrder;
const TwoDMatrix * ll_Incidence;
double qz_criterium;
vector<int> * JacobianIndices;
public:
KordpDynare(const char** endo, int num_endo,
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 nJcols, const int nSteps, const int ord, //const char* modName,
Journal& jr, DynamicModelDLL& dynamicDLL, double sstol,
const vector<int>* varOrder, const TwoDMatrix * ll_Incidence,
double qz_criterium );
/** Makes a deep copy of the object. */
KordpDynare(const KordpDynare& dyn);
virtual ~KordpDynare();
int nstat() const
{return nStat;}
int nboth() const
{return nBoth;}
int npred() const
{return nPred;}
int nforw() const
{return nForw;}
int nexog() const
{return nExog;}
int nys() const
{return nYs;}
int nyss() const
{return nYss;}
int ny() const
{return nY;}
int steps() const
{return nSteps;}
int order() const
{return nOrder;}
const NameList& getAllEndoNames() const
{return *dnl;}
const NameList& getStateNames() const
{return *dsnl;}
const NameList& getExogNames() const
{return *denl;}
const TwoDMatrix& getVcov() const
{return *vCov;}
Vector& getParams()
{return *params;}
const TensorContainer<FSSparseTensor>& getModelDerivatives() const
{return md;}
const Vector& getSteady() const
{return *ySteady;}
Vector& getSteady()
{return *ySteady;}
/// const ogdyn::DynareModel& getModel() const
/// {return *model;}
// here is true public interface
void solveDeterministicSteady(Vector& steady);
void solveDeterministicSteady()
{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);
void calcDerivatives(const Vector& yy, const Vector& xx);
//void calcDerivatives(const Vector& yy, TwoDMatrix& jj);
void calcDerivatives(const Vector& yy, ogu::Jacobian& jacob);
void calcDerivativesAtSteady();
DynamicModelDLL& dynamicDLL;
/// void writeMat4(FILE* fd, const char* prefix) const;
/// void writeDump(const std::string& basename) const;
DynamicModel* clone() const
{return new KordpDynare(*this);}
void ReorderCols(TwoDMatrix * tdx, const int * varOrder);
void ReorderCols(TwoDMatrix * tdx, const vector<int> * varOrder);
Vector * LLxSteady( const Vector& yS); // returns ySteady extended with leads and lags
private:
void writeModelInfo(Journal& jr) const;
int * ReorderDynareJacobianIndices( const int * varOrder);
vector<int> * ReorderDynareJacobianIndices( const vector<int> * varOrder);
void ReorderBlocks(TwoDMatrix * tdx, const int * varOrder);
void ReorderBlocks(TwoDMatrix * tdx, const vector<int> * vOrder);
void populateDerivativesContainer(TwoDMatrix*g, int ord, const vector<int>* vOrder);
};
/****************************
* ModelDerivativeContainer manages derivatives container
************************************/
class ModelDerivativeContainer //: public ogp::FormulaDerEvalLoader
{
protected:
// const ogp::FineAtoms& atoms;
TensorContainer<FSSparseTensor>& md;
public:
ModelDerivativeContainer(const KordpDynare& model, TensorContainer<FSSparseTensor>& mod_ders,
int order);
void load(int i, int iord, const int* vars, double res);
};
/****************************
* K-Order Perturbation instance of Jacobian:
************************************/
class KordpJacobian : public ogu::Jacobian ///, public ogp::FormulaDerEvalLoader
{
protected:
KordpDynare& dyn;
public:
KordpJacobian( KordpDynare& dyn);
virtual ~KordpJacobian() {}
// Load <mod>_dynamic.DLL
// void load(const char** modName);
void eval(const Vector& in);
};
/****************************
* K-Order Perturbation instance of VectorFunction:
************************************/
class KordpVectorFunction : public ogu::VectorFunction {
protected:
KordpDynare& d;
public:
KordpVectorFunction( KordpDynare& dyn)
: d(dyn) {}
virtual ~KordpVectorFunction() {}
int inDim() const
{return d.ny();}
int outDim() const
{return d.ny();}
void eval(const ConstVector& in, Vector& out);
};
#endif
// Local Variables:
// mode:C++
// End:
/*
* Copyright (C) 2005 Ondra Kamenik
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* 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/>.
*/
// based on: $Id: dynare3.h 1764 2008-03-31 14:30:55Z kamenik $
// by 2005, Ondra Kamenik
#ifndef K_ORD_DYNARE3_H
#define K_ORD_DYNARE3_H
#include <vector>
#include "t_container.h"
#include "sparse_tensor.h"
#include "decision_rule.h"
#include "dynamic_model.h"
#include "exception.h"
#include "dynare_exception.h"
#include "fs_tensor.h"
#include "SylvException.h"
#include "tl_exception.h"
#include "kord_exception.h"
#include "nlsolve.h"
#include "approximation.h"
#include "k_order_perturbation.h"
class KordpDynare;
/*////////////////////////////////////////////*/
// instantiations of pure abstract class NameList in dynamic_model.h:
/*////////////////////////////////////////////*/
class DynareNameList : public NameList {
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
* in the names. And returns the resulting vector of indices. If
* the name cannot be found, then an exception is raised. */
vector<int> selectIndices(const vector<const char*>& ns) const;
};
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];}
};
/*********************************************/
// The following only implements DynamicModel with help of ogdyn::DynareModel
// instantiation of pure abstract DynamicModel decl. in dynamic_model.h
//class DynamicModelDLL;
class KordpJacobian;
class KordpDynare : public DynamicModel {
friend class DynareNameList;
friend class DynareExogNameList;
friend class DynareStateNameList;
friend class KordpDynareJacobian;
friend class DynamicModelDLL;
//////////
const int nStat;
const int nBoth;
const int nPred;
const int nForw;
const int nExog;
const int nPar;
const int nYs; // ={npred + nboth ; }
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 int nSteps;
const int nOrder;
Journal& journal;
/// DynamicModel* model;
///const char* modName;
Vector* ySteady;
Vector* params;
TwoDMatrix* vCov;
TensorContainer<FSSparseTensor> md; // ModelDerivatives
DynareNameList* dnl;
DynareExogNameList* denl;
DynareStateNameList* dsnl;
const double ss_tol;
const vector<int>* varOrder;
const TwoDMatrix * ll_Incidence;
double qz_criterium;
vector<int> * JacobianIndices;
public:
KordpDynare(const char** endo, int num_endo,
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 nJcols, const int nSteps, const int ord, //const char* modName,
Journal& jr, DynamicModelDLL& dynamicDLL, double sstol,
const vector<int>* varOrder, const TwoDMatrix * ll_Incidence,
double qz_criterium );
/** Makes a deep copy of the object. */
KordpDynare(const KordpDynare& dyn);
virtual ~KordpDynare();
int nstat() const
{return nStat;}
int nboth() const
{return nBoth;}
int npred() const
{return nPred;}
int nforw() const
{return nForw;}
int nexog() const
{return nExog;}
int nys() const
{return nYs;}
int nyss() const
{return nYss;}
int ny() const
{return nY;}
int steps() const
{return nSteps;}
int order() const
{return nOrder;}
const NameList& getAllEndoNames() const
{return *dnl;}
const NameList& getStateNames() const
{return *dsnl;}
const NameList& getExogNames() const
{return *denl;}
const TwoDMatrix& getVcov() const
{return *vCov;}
Vector& getParams()
{return *params;}
const TensorContainer<FSSparseTensor>& getModelDerivatives() const
{return md;}
const Vector& getSteady() const
{return *ySteady;}
Vector& getSteady()
{return *ySteady;}
// here is true public interface
void solveDeterministicSteady(Vector& steady);
void solveDeterministicSteady()
{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);
void calcDerivatives(const Vector& yy, const Vector& xx);
//void calcDerivatives(const Vector& yy, TwoDMatrix& jj);
void calcDerivatives(const Vector& yy, ogu::Jacobian& jacob);
void calcDerivativesAtSteady();
DynamicModelDLL& dynamicDLL;
/// void writeMat4(FILE* fd, const char* prefix) const;
/// void writeDump(const std::string& basename) const;
DynamicModel* clone() const
{return new KordpDynare(*this);}
void ReorderCols(TwoDMatrix * tdx, const int * varOrder);
void ReorderCols(TwoDMatrix * tdx, const vector<int> * varOrder);
Vector * LLxSteady( const Vector& yS); // returns ySteady extended with leads and lags
private:
void writeModelInfo(Journal& jr) const;
int * ReorderDynareJacobianIndices( const int * varOrder);
vector<int> * ReorderDynareJacobianIndices( const vector<int> * varOrder);
void ReorderBlocks(TwoDMatrix * tdx, const int * varOrder);
void ReorderBlocks(TwoDMatrix * tdx, const vector<int> * vOrder);
void populateDerivativesContainer(TwoDMatrix*g, int ord, const vector<int>* vOrder);
};
/****************************
* ModelDerivativeContainer manages derivatives container
************************************/
class ModelDerivativeContainer //: public ogp::FormulaDerEvalLoader
{
protected:
// const ogp::FineAtoms& atoms;
TensorContainer<FSSparseTensor>& md;
public:
ModelDerivativeContainer(const KordpDynare& model, TensorContainer<FSSparseTensor>& mod_ders,
int order);
void load(int i, int iord, const int* vars, double res);
};
/****************************
* K-Order Perturbation instance of Jacobian:
************************************/
class KordpJacobian : public ogu::Jacobian ///, public ogp::FormulaDerEvalLoader
{
protected:
KordpDynare& dyn;
public:
KordpJacobian( KordpDynare& dyn);
virtual ~KordpJacobian() {}
// Load <mod>_dynamic.DLL
// void load(const char** modName);
void eval(const Vector& in);
};
/****************************
* K-Order Perturbation instance of VectorFunction:
************************************/
class KordpVectorFunction : public ogu::VectorFunction {
protected:
KordpDynare& d;
public:
KordpVectorFunction( KordpDynare& dyn)
: d(dyn) {}
virtual ~KordpVectorFunction() {}
int inDim() const
{return d.ny();}
int outDim() const
{return d.ny();}
void eval(const ConstVector& in, Vector& out);
};
#endif

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

File diff suppressed because it is too large Load Diff

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

@ -1,110 +1,108 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* 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/>.
*/
// 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
// symbol defined on the command line. this symbol should not be defined on any project
// that uses this DLL. This way any other project whose source files include this file see
// 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
#include "stdafx.h"
#else
#include <dlfcn.h> // unix/linux DLL (.so) handling routines
#endif
#include <string>
#include "mex.h"
// <model>_Dynamic DLL pointer
#ifdef WINDOWS
typedef void *(DynamicFn)
#else // linux
typedef void (*DynamicFn)
#endif
(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2);
//DynamicFn Dynamic;
typedef void *(mexFunctionPtr)(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);
const int MAX_MODEL_NAME=100;
/**
* creates pointer to Dynamic function inside <model>_dynamic.dll
* and handles calls to it.
**/
class DynamicModelDLL
{
private:
#ifdef WINDOWS
DynamicFn *Dynamic;// pointer to the Dynamic function in DLL
#else
DynamicFn Dynamic;// pointer to the Dynamic function in DLL
#endif
const int length; // tot num vars = Num of Jacobian rows
const int jcols; // tot num var t-1, t and t+1 instances + exogs = Num of Jacobian columns
const int nMax_lag; // no of lags
const int nExog; // no of exogenous
// const char * sExt; // dynamic file extension: dll, mexw32, ...
#ifdef WINDOWS
HINSTANCE dynamicHinstance; // DLL instance pointer in Windows
# else // linux
void * dynamicHinstance ; // and in Linux
#endif
public:
// construct and load Dynamic model DLL
DynamicModelDLL(const char* fname, const int length,const int jcols,
const int nMax_lag, const int nExog, const char * sExt);
virtual ~DynamicModelDLL(){close();};
// 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){
Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);
};
void eval(const Vector&y, const Vector&x, const Vector* params,
Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
void eval(const Vector&y, const TwoDMatrix&x, const Vector* params,
int it_, Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
void eval(const Vector&y, const TwoDMatrix&x, const Vector* params,
Vector& residual, TwoDMatrix *g1, TwoDMatrix *g2);
// void eval(const Vector&y, const TwoDMatrix&x, const Vector* params,
// Vector& residual, double *g1, double *g2);
// close DLL: If the referenced object was successfully closed,
// close() returns 0, non 0 otherwise
int close();
};
// convert Matlab endo and exo names array to C type array of strings
const char ** DynareMxArrayToString(const mxArray * mxFldp, const int len, const int width );
const char ** DynareMxArrayToString(const char * cArray, const int len, const int width );
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* 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/>.
*/
// 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
// symbol defined on the command line. this symbol should not be defined on any project
// that uses this DLL. This way any other project whose source files include this file see
// 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
#include "stdafx.h"
#else
#include <dlfcn.h> // unix/linux DLL (.so) handling routines
#endif
#include <string>
#include "mex.h"
// <model>_Dynamic DLL pointer
#ifdef WINDOWS
typedef void *(DynamicFn)
#else // linux
typedef void (*DynamicFn)
#endif
(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2);
//DynamicFn Dynamic;
typedef void *(mexFunctionPtr)(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);
const int MAX_MODEL_NAME=100;
/**
* creates pointer to Dynamic function inside <model>_dynamic.dll
* and handles calls to it.
**/
class DynamicModelDLL
{
private:
#ifdef WINDOWS
DynamicFn *Dynamic;// pointer to the Dynamic function in DLL
#else
DynamicFn Dynamic;// pointer to the Dynamic function in DLL
#endif
const int length; // tot num vars = Num of Jacobian rows
const int jcols; // tot num var t-1, t and t+1 instances + exogs = Num of Jacobian columns
const int nMax_lag; // no of lags
const int nExog; // no of exogenous
// const char * sExt; // dynamic file extension: dll, mexw32, ...
#ifdef WINDOWS
HINSTANCE dynamicHinstance; // DLL instance pointer in Windows
# else // linux
void * dynamicHinstance ; // and in Linux
#endif
public:
// construct and load Dynamic model DLL
DynamicModelDLL(const char* fname, const int length,const int jcols,
const int nMax_lag, const int nExog, const char * sExt);
virtual ~DynamicModelDLL(){close();};
// evaluate Dynamic model DLL
void eval(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2){
Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);
};
void eval(const Vector&y, const Vector&x, const Vector* params,
Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
void eval(const Vector&y, const TwoDMatrix&x, const Vector* params,
int it_, Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
void eval(const Vector&y, const TwoDMatrix&x, const Vector* params,
Vector& residual, TwoDMatrix *g1, TwoDMatrix *g2);
// void eval(const Vector&y, const TwoDMatrix&x, const Vector* params,
// Vector& residual, double *g1, double *g2);
// close DLL: If the referenced object was successfully closed,
// close() returns 0, non 0 otherwise
int close();
};
// convert Matlab endo and exo names array to C type array of strings
const char ** DynareMxArrayToString(const mxArray * mxFldp, const int len, const int width );
const char ** DynareMxArrayToString(const char * cArray, const int len, const int width );

648
mex/sources/korderpert/src/k_order_test_main.cpp
View File

@ -1,324 +1,324 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* 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/>.
*/
/*************************************
* This main() is for testing k_order DLL entry point by linking to
* the k_ord library statically and passing its hard-coded data:
* parameters, covar, ysteady and the variable names from fs2000a.mod model
* The main has been derived from mxFunction used for K-Order DLL
***************************************/
//#include "stdafx.h"
#include "k_ord_dynare.h"
int main(int argc, char* argv[])
{
double qz_criterium = 1+1e-6;
const int check_flag = 0;
const char* fName = "fs2000k";//mxArrayToString(mFname);
const char * dfExt=NULL ;//Dyanamic file extension, e.g.".dll";
#ifdef DEBUG
mexPrintf("k_order_perturbation: check_flag = %d , fName = %s .\n", check_flag,fName);
#endif
int kOrder =2;
int npar = 7;//(int)mxGetM(mxFldp);
double dparams[7]={ 0.3300,
0.9900,
0.0030,
1.0110,
0.7000,
0.7870,
0.0200
};
Vector * modParams = new Vector(dparams, npar);
#ifdef DEBUG
mexPrintf("k_ord_perturbation: nParams=%d .\n",npar);
for (int i = 0; i < npar; i++) {
mexPrintf("k_ord_perturbation: dParams[%d]= %g.\n", i, dparams+i*(sizeof(double)) ); }
for (int i = 0; i < npar; i++) {
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); }
#endif
double d2Dparams[4] = { //(double *) mxGetData(mxFldp);
0.1960e-3, 0.0,
0.0, 0.0250e-3};
npar = 2;//(int)mxGetN(mxFldp);
TwoDMatrix * vCov = new TwoDMatrix(npar, npar, (d2Dparams));
double dYSparams [16]= { // 27 mxGetData(mxFldp);
// 1.0110, 2.2582, 5.8012, 0.5808,
1.0110, 2.2582, 0.4477, 1.0000
, 4.5959, 1.0212, 5.8012, 0.8494
, 0.1872, 0.8604, 1.0030, 1.0080
, 0.5808, 1.0030, 2.2582, 0.4477
//, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477
};
const int nSteady = 16;//27 //31;//29, 16 (int)mxGetM(mxFldp);
Vector * ySteady = new Vector(dYSparams, nSteady);
//mxFldp = mxGetField(dr, 0,"nstatic" );
const int nStat = 7;//(int)mxGetScalar(mxFldp);
// mxFldp = mxGetField(dr, 0,"npred" );
const int nPred = 2;//6 - nBoth (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nspred" );
const int nsPred = 4;//(int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nboth" );
const int nBoth = 2;// (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nfwrd" );
const int nForw = 5;// 3 (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nsfwrd" );
const int nsForw = 7;//(int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(M_, 0,"exo_nbr" );
const int nExog =2;// (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(M_, 0,"endo_nbr" );
const int nEndo = 16;//16(int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(M_, 0,"param_nbr" );
const int nPar = 7;//(int)mxGetScalar(mxFldp);
// it_ should be set to M_.maximum_lag
//mxFldp = mxGetField(M_, 0,"maximum_lag" );
const int nMax_lag = 1;//(int)mxGetScalar(mxFldp);
int var_order[]//[18]
= {
5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16
// 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18
};
//Vector * varOrder = new Vector(var_order, nEndo);
vector<int> * var_order_vp = new vector<int>(nEndo);//nEndo));
for (int v=0;v<nEndo;v++)
(*var_order_vp)[v] =var_order[v];
const double ll_incidence []//[3][18]
= {
1, 5, 21
, 2, 6, 22
, 0, 7, 23
, 0, 8, 24
, 0, 9, 0
, 0, 10, 0
, 3, 11, 0
, 0, 12, 0
, 0, 13, 25
, 0, 14, 0
, 0, 15, 0
, 0, 16, 0
, 4, 17, 0
, 0, 18, 0
, 0, 19, 26
, 0, 20, 27
};
TwoDMatrix * llincidence = new TwoDMatrix ( 3, nEndo, ll_incidence );
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
#ifdef DEBUG
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
#endif
//mxFldp= mxGetField(M_, 0,"endo_names" );
const int nendo = 16;//16(int)mxGetM(mxFldp);
const int widthEndo = 6;// (int)mxGetN(mxFldp);
const char * cNamesCharStr= "mPceWRkdnlggydPc yp A22 __ oo bb ss ";
// const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo);
const char** endoNamesMX= DynareMxArrayToString( cNamesCharStr, nendo, widthEndo);
#ifdef DEBUG
for (int i = 0; i < nEndo; i++) {
mexPrintf("k_ord_perturbation: EndoNameList[%d][0]= %s.\n", i, endoNamesMX[i] );
}
#endif
//mxFldp = mxGetField(M_, 0,"exo_names" );
const int nexo = 2;//(int)mxGetM(mxFldp);
const int widthExog = 3;//(int)mxGetN(mxFldp);
// const char** exoNamesMX= DynareMxArrayToString( mxFldp,nexo,widthExog);
const char * cExoNamesCharStr= "ee__am";
const char** exoNamesMX= DynareMxArrayToString( cExoNamesCharStr,nexo,widthExog);
#ifdef DEBUG
for (int i = 0; i < nexo; i++) {
mexPrintf("k_ord_perturbation: ExoNameList[%d][0]= %s.\n", i, exoNamesMX[i] );
}
#endif
if ((nEndo != nendo)||(nExog != nexo)) { //(nPar != npar)
mexErrMsgTxt("Incorrect number of input parameters.\n");
//return;
}
#ifdef DEBUG
for (int i = 0; i < nEndo; i++) {
mexPrintf("k_ord_perturbation: EndoNameList[%d]= %s.\n", i, endoNamesMX[i] ); }
// for (int i = 0; i < nPar; i++) {
//, , mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); }
for (int i = 0; i < nPar; i++) {
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); }
for (int i = 0; i < nSteady; i++) {
mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]); }
mexPrintf("k_order_perturbation: 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 = 1;//2 params.num_threads;
try {
// make journal name and journal
std::string jName(fName); //params.basename);
jName += ".jnl";
Journal journal(jName.c_str());
#ifdef DEBUG
mexPrintf("k_order_perturbation: Call tls init\n");
#endif
tls.init(kOrder, (nStat+2*nPred+3*nBoth+2*nForw+nExog));
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynamicDLL constructor.\n");
#endif
// DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
DynamicModelDLL dynamicDLL(fName, nEndo, jcols, nMax_lag, nExog, dfExt);
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynare constructor.\n");
#endif
// make KordpDynare object
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar // paramNames,
, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth
, jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp //var_order
, llincidence, qz_criterium);
#ifdef DEBUG
mexPrintf("k_order_perturbation: Call Approximation constructor \n");
#endif
Approximation app(dynare, journal, nSteps, false, qz_criterium);
// run stochastic steady
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling walkStochSteady.\n");
#endif
app.walkStochSteady();
// open mat file
std::string matfile(fName);//(params.basename);
matfile += ".mat";
FILE* matfd = NULL;
if (NULL == (matfd=fopen(matfile.c_str(), "wb"))) {
fprintf(stderr, "Couldn't open %s for writing.\n", matfile.c_str());
exit(1);
}
#ifdef DEBUG
mexPrintf("k_order_perturbation: Filling Mat file outputs.\n");
#endif
std::string ss_matrix_name(fName);//params.prefix);
ss_matrix_name += "_steady_states";
ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str());
// write the folded decision rule to the Mat-4 file
app.getFoldDecisionRule().writeMat4(matfd, fName);//params.prefix);
fclose(matfd);
map<string,ConstTwoDMatrix> mm;
app.getFoldDecisionRule().writeMMap(&mm);
#ifdef DEBUG
app.getFoldDecisionRule().print();
mexPrintf("k_order_perturbation: Map print: \n");
for (map<string,ConstTwoDMatrix>::const_iterator cit=mm.begin();
cit !=mm.end(); ++cit) {
mexPrintf("k_order_perturbation: Map print: string: %s , g:\n", (*cit).first.c_str());
(*cit).second.print();
}
#endif
// get latest ysteady
double * dYsteady = (dynare.getSteady().base());
ySteady = (Vector*)(&dynare.getSteady());
} catch (const KordException& e) {
printf("Caugth Kord exception: ");
e.print();
return 1;// e.code();
} catch (const TLException& e) {
printf("Caugth TL exception: ");
e.print();
return 2;// 255;
} catch (SylvException& e) {
printf("Caught Sylv exception: ");
e.printMessage();
return 3;// 255;
} catch (const DynareException& e) {
printf("Caught KordpDynare exception: %s\n", e.message());
return 4;// 255;
} catch (const ogu::Exception& e) {
printf("Caught ogu::Exception: ");
e.print();
return 5;// 255;
}
// bones for future developement of the Matlab output.
const int nrhs=5;
const int nlhs=2;
mxArray* prhs[nrhs];
mxArray* plhs[nlhs];
#ifdef DEBUG
mexPrintf("k_order_perturbation: Filling Matlab outputs.\n");
#endif
double *dgy, *dgu, *ysteady;
int nb_row_x;
ysteady = NULL;
if (nlhs >= 1)
{
/* Set the output pointer to the output matrix ysteady. */
plhs[0] = mxCreateDoubleMatrix(nEndo,1, mxREAL);
/* Create a C pointer to a copy of the output ysteady. */
ysteady = mxGetPr(plhs[0]);
}
dgy = NULL;
if (nlhs >= 2)
{
/* Set the output pointer to the output matrix gy. */
plhs[1] = mxCreateDoubleMatrix(nEndo, jcols, mxREAL);
// plhs[1] = (double*)(gy->getData())->base();
/* Create a C pointer to a copy of the output matrix gy. */
dgy = mxGetPr(plhs[1]);
}
dgu = NULL;
if (nlhs >= 3)
{
/* Set the output pointer to the output matrix gu. */
plhs[2] = mxCreateDoubleMatrix(nEndo, nExog, mxREAL);
// plhs[2] = (double*)((gu->getData())->base());
/* Create a C pointer to a copy of the output matrix gu. */
dgu = mxGetPr(plhs[2]);
}
return 0;
}
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* 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/>.
*/
/*************************************
* This main() is for testing k_order DLL entry point by linking to
* the k_ord library statically and passing its hard-coded data:
* parameters, covar, ysteady and the variable names from fs2000a.mod model
* The main has been derived from mxFunction used for K-Order DLL
***************************************/
//#include "stdafx.h"
#include "k_ord_dynare.h"
int main(int argc, char* argv[])
{
double qz_criterium = 1+1e-6;
const int check_flag = 0;
const char* fName = "fs2000k";//mxArrayToString(mFname);
const char * dfExt=NULL ;//Dyanamic file extension, e.g.".dll";
#ifdef DEBUG
mexPrintf("k_order_perturbation: check_flag = %d , fName = %s .\n", check_flag,fName);
#endif
int kOrder =2;
int npar = 7;//(int)mxGetM(mxFldp);
double dparams[7]={ 0.3300,
0.9900,
0.0030,
1.0110,
0.7000,
0.7870,
0.0200
};
Vector * modParams = new Vector(dparams, npar);
#ifdef DEBUG
mexPrintf("k_ord_perturbation: nParams=%d .\n",npar);
for (int i = 0; i < npar; i++) {
mexPrintf("k_ord_perturbation: dParams[%d]= %g.\n", i, dparams+i*(sizeof(double)) ); }
for (int i = 0; i < npar; i++) {
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); }
#endif
double d2Dparams[4] = { //(double *) mxGetData(mxFldp);
0.1960e-3, 0.0,
0.0, 0.0250e-3};
npar = 2;//(int)mxGetN(mxFldp);
TwoDMatrix * vCov = new TwoDMatrix(npar, npar, (d2Dparams));
double dYSparams [16]= { // 27 mxGetData(mxFldp);
// 1.0110, 2.2582, 5.8012, 0.5808,
1.0110, 2.2582, 0.4477, 1.0000
, 4.5959, 1.0212, 5.8012, 0.8494
, 0.1872, 0.8604, 1.0030, 1.0080
, 0.5808, 1.0030, 2.2582, 0.4477
//, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477
};
const int nSteady = 16;//27 //31;//29, 16 (int)mxGetM(mxFldp);
Vector * ySteady = new Vector(dYSparams, nSteady);
//mxFldp = mxGetField(dr, 0,"nstatic" );
const int nStat = 7;//(int)mxGetScalar(mxFldp);
// mxFldp = mxGetField(dr, 0,"npred" );
const int nPred = 2;//6 - nBoth (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nspred" );
const int nsPred = 4;//(int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nboth" );
const int nBoth = 2;// (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nfwrd" );
const int nForw = 5;// 3 (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(dr, 0,"nsfwrd" );
const int nsForw = 7;//(int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(M_, 0,"exo_nbr" );
const int nExog =2;// (int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(M_, 0,"endo_nbr" );
const int nEndo = 16;//16(int)mxGetScalar(mxFldp);
//mxFldp = mxGetField(M_, 0,"param_nbr" );
const int nPar = 7;//(int)mxGetScalar(mxFldp);
// it_ should be set to M_.maximum_lag
//mxFldp = mxGetField(M_, 0,"maximum_lag" );
const int nMax_lag = 1;//(int)mxGetScalar(mxFldp);
int var_order[]//[18]
= {
5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16
// 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18
};
//Vector * varOrder = new Vector(var_order, nEndo);
vector<int> * var_order_vp = new vector<int>(nEndo);//nEndo));
for (int v=0;v<nEndo;v++)
(*var_order_vp)[v] =var_order[v];
const double ll_incidence []//[3][18]
= {
1, 5, 21
, 2, 6, 22
, 0, 7, 23
, 0, 8, 24
, 0, 9, 0
, 0, 10, 0
, 3, 11, 0
, 0, 12, 0
, 0, 13, 25
, 0, 14, 0
, 0, 15, 0
, 0, 16, 0
, 4, 17, 0
, 0, 18, 0
, 0, 19, 26
, 0, 20, 27
};
TwoDMatrix * llincidence = new TwoDMatrix ( 3, nEndo, ll_incidence );
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
#ifdef DEBUG
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
#endif
//mxFldp= mxGetField(M_, 0,"endo_names" );
const int nendo = 16;//16(int)mxGetM(mxFldp);
const int widthEndo = 6;// (int)mxGetN(mxFldp);
const char * cNamesCharStr= "mPceWRkdnlggydPc yp A22 __ oo bb ss ";
// const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo);
const char** endoNamesMX= DynareMxArrayToString( cNamesCharStr, nendo, widthEndo);
#ifdef DEBUG
for (int i = 0; i < nEndo; i++) {
mexPrintf("k_ord_perturbation: EndoNameList[%d][0]= %s.\n", i, endoNamesMX[i] );
}
#endif
//mxFldp = mxGetField(M_, 0,"exo_names" );
const int nexo = 2;//(int)mxGetM(mxFldp);
const int widthExog = 3;//(int)mxGetN(mxFldp);
// const char** exoNamesMX= DynareMxArrayToString( mxFldp,nexo,widthExog);
const char * cExoNamesCharStr= "ee__am";
const char** exoNamesMX= DynareMxArrayToString( cExoNamesCharStr,nexo,widthExog);
#ifdef DEBUG
for (int i = 0; i < nexo; i++) {
mexPrintf("k_ord_perturbation: ExoNameList[%d][0]= %s.\n", i, exoNamesMX[i] );
}
#endif
if ((nEndo != nendo)||(nExog != nexo)) { //(nPar != npar)
mexErrMsgTxt("Incorrect number of input parameters.\n");
//return;
}
#ifdef DEBUG
for (int i = 0; i < nEndo; i++) {
mexPrintf("k_ord_perturbation: EndoNameList[%d]= %s.\n", i, endoNamesMX[i] ); }
// for (int i = 0; i < nPar; i++) {
//, , mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); }
for (int i = 0; i < nPar; i++) {
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); }
for (int i = 0; i < nSteady; i++) {
mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]); }
mexPrintf("k_order_perturbation: 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 = 1;//2 params.num_threads;
try {
// make journal name and journal
std::string jName(fName); //params.basename);
jName += ".jnl";
Journal journal(jName.c_str());
#ifdef DEBUG
mexPrintf("k_order_perturbation: Call tls init\n");
#endif
tls.init(kOrder, (nStat+2*nPred+3*nBoth+2*nForw+nExog));
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynamicDLL constructor.\n");
#endif
// DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
DynamicModelDLL dynamicDLL(fName, nEndo, jcols, nMax_lag, nExog, dfExt);
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynare constructor.\n");
#endif
// make KordpDynare object
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar // paramNames,
, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth
, jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp //var_order
, llincidence, qz_criterium);
#ifdef DEBUG
mexPrintf("k_order_perturbation: Call Approximation constructor \n");
#endif
Approximation app(dynare, journal, nSteps, false, qz_criterium);
// run stochastic steady
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling walkStochSteady.\n");
#endif
app.walkStochSteady();
// open mat file
std::string matfile(fName);//(params.basename);
matfile += ".mat";
FILE* matfd = NULL;
if (NULL == (matfd=fopen(matfile.c_str(), "wb"))) {
fprintf(stderr, "Couldn't open %s for writing.\n", matfile.c_str());
exit(1);
}
#ifdef DEBUG
mexPrintf("k_order_perturbation: Filling Mat file outputs.\n");
#endif
std::string ss_matrix_name(fName);//params.prefix);
ss_matrix_name += "_steady_states";
ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str());
// write the folded decision rule to the Mat-4 file
app.getFoldDecisionRule().writeMat4(matfd, fName);//params.prefix);
fclose(matfd);
map<string,ConstTwoDMatrix> mm;
app.getFoldDecisionRule().writeMMap(&mm);
#ifdef DEBUG
app.getFoldDecisionRule().print();
mexPrintf("k_order_perturbation: Map print: \n");
for (map<string,ConstTwoDMatrix>::const_iterator cit=mm.begin();
cit !=mm.end(); ++cit) {
mexPrintf("k_order_perturbation: Map print: string: %s , g:\n", (*cit).first.c_str());
(*cit).second.print();
}
#endif
// get latest ysteady
double * dYsteady = (dynare.getSteady().base());
ySteady = (Vector*)(&dynare.getSteady());
} catch (const KordException& e) {
printf("Caugth Kord exception: ");
e.print();
return 1;// e.code();
} catch (const TLException& e) {
printf("Caugth TL exception: ");
e.print();
return 2;// 255;
} catch (SylvException& e) {
printf("Caught Sylv exception: ");
e.printMessage();
return 3;// 255;
} catch (const DynareException& e) {
printf("Caught KordpDynare exception: %s\n", e.message());
return 4;// 255;
} catch (const ogu::Exception& e) {
printf("Caught ogu::Exception: ");
e.print();
return 5;// 255;
}
// bones for future developement of the Matlab output.
const int nrhs=5;
const int nlhs=2;
mxArray* prhs[nrhs];
mxArray* plhs[nlhs];
#ifdef DEBUG
mexPrintf("k_order_perturbation: Filling Matlab outputs.\n");
#endif
double *dgy, *dgu, *ysteady;
int nb_row_x;
ysteady = NULL;
if (nlhs >= 1)
{
/* Set the output pointer to the output matrix ysteady. */
plhs[0] = mxCreateDoubleMatrix(nEndo,1, mxREAL);
/* Create a C pointer to a copy of the output ysteady. */
ysteady = mxGetPr(plhs[0]);
}
dgy = NULL;
if (nlhs >= 2)
{
/* Set the output pointer to the output matrix gy. */
plhs[1] = mxCreateDoubleMatrix(nEndo, jcols, mxREAL);
// plhs[1] = (double*)(gy->getData())->base();
/* Create a C pointer to a copy of the output matrix gy. */
dgy = mxGetPr(plhs[1]);
}
dgu = NULL;
if (nlhs >= 3)
{
/* Set the output pointer to the output matrix gu. */
plhs[2] = mxCreateDoubleMatrix(nEndo, nExog, mxREAL);
// plhs[2] = (double*)((gu->getData())->base());
/* Create a C pointer to a copy of the output matrix gu. */
dgu = mxGetPr(plhs[2]);
}
return 0;
}