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

v4.1: fixing k_order_perturbation DLL for order == 3 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3170 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2009-11-29 20:50:39 +00:00
parent 494668a167
commit 5f05a14b17
4 changed files with 188 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

95
mex/sources/k_order_perturbation/dynamic_dll.cpp
View File

@ -23,60 +23,6 @@
#include "math.h" #include "math.h"
#include <cstring> #include <cstring>
//////////////////////////////////////////////////////
// Convert Matlab Dynare endo and exo names array to C type array of string pointers
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
// hence a rather low level approach is needed
///////////////////////////////////////////////////////
const char **
DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width)
{
char *cNamesCharStr = mxArrayToString(mxFldp);
const char **ret = DynareMxArrayToString(cNamesCharStr, len, width);
return ret;
}
const char **
DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width)
{
char **cNamesMX;
cNamesMX = (char **) mxCalloc(len, sizeof(char *));
for(int i = 0; i < len; i++)
cNamesMX[i] = (char *) mxCalloc(width+1, sizeof(char));
#ifdef DEBUG
mexPrintf("loop DynareMxArrayToString cNamesCharStr = %s \n", cNamesCharStr);
#endif
for (int i = 0; i < width; i++)
{
for (int j = 0; j < len; j++)
{
// Allow alphanumeric and underscores "_" only:
if (isalnum(cNamesCharStr[j+i*len]) || ('_' == cNamesCharStr[j+i*len]))
{
cNamesMX[j][i] = cNamesCharStr[j+i*len];
}
else cNamesMX[j][i] = '\0';
}
}
const char **ret = (const char **) mxCalloc(len, sizeof(char *));
for (int j = 0; j < len; j++)
{
cNamesMX[j][width] = '\0';
#ifdef DEBUG
// mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
#endif
char *token = (char *) mxCalloc(strlen(cNamesMX[j])+1, sizeof(char));
strcpy(token, cNamesMX[j]);
ret[j] = token;
#ifdef DEBUG
mexPrintf("ret [%d]= %s \n", j, ret[j]);
#endif
}
mxFree(cNamesMX);
return ret;
}
/*********************************** /***********************************
* Members of DynamicModelDLL for handling loading and calling * Members of DynamicModelDLL for handling loading and calling
* <model>_dynamic () function * <model>_dynamic () function
@ -86,7 +32,12 @@ DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const
: length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog) : length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog)
{ {
char fName[MAX_MODEL_NAME]; char fName[MAX_MODEL_NAME];
#if defined _WIN32
strcpy(fName, modName); strcpy(fName, modName);
#else
strcpy(fName, "./");
strcat(fName, modName);
#endif
using namespace std; using namespace std;
strcat(fName, "_dynamic"); strcat(fName, "_dynamic");
#ifdef DEBUG #ifdef DEBUG
@ -115,14 +66,14 @@ DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const
if ((dynamicHinstance == NULL) || dlerror()) if ((dynamicHinstance == NULL) || dlerror())
{ {
cerr << dlerror() << endl; cerr << dlerror() << endl;
mexPrintf("MexPrintf:Error loading DLL"); // mexPrintf("MexPrintf:Error loading DLL");
throw 1; throw 1;
} }
Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic"); Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic");
if ((Dynamic == NULL) || dlerror()) if ((Dynamic == NULL) || dlerror())
{ {
cerr << dlerror() << endl; cerr << dlerror() << endl;
mexPrintf("MexPrintf:Error finding DLL function"); // mexPrintf("MexPrintf:Error finding DLL function");
throw 2; throw 2;
} }
#endif #endif
@ -130,15 +81,15 @@ DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const
} }
catch (int i) catch (int i)
{ {
mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i); // mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i);
mexErrMsgTxt("Err: An error in Load and run DLL .\n"); // mexErrMsgTxt("Err: An error in Load and run DLL .\n");
return; return;
} }
catch (...) catch (...)
{ {
mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName); // mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName);
mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n"); // mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n");
return; return;
} }
} }
@ -163,15 +114,16 @@ DynamicModelDLL::close()
void void
DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams, DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams,
int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2) int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3)
{ {
double *dresidual, *dg1 = NULL, *dg2 = NULL; double *dresidual, *dg1 = NULL, *dg2 = NULL, *dg3 = NULL;
//int length=y.length(); // not! //int length=y.length(); // not!
//mexPrintf("y.length() = %d\n",y.length());
if ((jcols-nExog) != y.length()) if ((jcols-nExog) != y.length())
{ {
// throw DLL Error // throw DLL Error
mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n"); // mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n");
return; return;
} }
if (residual.length() < length) // dummy or insufficient if (residual.length() < length) // dummy or insufficient
@ -196,6 +148,11 @@ DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modPa
dg2 = const_cast<double *>(g2->base()); dg2 = const_cast<double *>(g2->base());
} }
dresidual = const_cast<double *>(residual.base()); dresidual = const_cast<double *>(residual.base());
if (g3 != NULL)
{
dg3 = const_cast<double *>(g3->base());
}
dresidual = const_cast<double *>(residual.base());
double *dy = const_cast<double *>(y.base()); double *dy = const_cast<double *>(y.base());
double *dx = const_cast<double *>(x.base()); double *dx = const_cast<double *>(x.base());
double *dbParams = const_cast<double *>(modParams->base()); double *dbParams = const_cast<double *>(modParams->base());
@ -216,25 +173,25 @@ DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modPa
#endif #endif
try try
{ {
Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2); Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3);
} }
catch (...) catch (...)
{ {
mexPrintf("MexPrintf: error in run Dynamic DLL \n"); // mexPrintf("MexPrintf: error in run Dynamic DLL \n");
} }
}; };
void void
DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams, DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2) Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3)
{ {
eval(y, x, modParams, nMax_lag, residual, g1, g2); eval(y, x, modParams, nMax_lag, residual, g1, g2, g3);
}; };
void void
DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams, DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2) Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3)
{ {
/** ignore given exogens and create new 2D x matrix since /** ignore given exogens and create new 2D x matrix since
@ -244,7 +201,7 @@ DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams,
TwoDMatrix &mx = *(new TwoDMatrix(nMax_lag+1, nExog)); TwoDMatrix &mx = *(new TwoDMatrix(nMax_lag+1, nExog));
mx.zeros(); // initialise shocks to 0s mx.zeros(); // initialise shocks to 0s
eval(y, mx, modParams, nMax_lag, residual, g1, g2); eval(y, mx, modParams, nMax_lag, residual, g1, g2, g3);
}; };

15
mex/sources/k_order_perturbation/dynamic_dll.h
View File

@ -43,7 +43,7 @@ typedef void *(DynamicFn)
typedef void (*DynamicFn) typedef void (*DynamicFn)
#endif #endif
(double *y, double *x, int nb_row_x, double *params, (double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2); int it_, double *residual, double *g1, double *g2, double *g3);
//DynamicFn Dynamic; //DynamicFn Dynamic;
@ -86,18 +86,17 @@ public:
// evaluate Dynamic model DLL // evaluate Dynamic model DLL
void void
eval(double *y, double *x, int nb_row_x, double *params, eval(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2) int it_, double *residual, double *g1, double *g2, double *g3)
{ {
Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2); Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2, g3);
}; };
void eval(const Vector &y, const Vector &x, const Vector *params, void eval(const Vector &y, const Vector &x, const Vector *params,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2); Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3);
void eval(const Vector &y, const TwoDMatrix &x, const Vector *params, void eval(const Vector &y, const TwoDMatrix &x, const Vector *params,
int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2); int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3);
void eval(const Vector &y, const TwoDMatrix &x, const Vector *params, void eval(const Vector &y, const TwoDMatrix &x, const Vector *params,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2); Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3);
// 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 DLL: If the referenced object was successfully closed,
// close() returns 0, non 0 otherwise // close() returns 0, non 0 otherwise
int close(); int close();

106
mex/sources/k_order_perturbation/k_ord_dynare.cpp
View File

@ -84,7 +84,7 @@ KordpDynare::KordpDynare(const char **endo, int num_endo,
} }
catch (...) catch (...)
{ {
mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n"); // mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n");
throw DynareException(__FILE__, __LINE__, string("Could not construct Name Lists. \n")); throw DynareException(__FILE__, __LINE__, string("Could not construct Name Lists. \n"));
} }
} }
@ -159,7 +159,7 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yy, const Vector &xx)
{ {
dynamicDLL.eval(yy, xx, //int nb_row_x, dynamicDLL.eval(yy, xx, //int nb_row_x,
params, //int it_, params, //int it_,
out, NULL, NULL); out, NULL, NULL, NULL);
} }
@ -176,7 +176,7 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
#endif #endif
dynamicDLL.eval(yy, xx, //int nb_row_x, dynamicDLL.eval(yy, xx, //int nb_row_x,
params, //int it_, params, //int it_,
out, NULL, NULL); out, NULL, NULL, NULL);
} }
/************************************************ /************************************************
* this is main derivative calculation functin that indirectly calls dynamic.dll * this is main derivative calculation functin that indirectly calls dynamic.dll
@ -188,13 +188,14 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx)
{ {
// Hessian TwoDMatrix *g2; // Hessian TwoDMatrix *g2;
TwoDMatrix *g2 = NULL; TwoDMatrix *g2 = NULL;
TwoDMatrix *g3 = NULL;
//Jacobian //Jacobian
TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian
g1->zeros(); g1->zeros();
if ((nJcols != g1->ncols()) && (nY != g1->nrows())) if ((nJcols != g1->ncols()) && (nY != g1->nrows()))
{ {
mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: Created wrong jacobian"); // mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: Created wrong jacobian");
return; return;
} }
@ -210,6 +211,11 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx)
//g2->print(); //g2->print();
#endif #endif
} }
if (nOrder > 2)
{
g3 = new TwoDMatrix((int) (*NNZD)[2],3);
g3->zeros();
}
Vector &out = *(new Vector(nY)); Vector &out = *(new Vector(nY));
out.zeros(); out.zeros();
const Vector *llxYYp; // getting around the constantness const Vector *llxYYp; // getting around the constantness
@ -228,34 +234,25 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx)
#endif #endif
try try
{ {
dynamicDLL.eval(llxYY, xx, //int nb_row_x, dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3);
params, //int it_,
out, g1, g2);
// }
} }
catch (...) catch (...)
{ {
mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives"); // mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives");
return; return;
} }
if ((nJcols != g1->ncols()) && (nY != g1->nrows())) if ((nJcols != g1->ncols()) && (nY != g1->nrows()))
{ {
mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian"); // mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian");
return; return;
} }
// ReorderCols(g1, JacobianIndices); and populate container
populateDerivativesContainer(g1, 1, JacobianIndices); populateDerivativesContainer(g1, 1, JacobianIndices);
if (nOrder > 1) if (nOrder > 1)
{
// ReorderBlocks(g2,JacobianIndices);
#ifdef DEBUG
mexPrintf(" post dll g2 cols %d rows %d \n", g2->numCols(), g2->numRows());
for (int ii=0;ii<g2->numRows(); ii++)
mexPrintf(" g2[%d]: %d %d %f \n", ii, (int)g2->get(ii,0),(int)g2->get(ii,1),g2->get(ii,2));
//g2->print();
#endif
populateDerivativesContainer(g2, 2, JacobianIndices); populateDerivativesContainer(g2, 2, JacobianIndices);
} if (nOrder > 2)
populateDerivativesContainer(g3, 3, JacobianIndices);
} }
/* This version is not currently in use */ /* This version is not currently in use */
@ -271,7 +268,7 @@ KordpDynare::calcDerivatives(const Vector &yy, ogu::Jacobian &jacob)
xx.zeros(); xx.zeros();
dynamicDLL.eval(yy, xx, //int nb_row_x, dynamicDLL.eval(yy, xx, //int nb_row_x,
params, //int it_, params, //int it_,
out, jj, NULL); out, jj, NULL, NULL);
// model derivatives FSSparseTensor instance // model derivatives FSSparseTensor instance
FSSparseTensor &mdTi = *(new FSSparseTensor(1, jj->ncols(), jj->nrows())); FSSparseTensor &mdTi = *(new FSSparseTensor(1, jj->ncols(), jj->nrows()));
for (int i = 0; i < jj->ncols(); i++) for (int i = 0; i < jj->ncols(); i++)
@ -309,7 +306,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
#endif #endif
// model derivatives FSSparseTensor instance // model derivatives FSSparseTensor instance
FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, g->nrows())); FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY));
IntSequence s(ord, 0); IntSequence s(ord, 0);
//s[0] = 0; //s[0] = 0;
@ -364,14 +361,47 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
} }
} }
} }
else if (ord == 3)
{
int nJcols1 = nJcols-nExog;
int nJcols2 = nJcols*nJcols;
vector<int> revOrder(nJcols1);
for (int i = 0; i < nJcols1; i++)
revOrder[(*vOrder)[i]] = i;
for (int i = 0; i < g->nrows(); i++)
{
int j = (int)g->get(i,0)-1;
int i1 = (int)g->get(i,1) -1;
int s0 = (int)floor(i1/nJcols2);
int i2 = i1 - nJcols2*s0;
int s1 = (int)floor(i2/nJcols);
int s2 = i2 - nJcols*s1;
if (s0 < nJcols1)
s[0] = revOrder[s0];
else
s[0] = s0;
if (s1 < nJcols1)
s[1] = revOrder[s1];
else
s[1] = s1;
if (s2 < nJcols1)
s[2] = revOrder[s2];
else
s[2] = s2;
if ((s[2] >= s[1]) && (s[1] >= s[0]))
{
double x = g->get(i,2);
mdTi->insert(s, j, x);
}
}
}
#ifdef DEBUG #ifdef DEBUG
mexPrintf("k_ord_dynare.cpp: END populate FSSparseTensor in calcDerivatives: ord=%d \n",ord); mexPrintf("k_ord_dynare.cpp: END populate FSSparseTensor in calcDerivatives: ord=%d \n",ord);
mdTi->print();
#endif #endif
// md container // md container
//md.clear();// this is to be used only for 1st order!! //md.clear();// this is to be used only for 1st order!!
md.remove(Symmetry(ord)); md.remove(Symmetry(ord));
md.insert(mdTi); //(&mdTi); md.insert(mdTi);
} }
void void
@ -406,7 +436,7 @@ KordpDynare::LLxSteady(const Vector &yS)
{ {
if ((nJcols-nExog) == yS.length()) if ((nJcols-nExog) == yS.length())
{ {
mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size"); // mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size");
return NULL; return NULL;
} }
// create temporary square 2D matrix size nEndo x nEndo (sparse) // create temporary square 2D matrix size nEndo x nEndo (sparse)
@ -426,13 +456,13 @@ KordpDynare::LLxSteady(const Vector &yS)
} }
catch (const TLException &e) catch (const TLException &e)
{ {
mexPrintf("Caugth TL exception in LLxSteady: "); // mexPrintf("Caugth TL exception in LLxSteady: ");
e.print(); e.print();
return NULL; // 255; return NULL; // 255;
} }
catch (...) catch (...)
{ {
mexPrintf(" Error in LLxSteady - wrong index?"); // mexPrintf(" Error in LLxSteady - wrong index?");
} }
#ifdef DEBUG #ifdef DEBUG
@ -514,13 +544,13 @@ KordpDynare::ReorderDynareJacobianIndices(const vector<int> *varOrder)
} }
catch (const TLException &e) catch (const TLException &e)
{ {
mexPrintf("Caugth TL exception in ReorderIndices: "); // mexPrintf("Caugth TL exception in ReorderIndices: ");
e.print(); e.print();
return NULL; // 255; return NULL; // 255;
} }
catch (...) catch (...)
{ {
mexPrintf(" Error in ReorderIndices - wrong index?"); // mexPrintf(" Error in ReorderIndices - wrong index?");
} }
//add the indices for the nExog exogenous jacobians //add the indices for the nExog exogenous jacobians
for (j = nJcols-nExog; j < nJcols; j++) for (j = nJcols-nExog; j < nJcols; j++)
@ -547,7 +577,7 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector<int> *vOrder)
if (tdx->ncols() > vOrder->size()) if (tdx->ncols() > vOrder->size())
{ {
mexPrintf(" Error in ReorderColumns - size of order var is too small"); // mexPrintf(" Error in ReorderColumns - size of order var is too small");
return; return;
} }
TwoDMatrix tmp(*tdx); // temporary 2D matrix TwoDMatrix tmp(*tdx); // temporary 2D matrix
@ -561,13 +591,13 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector<int> *vOrder)
} }
catch (const TLException &e) catch (const TLException &e)
{ {
printf("Caugth TL exception in ReorderColumns: "); //printf("Caugth TL exception in ReorderColumns: ");
e.print(); e.print();
return; // 255; return; // 255;
} }
catch (...) catch (...)
{ {
mexPrintf(" Error in ReorderColumns - wrong index?"); // mexPrintf(" Error in ReorderColumns - wrong index?");
} }
} }
void void
@ -585,13 +615,13 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const int *vOrder)
} }
catch (const TLException &e) catch (const TLException &e)
{ {
printf("Caugth TL exception in ReorderColumns: "); //printf("Caugth TL exception in ReorderColumns: ");
e.print(); e.print();
return; // 255; return; // 255;
} }
catch (...) catch (...)
{ {
mexPrintf(" Error in ReorderColumns - wrong index?"); // mexPrintf(" Error in ReorderColumns - wrong index?");
} }
} }
@ -610,7 +640,7 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector<int> *vOrder)
int ibOrder = (int) dbOrder; int ibOrder = (int) dbOrder;
if ((double) ibOrder != dbOrder || ibOrder > nOrder) if ((double) ibOrder != dbOrder || ibOrder > nOrder)
{ {
mexPrintf(" Error in ReorderBlocks - wrong order %d", dbOrder); // mexPrintf(" Error in ReorderBlocks - wrong order %d", dbOrder);
return; return;
} }
TwoDMatrix tmp(*tdx); // temporary 2D matrix TwoDMatrix tmp(*tdx); // temporary 2D matrix
@ -633,7 +663,7 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector<int> *vOrder)
//ReorderColumns(TwoDMatrix * subtdx, const vector<int> * vOrder) //ReorderColumns(TwoDMatrix * subtdx, const vector<int> * vOrder)
if (tdx->ncols() > vOrder->size()) if (tdx->ncols() > vOrder->size())
{ {
mexPrintf(" Error in ReorderColumns - size of order var is too small"); // mexPrintf(" Error in ReorderColumns - size of order var is too small");
return; return;
} }
// reorder the columns // reorder the columns
@ -644,13 +674,13 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector<int> *vOrder)
} }
catch (const TLException &e) catch (const TLException &e)
{ {
printf("Caugth TL exception in ReorderColumns: "); // printf("Caugth TL exception in ReorderColumns: ");
e.print(); e.print();
return; // 255; return; // 255;
} }
catch (...) catch (...)
{ {
mexPrintf(" Error in ReorderColumns - wrong index?"); // mexPrintf(" Error in ReorderColumns - wrong index?");
} }
} }
} }

105
mex/sources/k_order_perturbation/k_order_perturbation.cpp
View File

@ -86,10 +86,10 @@ extern "C" {
mexFunction(int nlhs, mxArray *plhs[], mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]) int nrhs, const mxArray *prhs[])
{ {
if (nrhs < 5) // if (nrhs < 5)
mexErrMsgTxt("Must have at least 5 input parameters.\n"); // mexErrMsgTxt("Must have at least 5 input parameters.\n");
if (nlhs == 0) //if (nlhs == 0)
mexErrMsgTxt("Must have at least 1 output parameter.\n"); // mexErrMsgTxt("Must have at least 1 output parameter.\n");
const mxArray *dr = prhs[0]; const mxArray *dr = prhs[0];
const int check_flag = (int) mxGetScalar(prhs[1]); const int check_flag = (int) mxGetScalar(prhs[1]);
@ -119,7 +119,14 @@ extern "C" {
kOrder = (int) mxGetScalar(mxFldp); kOrder = (int) mxGetScalar(mxFldp);
else else
kOrder = 1; kOrder = 1;
// if (kOrder == 1 && nlhs != 1)
// mexErrMsgTxt("k_order_perturbation at order 1 requires exactly 1 argument in output\'n");
// else if (kOrder > 1 && nlhs != kOrder+1)
// mexErrMsgTxt("k_order_perturbation at order > 1 requires exactly order + 1 argument in output\'n");
double qz_criterium = 1+1e-6; double qz_criterium = 1+1e-6;
mxFldp = mxGetField(options_, 0, "qz_criterium"); mxFldp = mxGetField(options_, 0, "qz_criterium");
if (mxIsNumeric(mxFldp)) if (mxIsNumeric(mxFldp))
@ -146,9 +153,11 @@ extern "C" {
// mxFldp = mxGetField(oo_, 0,"steady_state" ); // use in order of declaration // mxFldp = mxGetField(oo_, 0,"steady_state" ); // use in order of declaration
mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var
//mexPrintf("mxFldp %x",mxFldp);
// mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // and NOT extended ys // mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // and NOT extended ys
dparams = (double *) mxGetData(mxFldp); dparams = (double *) mxGetData(mxFldp);
const int nSteady = (int) mxGetM(mxFldp); const int nSteady = (int) mxGetM(mxFldp);
// mexPrintf("nSteady %d\n",nSteady);
Vector *ySteady = new Vector(dparams, nSteady); Vector *ySteady = new Vector(dparams, nSteady);
mxFldp = mxGetField(dr, 0, "nstatic"); mxFldp = mxGetField(dr, 0, "nstatic");
@ -198,6 +207,8 @@ extern "C" {
dparams = (double *) mxGetData(mxFldp); dparams = (double *) mxGetData(mxFldp);
npar = (int) mxGetN(mxFldp); npar = (int) mxGetN(mxFldp);
int nrows = (int) mxGetM(mxFldp); int nrows = (int) mxGetM(mxFldp);
#ifdef DEBUG #ifdef DEBUG
mexPrintf("ll_Incidence nrows=%d, ncols = %d .\n", nrows, npar); mexPrintf("ll_Incidence nrows=%d, ncols = %d .\n", nrows, npar);
#endif #endif
@ -227,10 +238,9 @@ extern "C" {
#endif #endif
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
mxFldp = mxGetField(M_, 0, "var_order_endo_names"); mxFldp = mxGetField(M_, 0, "var_order_endo_names");
mexPrintf("k_order_perturbation: Get nendo .\n"); // mexPrintf("k_order_perturbation: Get nendo .\n");
const int nendo = (int) mxGetM(mxFldp); const int nendo = (int) mxGetM(mxFldp);
const int widthEndo = (int) mxGetN(mxFldp); const int widthEndo = (int) mxGetN(mxFldp);
const char **endoNamesMX = DynareMxArrayToString(mxFldp, nendo, widthEndo); const char **endoNamesMX = DynareMxArrayToString(mxFldp, nendo, widthEndo);
@ -301,6 +311,7 @@ extern "C" {
#ifdef DEBUG #ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynare constructor .\n"); mexPrintf("k_order_perturbation: Calling dynare constructor .\n");
#endif #endif
// make KordpDynare object // make KordpDynare object
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames, KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames,
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
@ -311,6 +322,7 @@ extern "C" {
#ifdef DEBUG #ifdef DEBUG
mexPrintf("k_order_perturbation: Call Approximation constructor with qz_criterium=%f \n", qz_criterium); mexPrintf("k_order_perturbation: Call Approximation constructor with qz_criterium=%f \n", qz_criterium);
#endif #endif
Approximation app(dynare, journal, nSteps, false, qz_criterium); Approximation app(dynare, journal, nSteps, false, qz_criterium);
// run stochastic steady // run stochastic steady
#ifdef DEBUG #ifdef DEBUG
@ -367,29 +379,31 @@ extern "C" {
int nb_row_x; int nb_row_x;
ysteady = NULL; ysteady = NULL;
if (nlhs >= 1) if (kOrder == 1)
{ {
/* Set the output pointer to the output matrix ysteady. */ /* Set the output pointer to the output matrix ysteady. */
plhs[0] = mxCreateDoubleMatrix(nEndo, 1, mxREAL); map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
/* Create a C pointer to a copy of the output ysteady. */ ++cit;
TwoDMatrix tmp_ss(nEndo, 1, mxGetPr(plhs[0])); plhs[0] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
tmp_ss = (const TwoDMatrix &)ss; TwoDMatrix g((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[0]));
g = (const TwoDMatrix &)(*cit).second;
// /* Create a C pointer to a copy of the output ysteady. */
// TwoDMatrix tmp_ss(nEndo, 1, mxGetPr(plhs[0]));
// tmp_ss = (const TwoDMatrix &)ss;
#ifdef DEBUG #ifdef DEBUG
// tmp_ss.print(); !! This print Crashes??? // tmp_ss.print(); !! This print Crashes???
#endif #endif
} }
if (nlhs >= 2) if (kOrder >= 2)
{ {
/* Set the output pointer to the combined output matrix gyu = [gy gu]. */ int ii = 0;
int ii = 1;
for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin(); for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
((cit != mm.end()) && (ii < nlhs)); ++cit) ((cit != mm.end()) && (ii < nlhs)); ++cit)
{ {
//if ((*cit).first!="g_0")
{ {
plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL); plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
TwoDMatrix dgyu((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[ii])); TwoDMatrix g_ii((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[ii]));
dgyu = (const TwoDMatrix &)(*cit).second; g_ii = (const TwoDMatrix &)(*cit).second;
#ifdef DEBUG #ifdef DEBUG
mexPrintf("k_order_perturbation: cit %d print: \n", ii); mexPrintf("k_order_perturbation: cit %d print: \n", ii);
(*cit).second.print(); (*cit).second.print();
@ -462,9 +476,64 @@ extern "C" {
e.print(errfd); e.print(errfd);
fclose(errfd); fclose(errfd);
return; // 255; return; // 255;
} //catch } //catch
}; // end of mexFunction() }; // end of mexFunction()
}; // end of extern C }; // end of extern C
//////////////////////////////////////////////////////
// Convert Matlab Dynare endo and exo names array to C type array of string pointers
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
// hence a rather low level approach is needed
///////////////////////////////////////////////////////
const char **
DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width)
{
char *cNamesCharStr = mxArrayToString(mxFldp);
const char **ret = DynareMxArrayToString(cNamesCharStr, len, width);
return ret;
}
const char **
DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width)
{
char **cNamesMX;
cNamesMX = (char **) calloc(len, sizeof(char *));
for(int i = 0; i < len; i++)
cNamesMX[i] = (char *) calloc(width+1, sizeof(char));
#ifdef DEBUG
mexPrintf("loop DynareMxArrayToString cNamesCharStr = %s \n", cNamesCharStr);
#endif
for (int i = 0; i < width; i++)
{
for (int j = 0; j < len; j++)
{
// Allow alphanumeric and underscores "_" only:
if (isalnum(cNamesCharStr[j+i*len]) || ('_' == cNamesCharStr[j+i*len]))
{
cNamesMX[j][i] = cNamesCharStr[j+i*len];
}
else cNamesMX[j][i] = '\0';
}
}
const char **ret = (const char **) calloc(len, sizeof(char *));
for (int j = 0; j < len; j++)
{
cNamesMX[j][width] = '\0';
#ifdef DEBUG
// mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
#endif
char *token = (char *) calloc(strlen(cNamesMX[j])+1, sizeof(char));
strcpy(token, cNamesMX[j]);
ret[j] = token;
#ifdef DEBUG
mexPrintf("ret [%d]= %s \n", j, ret[j]);
#endif
}
mxFree(cNamesMX);
return ret;
}
#endif // ifdef MATLAB_MEX_FILE to exclude mexFunction for other applications #endif // ifdef MATLAB_MEX_FILE to exclude mexFunction for other applications