2009-05-06 12:10:27 +02:00
/*
2010-03-09 15:20:18 +01:00
* Copyright ( C ) 2008 - 2010 Dynare Team
2009-05-19 10:57:07 +02:00
*
* 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/>.
*/
2009-05-06 12:10:27 +02:00
2010-08-30 17:11:58 +02:00
/*
Defines the entry point for the k - order perturbation application DLL .
Inputs :
1 ) dr
2 ) M_
3 ) options
4 ) oo_
5 ) string containing the MEX extension ( with a dot at the beginning )
Outputs :
- if order = = 1 : only g_1
- if order = = 2 : g_0 , g_1 , g_2
- if order = = 3 : g_0 , g_1 , g_2 , g_3
*/
2009-05-06 12:10:27 +02:00
2009-12-17 12:00:50 +01:00
# include "k_ord_dynare.hh"
# include "dynamic_dll.hh"
2009-05-06 12:10:27 +02:00
2009-11-30 17:31:27 +01:00
# include <cmath>
# include <cstring>
2009-05-06 12:10:27 +02:00
# include <cctype>
2010-02-08 16:58:24 +01:00
# include <cassert>
2009-05-06 12:10:27 +02:00
2009-10-02 12:23:49 +02:00
# if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) // exclude mexFunction for other applications
2009-11-30 17:31:27 +01:00
//////////////////////////////////////////////////////
2010-03-09 18:22:33 +01:00
// Convert MATLAB Dynare endo and exo names array to a vector<string> array of string pointers
2009-11-30 17:31:27 +01:00
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
// hence a rather low level approach is needed
///////////////////////////////////////////////////////
2010-03-09 18:22:33 +01:00
void
DynareMxArrayToString ( const mxArray * mxFldp , const int len , const int width , vector < string > & out )
2009-11-30 17:31:27 +01:00
{
char * cNamesCharStr = mxArrayToString ( mxFldp ) ;
2009-12-16 18:18:38 +01:00
2010-03-09 18:22:33 +01:00
out . resize ( len ) ;
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 ] = = ' _ ' ) )
out [ j ] + = cNamesCharStr [ j + i * len ] ;
2009-11-30 17:31:27 +01:00
}
2009-05-06 12:10:27 +02:00
extern " C " {
2009-05-27 16:28:23 +02:00
void
mexFunction ( int nlhs , mxArray * plhs [ ] ,
int nrhs , const mxArray * prhs [ ] )
2009-05-06 12:10:27 +02:00
{
2009-11-30 17:31:27 +01:00
if ( nrhs < 5 )
2010-08-30 17:11:58 +02:00
mexErrMsgTxt ( " Must have exactly 5 input parameters. " ) ;
2009-12-16 18:18:38 +01:00
2009-05-27 16:28:23 +02:00
const mxArray * dr = prhs [ 0 ] ;
2010-08-30 17:11:58 +02:00
const mxArray * M_ = prhs [ 1 ] ;
const mxArray * options_ = prhs [ 2 ] ;
const mxArray * oo_ = prhs [ 3 ] ;
2009-05-27 16:28:23 +02:00
mxArray * mFname = mxGetField ( M_ , 0 , " fname " ) ;
if ( ! mxIsChar ( mFname ) )
2009-12-16 18:18:38 +01:00
mexErrMsgTxt ( " Input must be of type char. " ) ;
2009-11-30 17:31:27 +01:00
string fName = mxArrayToString ( mFname ) ;
2010-08-30 17:11:58 +02:00
const mxArray * mexExt = prhs [ 4 ] ;
string dfExt = mxArrayToString ( mexExt ) ; // Dynamic file extension, e.g. ".dll" or ".mexw32"
2009-05-27 16:28:23 +02:00
2009-05-06 12:10:27 +02:00
int kOrder ;
2009-05-27 16:28:23 +02:00
mxArray * mxFldp = mxGetField ( options_ , 0 , " order " ) ;
2009-05-06 12:10:27 +02:00
if ( mxIsNumeric ( mxFldp ) )
2009-05-27 16:28:23 +02:00
kOrder = ( int ) mxGetScalar ( mxFldp ) ;
2009-05-06 12:10:27 +02:00
else
kOrder = 1 ;
2009-12-16 18:18:38 +01:00
2009-11-30 17:31:27 +01:00
if ( kOrder = = 1 & & nlhs ! = 1 )
mexErrMsgTxt ( " k_order_perturbation at order 1 requires exactly 1 argument in output " ) ;
else if ( kOrder > 1 & & nlhs ! = kOrder + 1 )
2010-08-30 17:11:58 +02:00
mexErrMsgTxt ( " k_order_perturbation at order > 1 requires exactly order+1 arguments in output " ) ;
2009-12-16 18:18:38 +01:00
2009-05-06 12:10:27 +02:00
double qz_criterium = 1 + 1e-6 ;
2009-05-27 16:28:23 +02:00
mxFldp = mxGetField ( options_ , 0 , " qz_criterium " ) ;
2009-05-06 12:10:27 +02:00
if ( mxIsNumeric ( mxFldp ) )
2009-05-27 16:28:23 +02:00
qz_criterium = ( double ) mxGetScalar ( mxFldp ) ;
2010-08-30 17:11:58 +02:00
mxFldp = mxGetField ( M_ , 0 , " params " ) ;
2009-05-27 16:28:23 +02:00
double * dparams = ( double * ) mxGetData ( mxFldp ) ;
int npar = ( int ) mxGetM ( mxFldp ) ;
2010-03-09 17:09:18 +01:00
Vector modParams ( dparams , npar ) ;
2009-05-27 16:28:23 +02:00
2010-08-30 17:11:58 +02:00
mxFldp = mxGetField ( M_ , 0 , " Sigma_e " ) ;
2009-05-06 12:10:27 +02:00
dparams = ( double * ) mxGetData ( mxFldp ) ;
2009-05-27 16:28:23 +02:00
npar = ( int ) mxGetN ( mxFldp ) ;
2010-03-09 17:09:18 +01:00
TwoDMatrix vCov ( npar , npar , dparams ) ;
2009-05-27 16:28:23 +02:00
2010-08-30 17:11:58 +02:00
mxFldp = mxGetField ( dr , 0 , " ys " ) ; // and not in order of dr.order_var
2009-05-06 12:10:27 +02:00
dparams = ( double * ) mxGetData ( mxFldp ) ;
2009-05-27 16:28:23 +02:00
const int nSteady = ( int ) mxGetM ( mxFldp ) ;
2010-03-09 17:09:18 +01:00
Vector ySteady ( dparams , nSteady ) ;
2009-05-27 16:28:23 +02:00
mxFldp = mxGetField ( dr , 0 , " nstatic " ) ;
const int nStat = ( int ) mxGetScalar ( mxFldp ) ;
mxFldp = mxGetField ( dr , 0 , " npred " ) ;
int nPred = ( int ) mxGetScalar ( mxFldp ) ;
mxFldp = mxGetField ( dr , 0 , " nspred " ) ;
const int nsPred = ( 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 ( dr , 0 , " nsfwrd " ) ;
const int nsForw = ( 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 ) ;
2009-08-24 18:01:25 +02:00
2009-05-06 12:10:27 +02:00
nPred - = nBoth ; // correct nPred for nBoth.
2009-05-27 16:28:23 +02:00
2010-08-30 17:11:58 +02:00
mxFldp = mxGetField ( dr , 0 , " order_var " ) ;
2009-05-06 12:10:27 +02:00
dparams = ( double * ) mxGetData ( mxFldp ) ;
2009-05-27 16:28:23 +02:00
npar = ( int ) mxGetM ( mxFldp ) ;
2010-08-30 17:11:58 +02:00
if ( npar ! = nEndo )
2009-12-16 18:18:38 +01:00
mexErrMsgTxt ( " Incorrect number of input var_order vars. " ) ;
2010-03-09 17:09:18 +01:00
vector < int > var_order_vp ( nEndo ) ;
2009-05-27 16:28:23 +02:00
for ( int v = 0 ; v < nEndo ; v + + )
2010-03-09 17:09:18 +01:00
var_order_vp [ v ] = ( int ) ( * ( dparams + + ) ) ;
2009-05-27 16:28:23 +02:00
2009-05-06 12:10:27 +02:00
// the lag, current and lead blocks of the jacobian respectively
2010-08-30 17:11:58 +02:00
mxFldp = mxGetField ( M_ , 0 , " lead_lag_incidence " ) ;
2009-05-06 12:10:27 +02:00
dparams = ( double * ) mxGetData ( mxFldp ) ;
2009-05-27 16:28:23 +02:00
npar = ( int ) mxGetN ( mxFldp ) ;
int nrows = ( int ) mxGetM ( mxFldp ) ;
2009-11-29 21:50:39 +01:00
2010-03-09 17:09:18 +01:00
TwoDMatrix llincidence ( nrows , npar , dparams ) ;
2009-05-27 16:28:23 +02:00
if ( npar ! = nEndo )
2009-11-30 17:31:27 +01:00
mexErrMsgIdAndTxt ( " dynare:k_order_perturbation " , " Incorrect length of lead lag incidences: ncol=%d != nEndo=%d. " , npar , nEndo ) ;
2009-08-24 18:01:25 +02:00
2009-12-16 18:18:38 +01:00
//get NNZH =NNZD(2) = the total number of non-zero Hessian elements
2009-08-24 18:01:25 +02:00
mxFldp = mxGetField ( M_ , 0 , " NNZDerivatives " ) ;
dparams = ( double * ) mxGetData ( mxFldp ) ;
2010-03-09 17:09:18 +01:00
Vector NNZD ( dparams , ( int ) mxGetM ( mxFldp ) ) ;
2010-07-17 10:14:22 +02:00
if ( NNZD [ kOrder - 1 ] = = - 1 )
mexErrMsgTxt ( " The derivatives were not computed for the required order. Make sure that you used the right order option inside the stoch_simul command " ) ;
2009-05-06 12:10:27 +02:00
const int jcols = nExog + nEndo + nsPred + nsForw ; // Num of Jacobian columns
2009-05-27 16:28:23 +02:00
mxFldp = mxGetField ( M_ , 0 , " var_order_endo_names " ) ;
const int nendo = ( int ) mxGetM ( mxFldp ) ;
const int widthEndo = ( int ) mxGetN ( mxFldp ) ;
2010-03-09 18:22:33 +01:00
vector < string > endoNames ;
DynareMxArrayToString ( mxFldp , nendo , widthEndo , endoNames ) ;
2009-05-27 16:28:23 +02:00
2010-08-30 17:11:58 +02:00
mxFldp = mxGetField ( M_ , 0 , " exo_names " ) ;
2009-05-27 16:28:23 +02:00
const int nexo = ( int ) mxGetM ( mxFldp ) ;
const int widthExog = ( int ) mxGetN ( mxFldp ) ;
2010-03-09 18:22:33 +01:00
vector < string > exoNames ;
DynareMxArrayToString ( mxFldp , nexo , widthExog , exoNames ) ;
2009-05-27 16:28:23 +02:00
2010-03-09 15:20:18 +01:00
if ( ( nEndo ! = nendo ) | | ( nExog ! = nexo ) )
2009-11-30 17:31:27 +01:00
mexErrMsgTxt ( " Incorrect number of input parameters. " ) ;
2009-05-27 16:28:23 +02:00
2009-05-06 12:10:27 +02:00
/* Fetch time index */
2009-05-27 16:28:23 +02:00
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 ( ) ) ;
2010-09-20 19:24:23 +02:00
DynamicModelDLL dynamicDLL ( fName , nExog , dfExt ) ;
2009-05-27 16:28:23 +02:00
// intiate tensor library
tls . init ( kOrder , nStat + 2 * nPred + 3 * nBoth + 2 * nForw + nExog ) ;
// make KordpDynare object
2010-03-09 18:22:33 +01:00
KordpDynare dynare ( endoNames , nEndo , exoNames , nExog , nPar ,
2009-05-27 16:28:23 +02:00
ySteady , vCov , modParams , nStat , nPred , nForw , nBoth ,
2009-12-16 18:18:38 +01:00
jcols , NNZD , nSteps , kOrder , journal , dynamicDLL ,
2009-06-30 09:33:40 +02:00
sstol , var_order_vp , llincidence , qz_criterium ) ;
2009-05-27 16:28:23 +02:00
// construct main K-order approximation class
2009-12-16 18:18:38 +01:00
2009-05-27 16:28:23 +02:00
Approximation app ( dynare , journal , nSteps , false , qz_criterium ) ;
// run stochastic steady
app . walkStochSteady ( ) ;
/* Write derivative outputs into memory map */
map < string , ConstTwoDMatrix > mm ;
2009-09-23 15:14:05 +02:00
app . getFoldDecisionRule ( ) . writeMMap ( mm , string ( ) ) ;
2009-05-27 16:28:23 +02:00
// get latest ysteady
2010-03-09 17:09:18 +01:00
ySteady = dynare . getSteady ( ) ;
2009-05-27 16:28:23 +02:00
2009-11-29 21:50:39 +01:00
if ( kOrder = = 1 )
2009-05-27 16:28:23 +02:00
{
/* Set the output pointer to the output matrix ysteady. */
2009-12-16 18:18:38 +01:00
map < string , ConstTwoDMatrix > : : const_iterator cit = mm . begin ( ) ;
+ + cit ;
2009-11-29 21:50:39 +01:00
plhs [ 0 ] = mxCreateDoubleMatrix ( ( * cit ) . second . numRows ( ) , ( * cit ) . second . numCols ( ) , mxREAL ) ;
2010-02-08 16:58:24 +01:00
// Copy Dynare++ matrix into MATLAB matrix
const ConstVector & vec = ( * cit ) . second . getData ( ) ;
assert ( vec . skip ( ) = = 1 ) ;
memcpy ( mxGetPr ( plhs [ 0 ] ) , vec . base ( ) , vec . length ( ) * sizeof ( double ) ) ;
2009-05-27 16:28:23 +02:00
}
2009-11-29 21:50:39 +01:00
if ( kOrder > = 2 )
2009-05-27 16:28:23 +02:00
{
2009-11-29 21:50:39 +01:00
int ii = 0 ;
2009-05-27 16:28:23 +02:00
for ( map < string , ConstTwoDMatrix > : : const_iterator cit = mm . begin ( ) ;
( ( cit ! = mm . end ( ) ) & & ( ii < nlhs ) ) ; + + cit )
{
{
plhs [ ii ] = mxCreateDoubleMatrix ( ( * cit ) . second . numRows ( ) , ( * cit ) . second . numCols ( ) , mxREAL ) ;
2010-02-08 16:58:24 +01:00
// Copy Dynare++ matrix into MATLAB matrix
const ConstVector & vec = ( * cit ) . second . getData ( ) ;
assert ( vec . skip ( ) = = 1 ) ;
memcpy ( mxGetPr ( plhs [ ii ] ) , vec . base ( ) , vec . length ( ) * sizeof ( double ) ) ;
2009-05-27 16:28:23 +02:00
+ + ii ;
}
}
}
2009-05-06 12:10:27 +02:00
}
2009-05-27 16:28:23 +02:00
catch ( const KordException & e )
2009-05-06 12:10:27 +02:00
{
2009-05-27 16:28:23 +02:00
e . print ( ) ;
2009-11-30 17:31:27 +01:00
mexErrMsgIdAndTxt ( " dynare:k_order_perturbation " , " Caught Kord exception: %s " , e . get_message ( ) ) ;
2009-05-06 12:10:27 +02:00
}
2009-05-27 16:28:23 +02:00
catch ( const TLException & e )
{
e . print ( ) ;
2009-11-30 17:31:27 +01:00
mexErrMsgIdAndTxt ( " dynare:k_order_perturbation " , " Caught TL exception " ) ;
2009-05-06 12:10:27 +02:00
}
2009-05-27 16:28:23 +02:00
catch ( SylvException & e )
{
e . printMessage ( ) ;
2009-11-30 17:31:27 +01:00
mexErrMsgIdAndTxt ( " dynare:k_order_perturbation " , " Caught Sylv exception " ) ;
2009-05-06 12:10:27 +02:00
}
2009-05-27 16:28:23 +02:00
catch ( const DynareException & e )
{
2009-11-30 17:31:27 +01:00
mexErrMsgIdAndTxt ( " dynare:k_order_perturbation " , " Caught KordDynare exception: %s " , e . message ( ) ) ;
2009-05-06 12:10:27 +02:00
}
2009-05-27 16:28:23 +02:00
catch ( const ogu : : Exception & e )
{
2009-11-30 17:31:27 +01:00
mexErrMsgIdAndTxt ( " dynare:k_order_perturbation " , " Caught general exception: %s " , e . message ( ) ) ;
}
} // end of mexFunction()
} // end of extern C
2009-11-29 21:50:39 +01:00
2009-08-24 18:01:25 +02:00
# endif // ifdef MATLAB_MEX_FILE to exclude mexFunction for other applications