rejection sampling tests

Author: pachadotdev

R/vendor.R

@@ -63,7 +63,7 @@ vendor <- function(path = "./src/vendor") {
       cpp4r_version,
       Sys.Date()
     ),
-    file.path(path2, "00-cpp4r-vendoring-info.txt")
+    file.path(path, "00-cpp4r-vendoring-info.txt")
   )
 
   file.copy(

Rcppbenchgibbs/R/Rcppbenchgibbs-package.R

@@ -2,14 +2,10 @@ Package: cpp4rexamples
 Type: Package
 Title: ADD TITLE
 Version: 0.1
-Authors@R: c(
-    person(
-        given = "YOUR",
-        family = "NAME",
-        role = c("aut", "cre"),
-        email = "[email protected]",
-        comment = c(ORCID = "0000-0001-0002-0003"))
-    )
+Authors@R:
+    person("Mauricio", "Vargas Sepulveda", role = c("aut","cre","cph"), 
+        email = "[email protected]",
+        comment = c(ORCID = "0000-0003-1017-7574")
 Suggests: 
     knitr,
     rmarkdown

Rcppbenchgibbs/R/base.R

@@ -30,19 +30,32 @@ class local_rng {
   // Manage R's RNG state
   local_rng rng_state;
 
-  // Dynamic vector that grows as we accept samples
+  // Pre-calculate acceptance rate for better initial allocation
+  double z_lower = (lower - mu) / sigma;
+  double z_upper = (upper - mu) / sigma;
+  double acceptance_rate = Rf_pnorm5(z_upper, 0.0, 1.0, 1, 0) - Rf_pnorm5(z_lower, 0.0, 1.0, 1, 0);
+  
+  // Pre-allocate based on expected number of samples needed
+  // Add 20% buffer for variance + ensure minimum reasonable size
+  R_xlen_t estimated_needed = static_cast<R_xlen_t>(n_samples / acceptance_rate * 1.2);
+  estimated_needed = std::max(estimated_needed, static_cast<R_xlen_t>(n_samples));
+  
   writable::doubles accepted_samples;
+  accepted_samples.reserve(estimated_needed);
 
   // Keep sampling until we have enough accepted samples
-  while (accepted_samples.size() < n_samples) {
+  // Cast to int once for faster comparison in tight loop
+  int target_samples = static_cast<int>(n_samples);
+  
+  while (static_cast<int>(accepted_samples.size()) < target_samples) {
     // Generate candidate from normal distribution
     double candidate = Rf_rnorm(mu, sigma);
 
-    // Accept if within bounds
-    if (candidate >= lower && candidate <= upper) {
+    // Fast bounds check: single branch for most common case
+    if (__builtin_expect(candidate >= lower && candidate <= upper, 1)) {
       accepted_samples.push_back(candidate);
     }
-    // If rejected, we just continue sampling - this is why we need dynamic growth!
+    // If rejected, we just continue sampling
   }
 
   return accepted_samples;

Rcppbenchgibbs/src/main.cpp

@@ -1,39 +1,29 @@
-library(cpp11test)
-library(cpp4rtest)
+library(cpp11sampling)
+library(cpp4rsampling)
+library(Rcppsampling)
+library(bench)
+library(purrr)
 
-cases <- expand.grid(
-  len = 1e7,
-  vector = c("normal", "altrep"),
-  method = c("for", "foreach", "accumulate"),
-  pkg = c("cpp4r", "Rcpp"),
-  stringsAsFactors = FALSE
-)
-
-# Add special case
-cases <- rbind(list(len = 1e7, vector = "normal", method = "for2", pkg = "cpp4r"), cases)
+rejection_sizes <- seq(25000, 100000, by = 25000)
 
-b_sum <- bench::press(
-  .grid = cases,
-  {
-    seq_real <- function(x) as.numeric(seq_len(x))
-    funs <- c("normal" = rnorm, "altrep" = seq_real)
-    x <- funs[[vector]](len)
-    fun <- match.fun(sprintf("%ssum_dbl_%s_", ifelse(pkg == "cpp4r", "", paste0(pkg, "_")), method))
-    bench::mark(
-      fun(x)
+rejection_bench <- map_df(
+  rejection_sizes,
+  function(n_samples) {
+    cat("Target samples:", n_samples, "\n")
+    d <- bench::mark(
+      # "R" = { set.seed(123); cpp4rsampling::rejection_sampling(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5) },
+      "cpp4r" = { set.seed(123); cpp4rsampling::rejection_sampling_cpp4r(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5) },
+      "cpp11" = { set.seed(123); cpp11sampling::rejection_sampling_cpp11(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5) },
+      "Rcpp" = { set.seed(123); Rcppsampling::rejection_sampling_Rcpp(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5) },
+      iterations = 10
     )
-  }
-)[c("pkg", "method", "vector", "median", "mem_alloc")]
 
-saveRDS(b_sum, "sum.Rds", version = 2)
+    d$len <- n_samples
+
+    d
+  }
+)
 
-grid <- expand.grid(len = seq(10000, 50000, by = 10000), pkg = c("cpp4r", "cpp11", "Rcpp"), stringsAsFactors = FALSE)
-b_release <- bench::press(.grid = grid, {
-  fun <- match.fun(sprintf("%s_release_", pkg))
-  bench::mark(
-    fun(len),
-    min_iterations = 100
-  )
-})[c("len", "pkg", "median")]
+rejection_bench[, c("expression", "median", "mem_alloc", "len")]
 
-saveRDS(b_release, "release.Rds", version = 2)
+saveRDS(rejection_bench, "vignettes/rejection_bench.rds", compress = "xz")

cpp11benchgibbs/src/main.cpp

@@ -0,0 +1,67 @@
+library(cpp11benchgibbs)
+library(cpp4rbenchgibbs)
+library(Rcppbenchgibbs)
+library(bench)
+library(purrr)
+
+cat("=== Benchmarking Dynamic Vector Growth with push_back ===\n")
+cat("cpp4r/cpp11: Efficient - reserves extra space when growing vectors\n")
+cat("Rcpp: Inefficient - copies entire vector on each push_back operation\n\n")
+
+cat("=== 1. Rejection Sampling Performance ===\n")
+cat("Perfect use case: We don't know how many candidates we'll need to generate\n")
+cat("to get the desired number of accepted samples.\n\n")
+
+# Test rejection sampling with different target sample sizes
+rejection_sizes <- c(1000, 5000, 10000, 25000)
+
+rejection_bench <- map_df(
+  rejection_sizes,
+  function(n_samples) {
+    cat("Target samples:", n_samples, "\n")
+    bench::mark(
+      "cpp4r" = cpp4rbenchgibbs::rejection_sampling_cpp4r(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5),
+      "cpp11" = cpp11benchgibbs::rejection_sampling_cpp11(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5),
+      "Rcpp" = Rcppbenchgibbs::rejection_sampling_Rcpp(n_samples, mu = 0, sigma = 1, lower = -1.5, upper = 1.5),
+      min_iterations = 5,
+      max_iterations = 20,
+      check = FALSE  # Results will differ due to randomness
+    )
+  }
+)
+
+print(rejection_bench[c("expression", "median", "mem_alloc", "n_itr")])
+
+cat("\n=== 2. Bootstrap with Variable Sample Sizes ===\n")
+cat("Another natural use case: Bootstrap samples with varying sizes\n\n")
+
+# Create test data
+test_data <- rnorm(1000)
+
+bootstrap_bench <- bench::mark(
+  "cpp4r" = cpp4rbenchgibbs::bootstrap_variable_cpp4r(test_data, min_size = 50, max_size = 200, n_bootstrap = 1000),
+  "cpp11" = cpp11benchgibbs::bootstrap_variable_cpp11(test_data, min_size = 50, max_size = 200, n_bootstrap = 1000),
+  "Rcpp" = Rcppbenchgibbs::bootstrap_variable_Rcpp(test_data, min_size = 50, max_size = 200, n_bootstrap = 1000),
+  min_iterations = 3,
+  max_iterations = 10,
+  check = FALSE
+)
+
+print(bootstrap_bench[c("expression", "median", "mem_alloc", "n_itr")])
+
+cat("\n=== Key Performance Insights ===\n")
+cat("1. Rejection Sampling: Rcpp should show quadratic time complexity\n")
+cat("   due to vector copying, while cpp4r/cpp11 remain linear\n")
+cat("2. Memory Usage: Rcpp will allocate much more memory due to\n")
+cat("   repeated copying of growing vectors\n")
+cat("3. Real-world Impact: These patterns appear in Monte Carlo methods,\n")
+cat("   adaptive algorithms, and data processing pipelines\n\n")
+
+# Demonstrate the acceptance rate for rejection sampling
+cat("=== Rejection Sampling Statistics ===\n")
+set.seed(42)
+sample1 <- cpp4rbenchgibbs::rejection_sampling_cpp4r(1000, mu = 0, sigma = 1, lower = -1.5, upper = 1.5)
+acceptance_rate <- pnorm(1.5) - pnorm(-1.5)  # Theoretical acceptance rate
+cat("Theoretical acceptance rate for [-1.5, 1.5] bounds:", round(acceptance_rate, 3), "\n")
+cat("Average candidates needed per accepted sample:", round(1/acceptance_rate, 1), "\n")
+cat("For 25,000 samples, expect ~", round(25000/acceptance_rate), "candidates to be generated\n")