Skip to content
Snippets Groups Projects
Verified Commit fa64a258 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

qmc_sequence DLL: various modernizations and simplifications

parent d0b1e0a7
No related branches found
No related tags found
No related merge requests found
Pipeline #1090 passed
......@@ -22,18 +22,16 @@
** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
*/
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <limits>
#include <algorithm>
#include <dynblas.h>
using namespace std;
#define lb .02425
#define ub .97575
constexpr double lb = .02425;
constexpr double ub = .97575;
#ifdef USE_OMP
# include <omp.h>
......@@ -81,8 +79,8 @@ icdf(const T uniform)
2.445134137142996e+00,
3.754408661907416e+00
};
T gaussian = (T) 0.0;
if ((0 < uniform) && (uniform < lb))
T gaussian = static_cast<T>(0.0);
if (0 < uniform && uniform < lb)
{
T tmp;
tmp = sqrt(-2*log(uniform));
......@@ -90,7 +88,7 @@ icdf(const T uniform)
}
else
{
if ((lb <= uniform) && (uniform <= ub))
if (lb <= uniform && uniform <= ub)
{
T tmp, TMP;
tmp = uniform - .5;
......@@ -99,7 +97,7 @@ icdf(const T uniform)
}
else
{
if ((ub < uniform) && (uniform < 1))
if (ub < uniform && uniform < 1)
{
T tmp;
tmp = sqrt(-2*log(1-uniform));
......@@ -107,7 +105,7 @@ icdf(const T uniform)
}
}
}
if ((0 < uniform) && (uniform < 1))
if (0 < uniform && uniform < 1)
{
T tmp, tmp_;
tmp = .5*erfc(-gaussian/sqrt(2.0))-uniform;
......@@ -115,33 +113,29 @@ icdf(const T uniform)
gaussian = gaussian - tmp_/(1+.5*gaussian*tmp_);
}
if (uniform == 0)
{
gaussian = -INFINITY;
}
gaussian = -numeric_limits<T>::infinity();
if (uniform == 1)
{
gaussian = INFINITY;
}
return (gaussian);
gaussian = numeric_limits<T>::infinity();
return gaussian;
}
template<typename T>
void
icdfm(const int n, T *U)
icdfm(int n, T *U)
{
#if USE_OMP
# pragma omp parallel for num_threads(omp_get_num_threads())
#endif
for (int i = 0; i < n; i++)
{
U[i] = icdf(U[i]);
}
U[i] = icdf(U[i]);
return;
}
template<typename T>
void
icdfmSigma(const int d, const int n, T *U, const double *LowerCholSigma)
icdfmSigma(int d, int n, T *U, const double *LowerCholSigma)
{
double one = 1.0;
double zero = 0.0;
......@@ -150,13 +144,12 @@ icdfmSigma(const int d, const int n, T *U, const double *LowerCholSigma)
icdfm(n*d, U);
double tmp[n*d];
dgemm("N", "N", &dd, &nn, &dd, &one, LowerCholSigma, &dd, U, &dd, &zero, tmp, &dd);
memcpy(U, tmp, d*n*sizeof(double));
return;
copy_n(tmp, d*n, U);
}
template<typename T>
void
usphere(const int d, const int n, T *U)
usphere(int d, int n, T *U)
{
icdfm(n*d, U);
#if USE_OMP
......@@ -167,21 +160,17 @@ usphere(const int d, const int n, T *U)
int k = j*d;
double norm = 0.0;
for (int i = 0; i < d; i++)// dimension index.
{
norm = norm + U[k+i]*U[k+i];
}
norm = norm + U[k+i]*U[k+i];
norm = sqrt(norm);
for (int i = 0; i < d; i++)// dimension index.
{
U[k+i] = U[k+i]/norm;
}
U[k+i] = U[k+i]/norm;
}
return;
}
template<typename T>
void
usphereRadius(const int d, const int n, double radius, T *U)
usphereRadius(int d, int n, double radius, T *U)
{
icdfm(n*d, U);
#if USE_OMP
......@@ -192,14 +181,10 @@ usphereRadius(const int d, const int n, double radius, T *U)
int k = j*d;
double norm = 0.0;
for (int i = 0; i < d; i++)// dimension index.
{
norm = norm + U[k+i]*U[k+i];
}
norm = norm + U[k+i]*U[k+i];
norm = sqrt(norm);
for (int i = 0; i < d; i++)// dimension index.
{
U[k+i] = radius*U[k+i]/norm;
}
U[k+i] = radius*U[k+i]/norm;
}
return;
}
......@@ -19,17 +19,13 @@
** along with Dynare. If not, see <http://www.gnu.org/licenses/>.
**/
#include <sstream>
#include <string.h>
#include <stdint.h>
#include <string>
#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[])
{
......@@ -38,16 +34,16 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
** 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.
** 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),
** 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[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
**
......@@ -55,136 +51,104 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
/*
** 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 (nrhs < 3 || nrhs > 5)
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!");
}
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]);
if (!mxIsNumeric(prhs[0]))
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be a positive integer!");
int dimension = static_cast<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]);
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Second input (seed) has to be an integer [int64]!");
int64_T seed = static_cast<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 (!mxIsNumeric(prhs[2]))
error_flag_3 = 1;
int type = static_cast<int>(mxGetScalar(prhs[2]));
if (!(type == 0 || type == 1 || type == 2))
{
error_flag_3 = 1;
}
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!");
}
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
** Test dimension2 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!");
}
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());
}
DYN_MEX_FUNC_ERR_MSG_TXT(("qmc_sequence:: First input (dimension) has to be smaller than " + to_string(DIM_MAX) + " !").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!");
}
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]);
}
sequence_size = static_cast<int>(mxGetScalar(prhs[3]));
else
{
sequence_size = 1;
}
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;
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 && static_cast<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 && !(static_cast<int>(mxGetN(prhs[4])) == dimension
&& static_cast<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 a hypershere
DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a positive scalar!");
const double *lower_bounds = nullptr, *upper_bounds = nullptr;
int unit_hypercube_flag = 1;
if ((type == 0) && (nrhs > 4))
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];
lower_bounds = mxGetPr(prhs[4]);
upper_bounds = lower_bounds + dimension;
unit_hypercube_flag = 0;
}
double *cholcov;
const double *cholcov = nullptr;
int identity_covariance_matrix = 1;
if ((type == 1) && (nrhs > 4))
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];
cholcov = mxGetPr(prhs[4]);
identity_covariance_matrix = 0;
}
double radius = 1.0;
int unit_radius = 1;
if ((type == 2) && (nrhs > 4))
if (type == 2 && nrhs > 4)
{
double *tmp;
tmp = (double *) mxCalloc(1, sizeof(double));
memcpy(tmp, mxGetPr(prhs[4]), sizeof(double));
radius = tmp[0];
radius = *mxGetPr(prhs[4]);
unit_radius = 0;
}
/*
** Initialize outputs of the mex file.
*/
double *qmc_draws;
plhs[0] = mxCreateDoubleMatrix(dimension, sequence_size, mxREAL);
qmc_draws = mxGetPr(plhs[0]);
double *qmc_draws = mxGetPr(plhs[0]);
int64_T seed_out;
if (sequence_size == 1)
......@@ -195,16 +159,16 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
else
seed_out = sobol_block(dimension, sequence_size, seed, qmc_draws);
if (type == 0 && unit_hypercube_flag == 0) // Uniform QMC sequence in an hypercube.
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.
else if (type == 1) // Normal QMC sequence in ℝⁿ
{
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.
else if (type == 2) // Uniform QMC sequence on a hypershere
{
if (unit_radius == 1)
usphere(dimension, sequence_size, qmc_draws);
......@@ -215,7 +179,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (nlhs >= 2)
{
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT64_CLASS, mxREAL);
*((int64_T *) mxGetData(plhs[1])) = seed_out;
*(reinterpret_cast<int64_T *>(mxGetData(plhs[1]))) = seed_out;
}
if (nlhs >= 3)
plhs[2] = mxCreateDoubleScalar(0);
......
......@@ -29,7 +29,7 @@
using namespace std;
#define DIM_MAX 1111
constexpr int DIM_MAX = 1111;
template<typename T>
int
......@@ -130,9 +130,8 @@ bit_lo0(T n)
bit++;
T n2 = n/2;
if (n == 2*n2)
{
break;
}
break;
n = n2;
}
return bit;
......@@ -157,8 +156,7 @@ ixor(T i, T j)
{
T i2 = i / 2;
T j2 = j / 2;
if (
((i == 2 * i2) && (j != 2 * j2))
if (((i == 2 * i2) && (j != 2 * j2))
|| ((i != 2 * i2) && (j == 2 * j2)))
{
k = k + l;
......@@ -392,7 +390,7 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
/*
** Set ATMOST = 2^LOG_MAX - 1.
*/
atmost = (T1) 0;
atmost = static_cast<T1>(0);
for (int i = 1; i <= LOG_MAX; i++)
atmost = 2 * atmost + 1;
/*
......@@ -404,7 +402,7 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
*/
for (T1 j = 0; j < maxcol; j++)
{
v[0][j] = (T1) 1;
v[0][j] = static_cast<T1>(1);
}
/*
** Initialize the remaining rows of V.
......@@ -423,9 +421,8 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
{
j = j / 2;
if (j <= 0)
{
break;
}
break;
m = m + 1;
}
/*
......@@ -453,9 +450,7 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
{
l = 2 * l;
if (includ[k])
{
newv = (newv ^ (l * v[i][j-k-1]));
}
newv = newv ^ (l * v[i][j-k-1]);
}
v[i][j] = newv;
}
......@@ -468,14 +463,12 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
{
l = 2 * l;
for (int i = 0; i < dim_num; i++)
{
v[i][j] = v[i][j] * l;
}
v[i][j] = v[i][j] * l;
}
/*
** RECIPD is 1/(common denominator of the elements in V).
*/
recipd = 1.0E+00 / ((T2) (2 * l));
recipd = 1.0E+00 / static_cast<T2>(2 * l);
}
if (*seed < 0)
*seed = 0;
......@@ -484,14 +477,10 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
{
l = 1;
for (int i = 0; i < dim_num; i++)
{
lastq[i] = 0;
}
lastq[i] = 0;
}
else if (*seed == seed_save + 1)
{
l = bit_lo0(*seed);
}
l = bit_lo0(*seed);
else if (*seed <= seed_save)
{
seed_save = 0;
......@@ -502,9 +491,7 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
{
l = bit_lo0(seed_temp);
for (int i = 0; i < dim_num; i++)
{
lastq[i] = (lastq[i] ^ v[i][l-1]);
}
lastq[i] = (lastq[i] ^ v[i][l-1]);
}
l = bit_lo0(*seed);
}
......@@ -514,9 +501,7 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
{
l = bit_lo0(seed_temp);
for (int i = 0; i < dim_num; i++)
{
lastq[i] = (lastq[i] ^ v[i][l-1]);
}
lastq[i] = (lastq[i] ^ v[i][l-1]);
}
l = bit_lo0(*seed);
}
......@@ -539,8 +524,8 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi)
*/
for (int i = 0; i < dim_num; i++)
{
quasi[i] = ((T2) lastq[i]) * recipd;
lastq[i] = (lastq[i]^v[i][l-1]);
quasi[i] = static_cast<T2>(lastq[i]) * recipd;
lastq[i] = lastq[i]^v[i][l-1];
}
seed_save = *seed;
*seed = *seed + 1;
......@@ -552,31 +537,24 @@ T1
sobol_block(int dimension, int block_size, T1 seed, T2 *block)
{
for (int iter = 0; iter < block_size; iter++)
{
next_sobol(dimension, &seed, &block[iter*dimension]);
}
next_sobol(dimension, &seed, &block[iter*dimension]);
return seed;
}
template<typename T>
void
expand_unit_hypercube(int dimension, int block_size, T *block, T *lower_bound, T *upper_bound)
expand_unit_hypercube(int dimension, int block_size, T *block, const T *lower_bound, const T *upper_bound)
{
T *hypercube_length = new T[dimension];
for (int dim = 0; dim < dimension; dim++)
{
hypercube_length[dim] = upper_bound[dim]-lower_bound[dim];
}
hypercube_length[dim] = upper_bound[dim]-lower_bound[dim];
int base = 0;
for (int sim = 0; sim < block_size; sim++)
{
for (int dim = 0; dim < dimension; dim++)
{
block[base+dim] = lower_bound[dim] + hypercube_length[dim]*block[base+dim];
}
block[base+dim] = lower_bound[dim] + hypercube_length[dim]*block[base+dim];
base += dimension;
}
delete[] hypercube_length;
}
#undef DIM_MAX
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment