From c79a0618d1c1f2ed89ad3b02b72e07fa0f1e54ac Mon Sep 17 00:00:00 2001
From: seanthimons <thimons.sean@gmail.com>
Date: Tue, 20 Oct 2026 12:29:31 -0400
Subject: [PATCH] perf(bootstrap): draw all replicate samples in one vectorized
 call

generate_bootstrap_tasks() built one tibble per replicate (about 0.7 ms each),
which took 7.8 s for 10,000 replicates and blocked the calling R process; in a
Shiny app that freezes every session while a run starts.

Draw every replicate in a single rbinom() call and build each sample with
tibble::new_tibble(). rbinom() consumes the random stream in order, so the
output is identical to the old per-replicate draws for any seed, and existing
seeds keep producing the same replicates.

Add tests that compare against the old algorithm exactly, for observed and
fitted probabilities, all three models, and a single replicate.
---
 R/bootstrap.R                         | 18 ++++--
 tests/testthat/test-bootstrap-tasks.R | 82 +++++++++++++++++++++++++++
 2 files changed, 95 insertions(+), 5 deletions(-)
 create mode 100644 tests/testthat/test-bootstrap-tasks.R

diff --git a/R/bootstrap.R b/R/bootstrap.R
index 67c39d0..b2b0025 100644
--- a/R/bootstrap.R
+++ b/R/bootstrap.R
@@ -170,14 +170,22 @@ new_qdr_bootstrap <- function(estimates, object, resample, backend) {
 
 generate_bootstrap_tasks <- function(object, probability, times, seed) {
   generate <- function() {
+    groups <- nrow(object$data)
+    total <- object$data$total
+    # One vectorized draw for every replicate. rbinom() consumes the random
+    # stream in order, so this is identical to drawing each replicate in turn,
+    # and far cheaper than building a tibble per replicate with tibble().
+    draws <- matrix(
+      stats::rbinom(groups * times, size = rep(total, times), prob = rep(probability, times)),
+      nrow = groups
+    )
     lapply(seq_len(times), function(replicate) {
-      positive <- stats::rbinom(nrow(object$data), size = object$data$total, prob = probability)
+      positive <- draws[, replicate]
       list(
         replicate = replicate,
-        sample_data = tibble::tibble(
-          dose = object$data$dose,
-          positive = positive,
-          negative = object$data$total - positive
+        sample_data = tibble::new_tibble(
+          list(dose = object$data$dose, positive = positive, negative = total - positive),
+          nrow = groups
         )
       )
     })
diff --git a/tests/testthat/test-bootstrap-tasks.R b/tests/testthat/test-bootstrap-tasks.R
new file mode 100644
index 0000000..41a8b50
--- /dev/null
+++ b/tests/testthat/test-bootstrap-tasks.R
@@ -0,0 +1,82 @@
+# The one-draw-per-replicate algorithm that generate_bootstrap_tasks() replaced.
+# The vectorized version must reproduce it exactly, so existing seeds keep giving
+# the same replicates.
+reference_tasks <- function(object, probability, times, seed) {
+  generate <- function() {
+    lapply(seq_len(times), function(replicate) {
+      positive <- stats::rbinom(nrow(object$data), size = object$data$total, prob = probability)
+      list(
+        replicate = replicate,
+        sample_data = tibble::tibble(
+          dose = object$data$dose,
+          positive = positive,
+          negative = object$data$total - positive
+        )
+      )
+    })
+  }
+  if (is.null(seed)) generate() else withr::with_seed(seed, generate())
+}
+
+ward_fit <- function(model = "beta_poisson") {
+  data <- read_dose_response(system.file("extdata", "Ward_rotavirus.txt", package = "singlehit"))
+  fit_dose_response(data, model)
+}
+
+test_that("vectorized task generation reproduces the per-replicate algorithm exactly", {
+  fit <- ward_fit()
+  for (probability in list(fit$data$response, fit$fitted)) {
+    expect_identical(
+      generate_bootstrap_tasks(fit, probability, 60L, seed = 7L),
+      reference_tasks(fit, probability, 60L, seed = 7L)
+    )
+  }
+})
+
+test_that("every model type and a single replicate produce the same tasks", {
+  for (model in c("exponential", "exact_beta_poisson")) {
+    fit <- ward_fit(model)
+    expect_identical(
+      generate_bootstrap_tasks(fit, fit$fitted, 25L, seed = 3L),
+      reference_tasks(fit, fit$fitted, 25L, seed = 3L)
+    )
+  }
+  fit <- ward_fit()
+  expect_identical(
+    generate_bootstrap_tasks(fit, fit$data$response, 1L, seed = 11L),
+    reference_tasks(fit, fit$data$response, 1L, seed = 11L)
+  )
+})
+
+test_that("tasks are numbered, well formed samples of the original groups", {
+  fit <- ward_fit()
+  tasks <- generate_bootstrap_tasks(fit, fit$data$response, 40L, seed = NULL)
+  expect_length(tasks, 40L)
+  expect_equal(vapply(tasks, function(task) task$replicate, integer(1)), 1:40)
+  for (task in tasks) {
+    sample <- task$sample_data
+    expect_s3_class(sample, "tbl_df")
+    expect_equal(sample$dose, fit$data$dose)
+    expect_equal(sample$positive + sample$negative, fit$data$total)
+    expect_true(all(sample$positive >= 0 & sample$negative >= 0))
+  }
+})
+
+test_that("a seed makes the tasks reproducible and leaves the global random stream alone", {
+  fit <- ward_fit()
+  withr::local_seed(99)
+  before <- .Random.seed
+  first <- generate_bootstrap_tasks(fit, fit$data$response, 30L, seed = 5L)
+  expect_identical(.Random.seed, before)
+  expect_identical(first, generate_bootstrap_tasks(fit, fit$data$response, 30L, seed = 5L))
+  expect_false(identical(first, generate_bootstrap_tasks(fit, fit$data$response, 30L, seed = 6L)))
+})
+
+test_that("bootstrap replicates are unchanged end to end", {
+  fit <- ward_fit()
+  tasks <- reference_tasks(fit, fit$data$response, 12L, seed = 21L)
+  expected <- dplyr::bind_rows(lapply(tasks, function(task) bootstrap_refit(fit, task$sample_data, task$replicate)))
+  boot <- bootstrap_dose_response(fit, times = 12L, seed = 21L, backend = "sequential")
+  expect_equal(boot$deviance, expected$deviance)
+  expect_equal(boot$converged, expected$converged)
+})