support sampling with and without replacement

NEWS.md

@@ -12,6 +12,8 @@
 
 ### enhancements
 
+* Support sampling with and without replacement in the bootstrapped version of pattern causality (#58).
+
 * Enable method dispatch compatibility in `pc` and `ops` generics via `...` (#43).
 
 * Improve handling of large-scale inputs (#33).

R/RcppExports.R

@@ -9,8 +9,8 @@ RcppPC <- function(target, source, lib, pred, E, tau, style = 0L, num_neighbors
     .Call(`_pc_RcppPC`, target, source, lib, pred, E, tau, style, num_neighbors, zero_tolerance, dist_metric, relative, weighted, threads, h, nb, nrows)
 }
 
-RcppPCboot <- function(target, source, libsizes, lib, pred, E, tau, style = 0L, num_neighbors = 4L, zero_tolerance = 0L, dist_metric = "euclidean", boot = 99L, random_sample = TRUE, seed = 42L, relative = TRUE, weighted = TRUE, threads = 1L, parallel_level = 0L, verbose = FALSE, h = 0L, nb = NULL, nrows = NULL) {
-    .Call(`_pc_RcppPCboot`, target, source, libsizes, lib, pred, E, tau, style, num_neighbors, zero_tolerance, dist_metric, boot, random_sample, seed, relative, weighted, threads, parallel_level, verbose, h, nb, nrows)
+RcppPCboot <- function(target, source, libsizes, lib, pred, E, tau, style = 0L, num_neighbors = 4L, zero_tolerance = 0L, dist_metric = "euclidean", boot = 99L, replace_sampling = TRUE, seed = 42L, relative = TRUE, weighted = TRUE, threads = 1L, parallel_level = 0L, verbose = FALSE, h = 0L, nb = NULL, nrows = NULL) {
+    .Call(`_pc_RcppPCboot`, target, source, libsizes, lib, pred, E, tau, style, num_neighbors, zero_tolerance, dist_metric, boot, replace_sampling, seed, relative, weighted, threads, parallel_level, verbose, h, nb, nrows)
 }
 
 RcppPCops <- function(target, source, lib, pred, E, tau, k, maximize = "positive", style = 0L, zero_tolerance = 0L, dist_metric = "euclidean", relative = TRUE, weighted = TRUE, threads = 1L, parallel_level = 0L, h = 0L, nb = NULL, nrows = NULL) {

R/pc.R

@@ -1,5 +1,5 @@
 .pc_ts = \(data, target, source, libsizes = NULL, E = 3, k = E, tau = 1, style = 1, lib = NULL, pred = NULL, boot = 99,
-           random = TRUE, seed = 42L, dist.metric = c("euclidean", "manhattan", "maximum"), zero.tolerance = max(k),
+           replace = TRUE, seed = 42L, dist.metric = c("euclidean", "manhattan", "maximum"), zero.tolerance = max(k),
            relative = TRUE, weighted = TRUE, threads = length(libsizes), higher.parallel = TRUE, verbose = TRUE, h = 0, ...) {
   dist.metric = match.arg(dist.metric)
   dlist = .validate_var(data, target, source)
@@ -12,12 +12,12 @@
                   dist.metric, relative, weighted, threads, h, NULL, NULL))
   } else {
     return(RcppPCboot(tv, sv, libsizes, lib, pred, E, tau, style, k, zero.tolerance, dist.metric, boot,
-                      random, seed, relative, weighted, threads, higher.parallel, verbose, h, NULL, NULL))
+                      replace, seed, relative, weighted, threads, higher.parallel, verbose, h, NULL, NULL))
   }
 }
 
 .pc_lattice = \(data, target, source, libsizes = NULL, E = 3, k = E+1, tau = 1, style = 1, lib = NULL, pred = NULL, boot = 99,
-                random = TRUE, seed = 42L, dist.metric = c("euclidean", "manhattan", "maximum"), zero.tolerance = max(k), 
+                replace = TRUE, seed = 42L, dist.metric = c("euclidean", "manhattan", "maximum"), zero.tolerance = max(k), 
                 relative = TRUE, weighted = TRUE, threads = length(libsizes), higher.parallel = TRUE, verbose = TRUE, detrend = FALSE, nb = NULL, ...) {
   if (is.null(nb)) nb = sdsfun::spdep_nb(data)
   dist.metric = match.arg(dist.metric)
@@ -31,12 +31,12 @@
                   dist.metric, relative, weighted, threads, 0, nb, NULL))
   } else {
     return(RcppPCboot(tv, sv, libsizes, lib, pred, E, tau, style, k, zero.tolerance, dist.metric, boot,
-                      random, seed, relative, weighted, threads, higher.parallel, verbose, 0, nb, NULL))
+                      replace, seed, relative, weighted, threads, higher.parallel, verbose, 0, nb, NULL))
   }
 }
 
 .pc_grid = \(data, target, source, libsizes = NULL, E = 3, k = E+1, tau = 1, style = 1, lib = NULL, pred = NULL, boot = 99,
-             random = TRUE, seed = 42L, dist.metric = c("euclidean", "manhattan", "maximum"), zero.tolerance = max(k), 
+             replace = TRUE, seed = 42L, dist.metric = c("euclidean", "manhattan", "maximum"), zero.tolerance = max(k), 
              relative = TRUE, weighted = TRUE, threads = length(libsizes), higher.parallel = TRUE, verbose = TRUE, detrend = FALSE, ...) {
   dist.metric = match.arg(dist.metric)
   dlist = .validate_var(data, target, source, detrend)
@@ -48,7 +48,7 @@
     return(RcppPC(tv, sv, lib, pred, E, tau, style, k, zero.tolerance, dist.metric,
                   relative, weighted, threads, 0, NULL, terra::nrow(data)))
   } else {
-    return(RcppPCboot(tv, sv, libsizes, lib, pred, E, tau, style, k, zero.tolerance, dist.metric, boot, random,
+    return(RcppPCboot(tv, sv, libsizes, lib, pred, E, tau, style, k, zero.tolerance, dist.metric, boot, replace,
                       seed, relative, weighted, threads, higher.parallel, verbose, 0, NULL, terra::nrow(data)))
   }
 }
@@ -60,7 +60,7 @@
 #' @param libsizes (optional) Number of observations used.
 #' @param k (optional) Number of nearest neighbors used for projection.
 #' @param boot (optional) Number of bootstraps to perform.
-#' @param random (optional) Whether to use random sampling.
+#' @param replace (optional) Should sampling be with replacement?
 #' @param seed (optional) Random seed.
 #' @param zero.tolerance (optional) Maximum number of zeros tolerated in signature space.
 #' @param relative (optional) Whether to calculate relative changes in embedding.

inst/include/pc/patcaus.hpp

@@ -224,7 +224,7 @@ namespace patcaus
         const size_t& h = 0,
         const std::string& dist_metric = "euclidean",
         size_t boot = 99,
-        bool random_sample = true,
+        bool replace_sampling = true,
         unsigned long long seed = 42,
         bool relative = true,
         bool weighted = true,
@@ -237,9 +237,6 @@ namespace patcaus
     // --------------------------------------------------------------------------
     threads = std::min(static_cast<size_t>(std::thread::hardware_concurrency()), threads);
 
-    // Enforce boot = 1 for deterministic sampling
-    if (!random_sample) boot = 1;
-
     // Prebuild 64-bit RNG pool for reproducibility
     std::vector<std::mt19937_64> rng_pool(boot);
     for (size_t i = 0; i < boot; ++i) {
@@ -305,17 +302,31 @@ namespace patcaus
 
         auto process_boot = [&](size_t b) 
         {
-            std::vector<size_t> sampled_lib, sampled_pred;
+            std::vector<size_t> sampled_lib;
+            // std::vector<size_t> sampled_pred;
+
+            if (boot == 1)
+            {
+                sampled_lib.assign(lib_indices.begin(), lib_indices.begin() + L);
+                // sampled_pred = sampled_lib;
+            }
+            else if (replace_sampling) 
+            {   
+                sampled_lib.resize(L);
+                std::uniform_int_distribution<size_t> dist(0, lib_indices.size() - 1);
 
-            if (random_sample) 
+                for (size_t i = 0; i < L; ++i)
+                {
+                    sampled_lib[i] = lib_indices[dist(rng_pool[b])];
+                }
+                // sampled_pred = sampled_lib;
+            } 
+            else 
             {
                 std::vector<size_t> shuffled_lib = lib_indices;
                 std::shuffle(shuffled_lib.begin(), shuffled_lib.end(), rng_pool[b]);
                 sampled_lib.assign(shuffled_lib.begin(), shuffled_lib.begin() + L);
                 // sampled_pred = sampled_lib;
-            } else {
-                sampled_lib.assign(lib_indices.begin(), lib_indices.begin() + L);
-                // sampled_pred = sampled_lib;
             }
 
             std::vector<std::vector<double>> PredSMy;

src/PC.cpp

@@ -344,7 +344,7 @@ Rcpp::List RcppPCboot(
     int zero_tolerance = 0,
     const std::string& dist_metric = "euclidean",
     int boot = 99,
-    bool random_sample = true,
+    bool replace_sampling = true,
     int seed = 42,
     bool relative = true,
     bool weighted = true,
@@ -554,7 +554,7 @@ Rcpp::List RcppPCboot(
             static_cast<size_t>(std::abs(num_neighbors)),
             static_cast<size_t>(std::abs(zero_tolerance)),
             static_cast<size_t>(std::abs(h)), dist_metric, 
-            static_cast<size_t>(std::abs(boot)), random_sample, 
+            static_cast<size_t>(std::abs(boot)), replace_sampling, 
             static_cast<unsigned long long>(std::abs(seed)),
             relative, weighted, static_cast<size_t>(std::abs(threads)),
             static_cast<size_t>(std::abs(parallel_level)), verbose);
@@ -602,7 +602,7 @@ Rcpp::List RcppPCboot(
             static_cast<size_t>(std::abs(num_neighbors)),
             static_cast<size_t>(std::abs(zero_tolerance)),
             static_cast<size_t>(std::abs(h)), dist_metric, 
-            static_cast<size_t>(std::abs(boot)), random_sample, 
+            static_cast<size_t>(std::abs(boot)), replace_sampling, 
             static_cast<unsigned long long>(std::abs(seed)),
             relative, weighted, static_cast<size_t>(std::abs(threads)),
             static_cast<size_t>(std::abs(parallel_level)), verbose);

src/RcppExports.cpp

@@ -61,8 +61,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // RcppPCboot
-Rcpp::List RcppPCboot(const Rcpp::NumericVector& target, const Rcpp::NumericVector& source, const Rcpp::IntegerVector& libsizes, const Rcpp::IntegerVector& lib, const Rcpp::IntegerVector& pred, const Rcpp::IntegerVector& E, const Rcpp::IntegerVector& tau, int style, int num_neighbors, int zero_tolerance, const std::string& dist_metric, int boot, bool random_sample, int seed, bool relative, bool weighted, int threads, int parallel_level, bool verbose, int h, Rcpp::Nullable<Rcpp::List> nb, Rcpp::Nullable<int> nrows);
-RcppExport SEXP _pc_RcppPCboot(SEXP targetSEXP, SEXP sourceSEXP, SEXP libsizesSEXP, SEXP libSEXP, SEXP predSEXP, SEXP ESEXP, SEXP tauSEXP, SEXP styleSEXP, SEXP num_neighborsSEXP, SEXP zero_toleranceSEXP, SEXP dist_metricSEXP, SEXP bootSEXP, SEXP random_sampleSEXP, SEXP seedSEXP, SEXP relativeSEXP, SEXP weightedSEXP, SEXP threadsSEXP, SEXP parallel_levelSEXP, SEXP verboseSEXP, SEXP hSEXP, SEXP nbSEXP, SEXP nrowsSEXP) {
+Rcpp::List RcppPCboot(const Rcpp::NumericVector& target, const Rcpp::NumericVector& source, const Rcpp::IntegerVector& libsizes, const Rcpp::IntegerVector& lib, const Rcpp::IntegerVector& pred, const Rcpp::IntegerVector& E, const Rcpp::IntegerVector& tau, int style, int num_neighbors, int zero_tolerance, const std::string& dist_metric, int boot, bool replace_sampling, int seed, bool relative, bool weighted, int threads, int parallel_level, bool verbose, int h, Rcpp::Nullable<Rcpp::List> nb, Rcpp::Nullable<int> nrows);
+RcppExport SEXP _pc_RcppPCboot(SEXP targetSEXP, SEXP sourceSEXP, SEXP libsizesSEXP, SEXP libSEXP, SEXP predSEXP, SEXP ESEXP, SEXP tauSEXP, SEXP styleSEXP, SEXP num_neighborsSEXP, SEXP zero_toleranceSEXP, SEXP dist_metricSEXP, SEXP bootSEXP, SEXP replace_samplingSEXP, SEXP seedSEXP, SEXP relativeSEXP, SEXP weightedSEXP, SEXP threadsSEXP, SEXP parallel_levelSEXP, SEXP verboseSEXP, SEXP hSEXP, SEXP nbSEXP, SEXP nrowsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type target(targetSEXP);
@@ -77,7 +77,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type zero_tolerance(zero_toleranceSEXP);
     Rcpp::traits::input_parameter< const std::string& >::type dist_metric(dist_metricSEXP);
     Rcpp::traits::input_parameter< int >::type boot(bootSEXP);
-    Rcpp::traits::input_parameter< bool >::type random_sample(random_sampleSEXP);
+    Rcpp::traits::input_parameter< bool >::type replace_sampling(replace_samplingSEXP);
     Rcpp::traits::input_parameter< int >::type seed(seedSEXP);
     Rcpp::traits::input_parameter< bool >::type relative(relativeSEXP);
     Rcpp::traits::input_parameter< bool >::type weighted(weightedSEXP);
@@ -87,7 +87,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type h(hSEXP);
     Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::List> >::type nb(nbSEXP);
     Rcpp::traits::input_parameter< Rcpp::Nullable<int> >::type nrows(nrowsSEXP);
-    rcpp_result_gen = Rcpp::wrap(RcppPCboot(target, source, libsizes, lib, pred, E, tau, style, num_neighbors, zero_tolerance, dist_metric, boot, random_sample, seed, relative, weighted, threads, parallel_level, verbose, h, nb, nrows));
+    rcpp_result_gen = Rcpp::wrap(RcppPCboot(target, source, libsizes, lib, pred, E, tau, style, num_neighbors, zero_tolerance, dist_metric, boot, replace_sampling, seed, relative, weighted, threads, parallel_level, verbose, h, nb, nrows));
     return rcpp_result_gen;
 END_RCPP
 }