minor bug fix and removed a dependency

- resolved bug in boot.mex which caused the sampling to not be deterministic for a given seed (when u = false only)
- removed C+11 dependency for compiling boot.cpp to make boot.mex

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
 

src/boot.cpp

@@ -10,32 +10,29 @@
 // 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
 // divisible by the sample size (n), then the sample index omitted is
 // 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 <vector>
-#include <random>
 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<long long int> 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<long long int> 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<int> 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<int> 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 ) {