diff --git a/inst/boot.m b/inst/boot.m index d91b4c9d..7e1c283c 100644 --- a/inst/boot.m +++ b/inst/boot.m @@ -149,6 +149,6 @@ %! %! % N as input; balanced resampling with replacement; setting the random seed %! boot(3,20,false,[],1) % Set random seed -%! boot(3,20,false,[],1) % Reset random seed, BOOTSAM is the same (if running on the same core) +%! boot(3,20,false,[],1) % Reset random seed, BOOTSAM is the same %! boot(3,20,false,[]) % Without setting random seed, BOOTSAM is different diff --git a/precompiled_binaries/matlab/linux/glnxa64/boot.mexa64 b/precompiled_binaries/matlab/linux/glnxa64/boot.mexa64 index ff05a506..d08aa9af 100755 Binary files a/precompiled_binaries/matlab/linux/glnxa64/boot.mexa64 and b/precompiled_binaries/matlab/linux/glnxa64/boot.mexa64 differ diff --git a/precompiled_binaries/matlab/mac/maci64/boot.mexmaci64 b/precompiled_binaries/matlab/mac/maci64/boot.mexmaci64 index 8d733ac2..ac92f310 100755 Binary files a/precompiled_binaries/matlab/mac/maci64/boot.mexmaci64 and b/precompiled_binaries/matlab/mac/maci64/boot.mexmaci64 differ diff --git a/precompiled_binaries/matlab/win/win32/boot.mex b/precompiled_binaries/matlab/win/win32/boot.mex new file mode 100644 index 00000000..06376899 Binary files /dev/null and b/precompiled_binaries/matlab/win/win32/boot.mex differ diff --git a/precompiled_binaries/matlab/win/win32/boot.mexw32 b/precompiled_binaries/matlab/win/win32/boot.mexw32 index e29561d9..a7cf7064 100644 Binary files a/precompiled_binaries/matlab/win/win32/boot.mexw32 and b/precompiled_binaries/matlab/win/win32/boot.mexw32 differ diff --git a/precompiled_binaries/matlab/win/win32/smoothmedian.mexw32 b/precompiled_binaries/matlab/win/win32/smoothmedian.mexw32 index 648b96f7..4cf68a83 100644 Binary files a/precompiled_binaries/matlab/win/win32/smoothmedian.mexw32 and b/precompiled_binaries/matlab/win/win32/smoothmedian.mexw32 differ diff --git a/precompiled_binaries/matlab/win/win64/boot.mexw64 b/precompiled_binaries/matlab/win/win64/boot.mexw64 index 6cceff1b..64cf7095 100644 Binary files a/precompiled_binaries/matlab/win/win64/boot.mexw64 and b/precompiled_binaries/matlab/win/win64/boot.mexw64 differ diff --git a/precompiled_binaries/matlab/win/win64/smoothmedian.mexw64 b/precompiled_binaries/matlab/win/win64/smoothmedian.mexw64 index 1813f39b..7688f137 100644 Binary files a/precompiled_binaries/matlab/win/win64/smoothmedian.mexw64 and b/precompiled_binaries/matlab/win/win64/smoothmedian.mexw64 differ diff --git a/precompiled_binaries/octave/linux/glnxa64/boot.mex b/precompiled_binaries/octave/linux/glnxa64/boot.mex index 2f4f876f..7c862621 100755 Binary files a/precompiled_binaries/octave/linux/glnxa64/boot.mex and b/precompiled_binaries/octave/linux/glnxa64/boot.mex differ diff --git a/precompiled_binaries/octave/mac/maci64/boot.mex b/precompiled_binaries/octave/mac/maci64/boot.mex index b44e3af3..d4fcf675 100755 Binary files a/precompiled_binaries/octave/mac/maci64/boot.mex and b/precompiled_binaries/octave/mac/maci64/boot.mex differ diff --git a/precompiled_binaries/octave/mac/maci64/smoothmedian.mex b/precompiled_binaries/octave/mac/maci64/smoothmedian.mex index d8b09adc..edb1851c 100755 Binary files a/precompiled_binaries/octave/mac/maci64/smoothmedian.mex and b/precompiled_binaries/octave/mac/maci64/smoothmedian.mex differ diff --git a/precompiled_binaries/octave/win/win64/boot.mex b/precompiled_binaries/octave/win/win64/boot.mex index 12bd2cce..54e3d39b 100644 Binary files a/precompiled_binaries/octave/win/win64/boot.mex and b/precompiled_binaries/octave/win/win64/boot.mex differ diff --git a/precompiled_binaries/octave/win/win64/smoothmedian.mex b/precompiled_binaries/octave/win/win64/smoothmedian.mex index 96e64f96..732fff14 100644 Binary files a/precompiled_binaries/octave/win/win64/smoothmedian.mex and b/precompiled_binaries/octave/win/win64/smoothmedian.mex differ diff --git a/src/boot.cpp b/src/boot.cpp index 23355e7a..9657928a 100644 --- a/src/boot.cpp +++ b/src/boot.cpp @@ -10,21 +10,20 @@ // bootsam = boot (x, nboot) // bootsam = boot (..., nboot, u) // bootsam = boot (..., nboot, u, w) -// bootsam = boot (..., nboot, u, w, s) +// bootsam = boot (..., nboot, u, w, seed) // // INPUT VARIABLES // n (double) is the number of rows (of the data vector) // x (double) is a data vector intended for resampling // nboot (double) is the number of bootstrap resamples // u (boolean) for unbiased: false (for bootstrap) or true (for bootknife) -// w (double) is a weight vector of length n. -// s (double) is a seed for the pseudo-random number generator. +// w (double) is a weight vector of length n +// seed (double) is a seed used to initialise the pseudo-random number generator // // OUTPUT VARIABLE // bootsam (double) is an n x nboot matrix of resampled data or indices // // NOTES -// Uniform random numbers are generated by the Mersenne Twister 19937 generator. // u is an optional input argument. The default is false. If u is true then // the sample index for omission in each bootknife resample is selected // systematically. If the remaining number of bootknife resamples is not @@ -32,10 +31,8 @@ // selected randomly. // w is an optional input argument. If w is empty or not provided, the default // is a vector of each element equal to nboot (i.e. uniform weighting). Each -// element of w is the number of times that the corresponding index is represented -// in bootsam. For example, if the second element is 500, then the value 2 will -// be assigned to 500 elements within bootsam. Therefore, the sum of w should -// equal n * nboot. +// element of w is the number of times that the corresponding index is +// represented in bootsam. Therefore, the sum of w should equal n * nboot. // Note that the mex function compiled from this source code is not thread // safe. Below is an example of a line of code one can run in Octave/Matlab // before attempting parallel operation of boot.mex in order to ensure that @@ -45,14 +42,12 @@ // In Matlab: // >> ncpus = feature('numcores'); parfor i = 1:ncpus; boot (1, 1, false, [], i); end; // -// Requirements: Compilation requires C++11 // // Author: Andrew Charles Penn (2022) #include "mex.h" #include -#include using namespace std; @@ -112,20 +107,20 @@ void mexFunction (int nlhs, mxArray* plhs[], if (mxGetNumberOfElements (prhs[2]) > 1 || !mxIsClass (prhs[2], "logical")) { mexErrMsgTxt ("the third input argument (u) must be a logical scalar value"); } - u = *(mxGetPr (prhs[2])); + u = *(mxGetLogicals (prhs[2])); } // Fourth input argument (w) - error checking is handled later (see below) // Fifth input argument (s) - unsigned int s; + unsigned int seed; if ( nrhs > 4 ) { if ( mxGetNumberOfElements (prhs[4]) > 1 ) { mexErrMsgTxt ("the fifth input argument (s) must be a scalar value"); } - s = *(mxGetPr(prhs[4])); - if ( !mxIsFinite (s) ) { + seed = *(mxGetPr(prhs[4])); + if ( !mxIsFinite (seed) ) { mexErrMsgTxt ("the fifth input argument (s) cannot be NaN or Inf"); } - srand (s); + srand (seed); } // Output variables @@ -137,11 +132,11 @@ void mexFunction (int nlhs, mxArray* plhs[], mwSize dims[2] = {n, nboot}; plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, - mxREAL); // Prepare array for bootstrap sample indices - long long int N = n * nboot; // Total counts of all sample indices - long long int k; // Variable to store random number - long long int d; // Counter for cumulative sum calculations - vector c; // Counter for each of the sample indices + mxREAL); // Prepare array for bootstrap sample indices + long long int N = n * nboot; // Total counts of all sample indices + long long int k; // Variable to store random number + long long int d; // Counter for cumulative sum calculations + vector c; // Counter for each of the sample indices c.reserve (n); if ( nrhs > 3 && !mxIsEmpty (prhs[3]) ) { // Assign user defined weights (counts) @@ -163,7 +158,7 @@ void mexFunction (int nlhs, mxArray* plhs[], if ( w[i] < 0 ) { mexErrMsgTxt ("the fourth input argument (weights) must contain only positive integers"); } - c.push_back (w[i]); // Set each element in c to the specified weight + c.push_back (w[i]); // Set each element in c to the specified weight s += c[i]; } if ( s != N ) { @@ -175,36 +170,33 @@ void mexFunction (int nlhs, mxArray* plhs[], c.push_back (nboot); // Set each element in c to nboot } } - bool LOO = false; // Leave-one-out (LOO) flag for the current bootstrap iteration (remains false if u is false) - long long int m = 0; // Counter for LOO sample index r (remains 0 if u is false) - int r = -1; // Sample index for LOO (remains -1 and is ignored if u is false) + bool LOO = false; // Leave-one-out (LOO) flag + long long int m = 0; // Counter for LOO sample index r + int r = -1; // Sample index for LOO // Create pointer so that we can access elements of bootsam (i.e. plhs[0]) double *ptr = (double *) mxGetData(plhs[0]); - - // Initialize pseudo-random number generator (Mersenne Twister 19937) - mt19937 rng (rand()); - uniform_int_distribution distr (0, n - 1); // Perform balanced sampling for ( int b = 0; b < nboot ; b++ ) { if (u) { if ( (b / n) == (nboot / n) ) { - r = distr (rng); // random + r = rand () / (RAND_MAX / n + 1); // random } else { - r = b - (b / n) * n; // systematic + r = b - (b / n) * n; // systematic } } for ( int i = 0; i < n ; i++ ) { if (u) { - if (c[r] < N) { // Only LOO if sample index r doesn't account for all remaining sampling counts + // Only LOO if sample index r doesn't account for all remaining + // sampling counts + if (c[r] < N) { m = c[r]; c[r] = 0; LOO = true; } } - uniform_int_distribution distk (0, N - m - 1); - k = distk (rng); + k = rand () / (RAND_MAX / (N - m) + 1); d = c[0]; for ( int j = 0; j < n ; j++ ) { if ( k < d ) {