Updates for prototype based on reading output mat file.

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2388 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
george 2009-02-02 12:57:07 +00:00
parent 60feef4a0a
commit c98d175e4b
4 changed files with 99 additions and 44 deletions

View File

@ -86,12 +86,13 @@ KordpDynare::KordpDynare(const char** endo, int num_endo,
const char** exo, int nexog, int npar, //const char** par,
Vector* ysteady, TwoDMatrix* vcov, Vector* inParams, int nstat,
int npred, int nforw, int nboth, const int jcols, const int nsteps, int norder, //const char* modName,
Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const int* var_order,
const TwoDMatrix * llincidence )
Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const vector<int>* var_order,
const TwoDMatrix * llincidence, double criterium )
: nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar),
nYs(npred + nboth), nYss(nboth + nforw),nY(num_endo), nJcols(jcols), nSteps(nsteps), nOrder(norder),
journal(jr), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params (inParams),
md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order), ll_Incidence(llincidence)
md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order),
ll_Incidence(llincidence), qz_criterium(criterium)
{
#ifdef DEBUG
mexPrintf("k_ord_dynare Dynare constructor: ny=%d, order=%d, nPar=%d .\n", nY,nOrder,nPar);
@ -138,7 +139,7 @@ KordpDynare::KordpDynare(const KordpDynare& dynare)
ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md),
dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol),
varOrder(dynare.varOrder), ll_Incidence(dynare.ll_Incidence),
JacobianIndices(dynare.JacobianIndices)
JacobianIndices(dynare.JacobianIndices), qz_criterium(dynare.qz_criterium)
{
/// model = dynare.model->clone();
ySteady = new Vector(*(dynare.ySteady));
@ -391,7 +392,7 @@ Vector * KordpDynare::LLxSteady( const Vector& yS){
// populate (non-sparse) vector with ysteady values
for (int i=0;i<nY;i++){
if(ll_Incidence->get(ll_row,i))
(*llxSteady)[ll_Incidence->get(ll_row,i)-1] = yS[i];
(*llxSteady)[((int)ll_Incidence->get(ll_row,i))-1] = yS[i];
}
}
} catch (const TLException& e) {
@ -435,16 +436,20 @@ Vector * KordpDynare::LLxSteady( const Vector& yS){
exogen
************************************/
int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){
vector<int> * KordpDynare::ReorderDynareJacobianIndices( const vector<int> * varOrder){
// if ((nJcols != tdx->ncols()) && ( nY != tdx->nrows())) {
// mexPrintf("k_ord_dynare.cpp: Error in ReorderBlocks: wrong size of jacobian");
// return;
// }
// create temporary square 2D matrix size nEndo x nEndo (sparse)
// for the lag, current and lead blocks of the jacobian
int * JacobianIndices = (int*) calloc(nJcols+1, sizeof(int));
// int * JacobianIndices = (int*) calloc(nJcols+1, sizeof(int));
vector<int> * JacobianIndices = new vector<int>(nJcols);
vector <int> tmp(nY);
int i,j, rjoff=nJcols-nExog-1; //, ll_off, j;
#ifdef DEBUG
mexPrintf("ReorderDynareJacobianIndice:ll_Incidence->nrows() =%d .\n", ll_Incidence->nrows());
#endif
try{
for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++)
{
@ -452,13 +457,17 @@ int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){
// the lag, current and lead blocks of the jacobian respectively
for (i=0;i<nY;i++){
// tmp[varOrder[j]]=(int)ll_Incidence->get(ll_row,j);
tmp[i]=((int)ll_Incidence->get(ll_row,varOrder[i]-1));
tmp[i]=((int)ll_Incidence->get(ll_row,(*varOrder)[i]-1));
#ifdef DEBUG
mexPrintf("get(ll_row,(*varOrder)[%d]-1)) = tmp[%d]=%d .\n",
i, i, (int)ll_Incidence->get(ll_row,(*varOrder)[i]-1));
#endif
}
// write the reordered blocks back to the jacobian
// in reverse order
for (j=nY-1;j>=0;j--){
if (tmp[j]){
JacobianIndices[rjoff]=tmp[j] -1;
(*JacobianIndices)[rjoff]=tmp[j] -1;
rjoff--;
if (rjoff<0){
// mexPrintf(" Error in ReorderIndices - negative rjoff index?");
@ -477,11 +486,11 @@ int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){
}
//add the indices for the nExog exogenous jacobians
for (j=nJcols-nExog;j<nJcols;j++){
JacobianIndices[j]=j;
(*JacobianIndices)[j]=j;
}
#ifdef DEBUG
for (j=0;j<nJcols;j++)
mexPrintf("ReorderDynareJacobianIndice: [%d] =%d .\n", j, JacobianIndices[j]);
mexPrintf("ReorderDynareJacobianIndice: [%d] =%d .\n", j, (*JacobianIndices)[j]);
#endif
return JacobianIndices;
}
@ -494,7 +503,28 @@ int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){
* of any of its elements too.
************************************/
void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * varsOrder){
void KordpDynare::ReorderCols(TwoDMatrix * tdx, const vector<int> * vOrder){
if (tdx->ncols() > vOrder->size()){
mexPrintf(" Error in ReorderColumns - size of order var is too small");
return;
}
TwoDMatrix tmp(*tdx); // temporary 2D matrix
TwoDMatrix &tmpR=tmp;
tdx->zeros();// empty original matrix
// reorder the columns
try{
for (int i =0; i<tdx->ncols() ; i++)
tdx->copyColumn(tmpR,(*vOrder)[i],i);
} catch (const TLException& e) {
printf("Caugth TL exception in ReorderColumns: ");
e.print();
return;// 255;
}catch (...){
mexPrintf(" Error in ReorderColumns - wrong index?");
}
}
void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * vOrder){
TwoDMatrix tmp(*tdx); // temporary 2D matrix
TwoDMatrix &tmpR=tmp;
@ -502,7 +532,7 @@ void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * varsOrder){
// reorder the columns
try{
for (int i =0; i<tdx->ncols() ; i++)
tdx->copyColumn(tmpR,varsOrder[i],i);
tdx->copyColumn(tmpR,vOrder[i],i);
} catch (const TLException& e) {
printf("Caugth TL exception in ReorderColumns: ");
e.print();

View File

@ -141,16 +141,18 @@ friend class DynamicModelDLL;
DynareExogNameList* denl;
DynareStateNameList* dsnl;
const double ss_tol;
const int* varOrder;
const vector<int>* varOrder;
const TwoDMatrix * ll_Incidence;
int * JacobianIndices;
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 int* varOrder, const TwoDMatrix * ll_Incidence );
const vector<int>* varOrder, const TwoDMatrix * ll_Incidence,
double qz_criterium );
/** Makes a deep copy of the object. */
KordpDynare(const KordpDynare& dyn);
@ -212,11 +214,13 @@ public:
DynamicModel* clone() const
{return new KordpDynare(*this);}
void ReorderCols(TwoDMatrix * tdx, const int * varOrder);
void ReorderCols(TwoDMatrix * tdx, const vector<int> * varOrder);
Vector * KordpDynare::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);
};

View File

@ -110,6 +110,11 @@ extern "C" {
kOrder = (int)mxGetScalar(mxFldp);
else
kOrder = 1;
double qz_criterium = 1+1e-6;
mxFldp = mxGetField(options_, 0,"qz_criterium" );
if (mxIsNumeric(mxFldp))
qz_criterium = (int)mxGetScalar(mxFldp);
mxFldp = mxGetField(M_, 0,"params" );
double * dparams = (double *) mxGetData(mxFldp);
@ -117,30 +122,25 @@ extern "C" {
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
const mxArray* const mxParFldp = mxGetField(M_, 0,"params" );
Vector params_vec((const double*)mxGetPr(mxParFldp), npar);
//const mxArray* const mxParFldp = mxGetField(M_, 0,"params" );
//Vector params_vec((const double*)mxGetPr(mxParFldp), npar);
#ifdef DEBUG
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_vec[%d]= %g.\n", i, params_vec[i] ); }
#endif
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(oo_, 0,"steady_state" ); // use in order of declaration
// mxFldp = mxGetField(dr, 0,"ys" ); // and not in order of dr.order_var
mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // extended ys
// mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // and NOT extended ys
dparams = (double *) mxGetData(mxFldp);
const int nSteady = (int)mxGetM(mxFldp);
Vector * ySteady = new Vector(dparams, nSteady);
@ -172,18 +172,39 @@ extern "C" {
nPred -= nBoth; // correct nPred for nBoth.
mxFldp = mxGetField(dr, 0,"order_var" );
int * var_order = (int *) mxGetData(mxFldp);
dparams = (double *) mxGetData(mxFldp);
npar = (int)mxGetM(mxFldp);
if (npar != nEndo) { //(nPar != npar)
mexErrMsgTxt("Incorrect number of input var_order vars.\n");
//return;
}
// 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] =(int)(*(dparams++));//[v];
#ifdef DEBUG
mexPrintf("var_order_vp)[%d] = %d .\n", v, (*var_order_vp)[v] );
#endif
}
// the lag, current and lead blocks of the jacobian respectively
mxFldp = mxGetField(M_, 0,"lead_lag_incidence" );
dparams = (double *) mxGetData(mxFldp);
npar = (int)mxGetN(mxFldp);
TwoDMatrix * llincidence = new TwoDMatrix(npar, nEndo, dparams);
int nrows = (int)mxGetM(mxFldp);
#ifdef DEBUG
mexPrintf("ll_Incidence nrows=%d, ncols = %d .\n", nrows, npar);
#endif
TwoDMatrix * llincidence = new TwoDMatrix ( nrows, npar,dparams);
if (npar != nEndo){
mexPrintf("Incorrect length of lead lag incidences: ncol=%d != nEndo =%d .\n", npar, nEndo);
return;
}
#ifdef DEBUG
for (int j=0;j<nrows;j++){
for (int i=0;i<nEndo;i++){
mexPrintf("llincidence->get(%d,%d) =%d .\n",
j, i, (int)llincidence->get(j,i));}}
#endif
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
@ -218,16 +239,14 @@ extern "C" {
************/
if ((nEndo != nendo)||(nExog != nexo)) { //(nPar != npar)
mexErrMsgTxt("Incorrect number of input parameters.\n");
//return;
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 < 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]); }
@ -257,8 +276,8 @@ extern "C" {
// 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,
llincidence );
jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp,
llincidence,qz_criterium );
/************
// make list of shocks for which we will do IRFs
vector<int> irf_list_ind;
@ -650,8 +669,8 @@ void DynamicModelDLL::eval(const Vector&y, const TwoDMatrix&x, const Vector* mo
, g1->ncols(),g1->nrows());
for (int i = 0; i < modParams->length(); i++) {
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); }
for (int i = 0; i < length; i++) {
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, y[i]);}
for (int i = 0; i < jcols-nExog; i++) {
mexPrintf("k_ord_perturbation: Ys[%d]= %g.\n", i, y[i]);}
mexPrintf("k_order_perturbation: call <model> Dynamic dParams= %d , , dy = %d dx = %d .\n"
,dbParams[0],dy[0],dx[0]);

View File

@ -31,7 +31,7 @@
int main(int argc, char* argv[])
{
double qz_criterium = 1+1e-6;
const int check_flag = 0;
const char* fName = "fs2000k";//mxArrayToString(mFname);
@ -48,7 +48,6 @@ int main(int argc, char* argv[])
0.7870,
0.0200
};
Vector * modParams = new Vector(dparams, npar);
#ifdef DEBUG
@ -106,7 +105,10 @@ int main(int argc, char* argv[])
// 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
@ -195,8 +197,8 @@ int main(int argc, char* argv[])
// 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
, llincidence );
, jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp //var_order
, llincidence, qz_criterium);
// intiate tensor library
#ifdef DEBUG
mexPrintf("k_order_perturbation: Call tls init\n");