224 lines
7.8 KiB
C++
224 lines
7.8 KiB
C++
/*
|
|
** Computes Quasi Monte-Carlo sequence.
|
|
**
|
|
** Copyright (C) 2010-2012 Dynare Team
|
|
**
|
|
** This file is part of Dynare (can be used outside 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/>.
|
|
**
|
|
** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
|
|
**/
|
|
|
|
#include <sstream>
|
|
#include <string.h>
|
|
#include <stdint.h>
|
|
#include <dynmex.h>
|
|
|
|
#include "sobol.hh"
|
|
#include "gaussian.hh"
|
|
|
|
// the maximum dimension defined in sobol.ff (but undef at the end)
|
|
#define DIM_MAX 1111
|
|
|
|
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|
{
|
|
/*
|
|
** INPUTS:
|
|
** prhs[0] [integer] scalar, dimension.
|
|
** prhs[1] [integer] scalar, seed.
|
|
** prhs[2] [integer] scalar, sequence type:
|
|
** 0 ==> uniform,
|
|
** 1 ==> gaussian,
|
|
** 2 ==> uniform on an hypershere.
|
|
** prhs[3] [integer] scalar, sequence size.
|
|
** prhs[4] [double] dimension*2 array, lower and upper bounds of the hypercube (default is 0-1 in all dimensions) if prhs[2]==0,
|
|
** dimension*dimension array, covariance of the multivariate gaussian distribution of prhs[2]==1 (default is the identity matrix),
|
|
** scalar, radius of the hypershere if prhs[2]==2 (default is one).
|
|
**
|
|
** OUTPUTS:
|
|
** plhs[0] [double] sequence_size*dimension array, the Sobol sequence.
|
|
** plhs[1] [integer] scalar, seed.
|
|
** plhs[2] [integer] zero in case of success, one in case of error
|
|
**
|
|
*/
|
|
/*
|
|
** Check the number of input and output arguments.
|
|
*/
|
|
if ( !( (nrhs==5) | (nrhs==4) | (nrhs==3) ) )
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Five, four or three input arguments are required!");
|
|
}
|
|
if (nlhs == 0)
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: At least one output argument is required!");
|
|
}
|
|
/*
|
|
** Test the first input argument and assign it to dimension.
|
|
*/
|
|
if ( !( mxIsNumeric(prhs[0]) ) )
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be a positive integer!");
|
|
}
|
|
int dimension = ( int ) mxGetScalar(prhs[0]);
|
|
/*
|
|
** Test the second input argument and assign it to seed.
|
|
*/
|
|
if ( !( mxIsNumeric(prhs[1]) && mxIsClass(prhs[1],"int64") ) )
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Second input (seed) has to be an integer [int64]!");
|
|
}
|
|
int64_T seed = ( int64_T ) mxGetScalar( prhs[1] );
|
|
/*
|
|
** Test the third input argument and assign it to type (kind of QMC sequence).
|
|
*/
|
|
int error_flag_3 = 0;
|
|
if ( !(mxIsNumeric(prhs[2])) )
|
|
{
|
|
error_flag_3 = 1;
|
|
}
|
|
int type = ( int ) mxGetScalar(prhs[2]);
|
|
if ( !(type==0 || type==1 || type==2) )
|
|
{
|
|
error_flag_3 = 1;
|
|
}
|
|
if (error_flag_3==1)
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Third input (type of QMC sequence) has to be an integer equal to 0, 1 or 2!");
|
|
}
|
|
/*
|
|
** Test dimension>=2 when type==2
|
|
*/
|
|
if ( (type==2) && (dimension<2) )
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be greater than 1 for a uniform QMC on an hypershere!");
|
|
}
|
|
else if (dimension > DIM_MAX)
|
|
{
|
|
stringstream msg;
|
|
msg << "qmc_sequence:: First input (dimension) has to be smaller than " << DIM_MAX << " !";
|
|
DYN_MEX_FUNC_ERR_MSG_TXT(msg.str().c_str());
|
|
}
|
|
|
|
/*
|
|
** Test the optional fourth input argument and assign it to sequence_size.
|
|
*/
|
|
if ( ( nrhs>3 ) && !mxIsNumeric(prhs[3]) )
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Fourth input (qmc sequence size) has to be a positive integer!");
|
|
}
|
|
int sequence_size;
|
|
if ( nrhs>3)
|
|
{
|
|
sequence_size = ( int ) mxGetScalar( prhs[3] );
|
|
}
|
|
else
|
|
{
|
|
sequence_size = 1;
|
|
}
|
|
/*
|
|
** Test the optional fifth input argument and assign it to lower_and_upper_bounds.
|
|
*/
|
|
if ( ( nrhs>4 ) && (type==0) && ( !( mxGetN(prhs[4])==2) ) )// Sequence of uniformly distributed numbers in an hypercube.
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be an array with two columns!");
|
|
}
|
|
if ( (nrhs>4) && (type==0) && ( !( (int)mxGetM(prhs[4])==dimension) ) )
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fourth input argument must be an array with a number of lines equal to dimension (first input argument)!");
|
|
}
|
|
if ( ( nrhs>4 ) && (type==1) && ( !( ((int)mxGetN(prhs[4])==dimension) && ((int)mxGetM(prhs[4])==dimension) ) ) )// Sequence of normally distributed numbers.
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a squared matrix (whose dimension is given by the first input argument)!");
|
|
}
|
|
if ( ( nrhs>4 ) && (type==2) && ( !( (mxGetN(prhs[4])==1) && (mxGetM(prhs[4])==1) ) ) )// Sequence of uniformly distributed numbers on an hypershere.
|
|
{
|
|
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a positive scalar!");
|
|
}
|
|
double *lower_bounds = NULL, *upper_bounds = NULL;
|
|
int unit_hypercube_flag = 1;
|
|
if ( (type==0) && (nrhs>4) )
|
|
{
|
|
lower_bounds = (double *) mxCalloc(dimension,sizeof(double));
|
|
upper_bounds = (double *) mxCalloc(dimension,sizeof(double));
|
|
double *tmp;
|
|
tmp = (double *) mxCalloc(dimension*2,sizeof(double));
|
|
memcpy(tmp,mxGetPr(prhs[4]),dimension*2*sizeof(double));
|
|
lower_bounds = &tmp[0];
|
|
upper_bounds = &tmp[dimension];
|
|
unit_hypercube_flag = 0;
|
|
}
|
|
double *cholcov;
|
|
int identity_covariance_matrix = 1;
|
|
if ( (type==1) && (nrhs>4) )
|
|
{
|
|
cholcov = (double *) mxCalloc(dimension*dimension,sizeof(double));
|
|
double *tmp;
|
|
tmp = (double *) mxCalloc(dimension*dimension,sizeof(double));
|
|
memcpy(tmp,mxGetPr(prhs[4]),dimension*dimension*sizeof(double));
|
|
cholcov = &tmp[0];
|
|
identity_covariance_matrix = 0;
|
|
}
|
|
double radius = 1.0;
|
|
int unit_radius = 1;
|
|
if ( (type==2) && (nrhs>4) )
|
|
{
|
|
double *tmp;
|
|
tmp = (double *) mxCalloc(1,sizeof(double));
|
|
memcpy(tmp,mxGetPr(prhs[4]),sizeof(double));
|
|
radius = tmp[0];
|
|
unit_radius = 0;
|
|
}
|
|
/*
|
|
** Initialize outputs of the mex file.
|
|
*/
|
|
double *qmc_draws;
|
|
plhs[0] = mxCreateDoubleMatrix(dimension,sequence_size,mxREAL);
|
|
qmc_draws = mxGetPr(plhs[0]);
|
|
int64_T seed_out;
|
|
|
|
if (sequence_size==1)
|
|
{
|
|
next_sobol ( dimension, &seed, qmc_draws );
|
|
seed_out = seed;
|
|
}
|
|
else
|
|
seed_out = sobol_block( dimension, sequence_size, seed, qmc_draws);
|
|
|
|
if (type==0 && unit_hypercube_flag==0)// Uniform QMC sequence in an hypercube.
|
|
expand_unit_hypercube( dimension, sequence_size, qmc_draws, lower_bounds, upper_bounds);
|
|
else if (type==1)// Normal QMC sequance in R^n.
|
|
{
|
|
if (identity_covariance_matrix==1)
|
|
icdfm(dimension*sequence_size, qmc_draws);
|
|
else
|
|
icdfmSigma(dimension,sequence_size, qmc_draws, cholcov);
|
|
}
|
|
else if (type==2)// Uniform QMC sequence on an hypershere.
|
|
{
|
|
if (unit_radius==1)
|
|
usphere(dimension, sequence_size, qmc_draws);
|
|
else
|
|
usphereRadius(dimension, sequence_size, radius, qmc_draws);
|
|
}
|
|
|
|
if (nlhs >= 2)
|
|
{
|
|
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT64_CLASS, mxREAL);
|
|
*((int64_T *) mxGetData(plhs[1])) = seed_out;
|
|
}
|
|
if (nlhs >= 3)
|
|
plhs[2] = mxCreateDoubleScalar(0);
|
|
}
|