From fa64a25825aaa059f9b5f0d9fb57343c5c8989a2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 3 Apr 2019 18:49:32 +0200
Subject: [PATCH] qmc_sequence DLL: various modernizations and simplifications
---
mex/sources/sobol/gaussian.hh | 67 +++++-------
mex/sources/sobol/qmc_sequence.cc | 168 ++++++++++++------------------
mex/sources/sobol/sobol.hh | 66 ++++--------
3 files changed, 114 insertions(+), 187 deletions(-)
diff --git a/mex/sources/sobol/gaussian.hh b/mex/sources/sobol/gaussian.hh
index d1fcd1b88a..fb44f589e1 100644
--- a/mex/sources/sobol/gaussian.hh
+++ b/mex/sources/sobol/gaussian.hh
@@ -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;
}
diff --git a/mex/sources/sobol/qmc_sequence.cc b/mex/sources/sobol/qmc_sequence.cc
index e99d805945..5b4d1e6f4d 100644
--- a/mex/sources/sobol/qmc_sequence.cc
+++ b/mex/sources/sobol/qmc_sequence.cc
@@ -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 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!");
- }
+ 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);
diff --git a/mex/sources/sobol/sobol.hh b/mex/sources/sobol/sobol.hh
index 8926daca92..b39795314c 100644
--- a/mex/sources/sobol/sobol.hh
+++ b/mex/sources/sobol/sobol.hh
@@ -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
--
GitLab