feat(app): add a Bootstrap tab and use the Ward rotavirus data as the example

app/DESCRIPTION

@@ -18,17 +18,19 @@ Depends:
 Imports:
     bslib,
     config,
-    esquisse (>= 2.0.0),
+    esquisse (>= 2.1.0),
     generics,
     ggplot2,
     golem,
+    mirai (>= 2.5.0),
+    parallel,
     plotly,
     readr,
     shiny,
-    singlehit (>= 0.2.0)
+    singlehit (>= 0.3.0)
 Suggests:
     pkgload,
-    testthat (>= 3.0.0),
+    testthat (>= 3.1.7),
     withr
 Remotes:
     seanthimons/singlehit

app/R/app_server.R

@@ -8,6 +8,9 @@ app_server <- function(input, output, session) {
   data <- mod_upload_server("upload")
   analysis <- mod_fit_server("fit", data)
   mod_compare_server("compare", analysis)
-  mod_figures_server("figures", analysis)
-  mod_explore_server("explore", analysis)
+  # Plots and Explore get the analysis with the bootstraps attached once a run
+  # finishes; Compare does not use them.
+  analysis_with_bands <- mod_bootstrap_server("bootstrap", analysis)
+  mod_figures_server("figures", analysis_with_bands)
+  mod_explore_server("explore", analysis_with_bands)
 }

app/R/app_ui.R

@@ -7,6 +7,10 @@
 app_ui <- function(request) {
   tagList(
     golem_add_external_resources(),
+    # A pulse bar across the top while the server is busy, and spinners on
+    # outputs that are recalculating. R is single threaded, so a long fit or a
+    # run starting can make the app unresponsive; this says why.
+    shiny::useBusyIndicators(),
     bslib::page_navbar(
       title = "singlehit",
       bslib::nav_panel(
@@ -18,6 +22,7 @@ app_ui <- function(request) {
         bslib::layout_sidebar(sidebar = bslib::sidebar(mod_fit_ui("fit")), mod_fit_results_ui("fit"))
       ),
       bslib::nav_panel("Compare", mod_compare_ui("compare")),
+      bslib::nav_panel("Bootstrap", mod_bootstrap_ui("bootstrap")),
       bslib::nav_panel("Plots", mod_figures_ui("figures")),
       bslib::nav_panel("Explore", mod_explore_ui("explore"))
     )

app/R/fct_bootstrap.R

@@ -0,0 +1,172 @@
+#' Replicate limits for the bootstrap controls
+#'
+#' The minimum is 100 because a 99% percentile interval from fewer draws is
+#' meaningless. The exact model is far slower per replicate, so it has a smaller
+#' ceiling: the shared pool is first-in first-out, and a very long exact run
+#' would hold up every other session.
+#'
+#' The ceiling matches singlehit's own default of 10,000. Starting a run prepares
+#' every replicate on the app's single R process, which blocked it for 17 s at
+#' 10,000 replicates (two models) before singlehit drew the replicates in one
+#' vectorized call and stopped shipping the whole task list with every batch; it
+#' is now about a second. See dev/benchmark_bootstrap.R.
+#' @noRd
+bootstrap_limits <- list(min = 100L, max = 10000L, exact_max = 1000L)
+
+#' Gentle defaults for the bootstrap controls
+#' @noRd
+bootstrap_defaults <- list(times = 1000L, exact_times = 200L)
+
+check_count <- function(x, label, lower, upper) {
+  if (length(x) != 1L || !is.numeric(x) || is.na(x) || x != floor(x) || x < lower || x > upper) {
+    stop(
+      sprintf(
+        "%s must be a whole number between %s and %s.",
+        label,
+        format(lower, big.mark = ","),
+        format(upper, big.mark = ",")
+      ),
+      call. = FALSE
+    )
+  }
+  as.integer(x)
+}
+
+#' Work out what to run for each fitted model
+#'
+#' Seeds follow `analyze_dose_response()`: `seed + index - 1` over the fitted
+#' models, so a run can be reproduced from the exported script.
+#'
+#' @param analysis A `qdr_analysis`.
+#' @param times Replicates for the exponential and beta-Poisson models.
+#' @param exact_times Replicates for the exact model; used only if it was fitted.
+#' @param seed Optional whole number; `NULL` or `NA` for none.
+#' @return A named list (by model) of `list(model, times, seed)`.
+#' @noRd
+bootstrap_plan <- function(analysis, times, exact_times = NULL, seed = NULL) {
+  models <- names(analysis$fits)
+  times <- check_count(times, "Replicates", bootstrap_limits$min, bootstrap_limits$max)
+  if ("exact_beta_poisson" %in% models) {
+    if (is.null(exact_times)) {
+      exact_times <- bootstrap_defaults$exact_times
+    }
+    exact_times <- check_count(exact_times, "Exact-model replicates", bootstrap_limits$min, bootstrap_limits$exact_max)
+  }
+  if (is.null(seed) || length(seed) != 1L || is.na(seed)) {
+    seed <- NULL
+  } else {
+    seed <- check_count(seed, "The seed", 0L, .Machine$integer.max - length(models))
+  }
+  plan <- lapply(seq_along(models), function(i) {
+    list(
+      model = models[[i]],
+      times = if (models[[i]] == "exact_beta_poisson") exact_times else times,
+      seed = if (is.null(seed)) NULL else seed + i - 1L
+    )
+  })
+  names(plan) <- models
+  plan
+}
+
+#' Summarize the progress of several jobs
+#'
+#' @param progress Named list of `singlehit::bootstrap_progress()` results.
+#' @return A list with a per-model data frame and the overall fraction.
+#' @noRd
+progress_summary <- function(progress) {
+  done <- vapply(progress, function(p) p$done, integer(1))
+  total <- vapply(progress, function(p) p$total, integer(1))
+  list(
+    models = data.frame(
+      model = names(progress),
+      done = unname(done),
+      total = unname(total),
+      fraction = ifelse(total > 0L, unname(done) / unname(total), 1),
+      stringsAsFactors = FALSE
+    ),
+    fraction = if (sum(total) > 0L) sum(done) / sum(total) else 1,
+    finished = all(done == total)
+  )
+}
+
+#' Split bootstraps into those that can draw bands and those that cannot
+#'
+#' `prediction_curve()` stops when fewer than two replicates converged, which
+#' would break the Plots tab, so such a bootstrap is reported and not attached.
+#'
+#' @param boots Named list of `qdr_bootstrap` objects.
+#' @return A list with `attached` (named list) and `dropped` (model names).
+#' @noRd
+attachable_bootstraps <- function(boots) {
+  ok <- vapply(boots, function(boot) sum(boot$converged, na.rm = TRUE) >= 2L, logical(1))
+  list(attached = boots[ok], dropped = names(boots)[!ok])
+}
+
+#' How many replicates converged, per model
+#' @noRd
+convergence_table <- function(boots) {
+  data.frame(
+    Model = model_label_ui(names(boots)),
+    `Converged replicates` = vapply(
+      boots,
+      function(boot) sprintf("%d of %d", sum(boot$converged, na.rm = TRUE), nrow(boot)),
+      character(1)
+    ),
+    check.names = FALSE,
+    row.names = NULL
+  )
+}
+
+#' Percentile intervals for one bootstrap, one row per quantity
+#'
+#' Reshapes `singlehit::bootstrap_confint()` (long, one row per quantity and
+#' level) into one row per quantity with a column per confidence level.
+#' @noRd
+bootstrap_interval_table <- function(boot) {
+  ci <- as.data.frame(singlehit::bootstrap_confint(boot))
+  quantities <- unique(ci$term)
+  labels <- c(n50 = "N50", ed10 = "ED10", ed50 = "ED50")
+  shown <- ifelse(quantities %in% names(labels), labels[quantities], quantities)
+  out <- data.frame(
+    Model = model_label_ui(ci$model[[1L]]),
+    Quantity = unname(shown),
+    stringsAsFactors = FALSE
+  )
+  for (level in sort(unique(ci$level))) {
+    part <- ci[ci$level == level, , drop = FALSE]
+    idx <- match(quantities, part$term)
+    out[[paste0(format(100 * level, trim = TRUE), "% interval")]] <- paste0(
+      "(",
+      format_number(part$lower[idx]),
+      ", ",
+      format_number(part$upper[idx]),
+      ")"
+    )
+  }
+  out
+}
+
+#' Interval tables for all models, stacked
+#' @noRd
+bootstrap_tables <- function(boots) {
+  do.call(rbind, unname(lapply(boots, bootstrap_interval_table)))
+}
+
+#' Attach finished bootstraps to the analysis they were computed from
+#'
+#' Bands must never describe a different fit than the one on screen, so the
+#' bootstraps are attached only when `base` is identical to the analysis they
+#' were run against.
+#'
+#' @param base The current `qdr_analysis`, or `NULL`.
+#' @param boots Named list of finished `qdr_bootstrap` objects (possibly empty).
+#' @param bootstrapped_base The analysis the bootstraps were run against.
+#' @return `base` with the bootstraps attached, or `base` unchanged.
+#' @noRd
+attach_bands <- function(base, boots, bootstrapped_base) {
+  if (!is.null(base) && length(boots) > 0L && identical(bootstrapped_base, base)) {
+    singlehit::add_bootstrap(base, boots)
+  } else {
+    base
+  }
+}

app/R/fct_explore.R

@@ -57,3 +57,31 @@ explore_datasets <- function(analysis) {
     model_residuals = model_residuals
   )
 }
+
+#' Input ids inside esquisse's module, relative to the Explore module
+#' @noRd
+esquisse_ids <- list(
+  aesthetics = "esquisse-geomaes-aes_1-dragvars",
+  geom = "esquisse-geomaes-geom_1",
+  x_transform = "esquisse-controls-axes-transX"
+)
+
+#' What the "Draw confidence bands" button fills in on the Fitted curves data
+#' @noRd
+bands_preset <- list(
+  aesthetics = list(xvar = "dose", yvar = "estimate", ymin = "lower", ymax = "upper", fill = "level", facet = "model"),
+  geom = "ribbon",
+  x_transform = "log10"
+)
+
+#' Does esquisse's reported drag-and-drop mapping contain the given variables?
+#'
+#' @param target The `target` element of the dragula input value: a list keyed
+#'   by slot name.
+#' @param wanted Named list: slot name to variable name.
+#' @noRd
+mapping_matches <- function(target, wanted) {
+  all(vapply(names(wanted), function(slot) {
+    identical(as.character(unlist(target[[slot]])), wanted[[slot]])
+  }, logical(1)))
+}

app/R/fct_figures.R

@@ -51,3 +51,16 @@ bounded_number <- function(x, lower, upper, default) {
   }
   min(max(x, lower), upper)
 }
+
+#' The view to switch to when confidence bands become available
+#'
+#' The overlay never draws bands, so someone looking at it right after a
+#' bootstrap would see no change and conclude the bands did not apply.
+#'
+#' @param type The currently selected view.
+#' @param has_bands Whether the analysis now carries bootstraps.
+#' @return The view to switch to, or `NULL` to leave the choice alone.
+#' @noRd
+view_for_bands <- function(type, has_bands) {
+  if (isTRUE(has_bands) && identical(type, "overlay")) "panels" else NULL
+}

app/R/fct_workers.R

@@ -0,0 +1,54 @@
+#' mirai compute profile shared by every bootstrap run in the app
+#'
+#' One named pool serves all sessions and all models. Without it each async
+#' bootstrap job would start its own pool sized to 75% of the cores, and only
+#' `collect_bootstrap()` shuts a job's own pool down.
+#' @noRd
+bootstrap_compute_profile <- "singlehitapp"
+
+#' Default number of bootstrap workers
+#'
+#' 75% of the physical cores, capped at 8 so a hosted deployment on a many-core
+#' machine does not claim all of its memory. Override with
+#' `options(singlehitapp.workers = n)`.
+#'
+#' @param cores Detected core count; injectable for tests.
+#' @return A whole number of at least 1.
+#' @noRd
+default_bootstrap_workers <- function(cores = parallel::detectCores(logical = FALSE)) {
+  chosen <- getOption("singlehitapp.workers")
+  if (!is.null(chosen)) {
+    return(max(1L, as.integer(chosen)))
+  }
+  if (length(cores) != 1L || is.na(cores) || cores < 1L) {
+    cores <- parallel::detectCores(logical = TRUE)
+  }
+  if (length(cores) != 1L || is.na(cores) || cores < 1L) {
+    return(1L)
+  }
+  max(1L, min(8L, as.integer(floor(0.75 * cores))))
+}
+
+#' Start the shared bootstrap workers if they are not running
+#'
+#' Safe to call before every run. Starting the pool ourselves means no job
+#' "owns" it, so one session finishing a job cannot shut the pool down under
+#' another session's job.
+#'
+#' @param workers Number of daemons to start.
+#' @noRd
+start_bootstrap_workers <- function(workers = default_bootstrap_workers()) {
+  if (!mirai::daemons_set(.compute = bootstrap_compute_profile)) {
+    mirai::daemons(workers, .compute = bootstrap_compute_profile)
+  }
+  invisible(TRUE)
+}
+
+#' Stop the shared bootstrap workers if they are running
+#' @noRd
+stop_bootstrap_workers <- function() {
+  if (mirai::daemons_set(.compute = bootstrap_compute_profile)) {
+    mirai::daemons(0L, .compute = bootstrap_compute_profile)
+  }
+  invisible(TRUE)
+}