Update benchmark: two-panel script with warm-up and column-count sweep

Previous single-panel benchmark was misleading: the first sequential run
paid one-time cache-warm-up cost (skewing n=20), and only 6 columns were
used so only 6 cores were ever active on a 16-core machine.

New two-panel design:
  Left  – vary rows with 16 fixed AB columns (shows memory-bandwidth
          saturation for large n)
  Right – vary columns with fixed rows (shows the real speedup profile:
          parallel wins when n_cols >> 1)

Also adds a warm-up pass before measurements to eliminate first-call bias.

https://claude.ai/code/session_012DXCXbZUC54Zij1z9bFiHR

tools/benchmark_parallel.R

@@ -1,71 +1,115 @@
-# Benchmark: sequential vs parallel as.sir() across data-set sizes
+# Benchmark: sequential vs parallel as.sir() across data-set shapes
 #
-# Run from the repo root with:
+# Run from the repo root:
 #   Rscript tools/benchmark_parallel.R
-# or from inside an R session:
+# or inside an R session:
 #   source("tools/benchmark_parallel.R")
 #
-# Requires ggplot2 for the output plot; uses devtools::load_all() so the
-# package does not need to be installed.
+# Two panels:
+#   Left  – fixed columns (n_ab_fixed), varying rows.
+#           Parallel wins at small n; sequential catches up at large n due to
+#           memory-bandwidth saturation (all workers compete for the same
+#           clinical_breakpoints lookup table in L3 cache / RAM).
+#   Right – fixed rows (n_rows_fixed), varying column count.
+#           This is the shape that actually benefits: each additional column
+#           keeps another core busy.  The "real world" gain for a 2854×65
+#           dataset lives here.
+#
+# Requires ggplot2; uses devtools::load_all() so the package need not be
+# installed.
 
 devtools::load_all(".", quiet = TRUE)
 
-sizes <- c(20, 200, 2000, 20000)
-n_ab  <- 6   # number of antibiotic columns
+# ── configuration ─────────────────────────────────────────────────────────────
+row_sizes   <- c(200, 1000, 5000, 20000)
+col_sizes   <- c(4, 8, 16, 32, 48)
+n_rows_fixed <- 1000
+n_ab_fixed   <- 16
+n_cores_avail <- AMR:::get_n_cores(Inf)
 
-make_df <- function(n) {
+all_abs <- c("AMC", "GEN", "CIP", "TZP", "IPM", "MEM",
+             "AMP", "TMP", "SXT", "NIT", "FOX", "CRO",
+             "FEP", "CAZ", "CTX", "TOB", "AMK", "ERY",
+             "AZM", "CLI", "VAN", "TEC", "RIF", "MTR",
+             "MFX", "LNZ", "TGC", "DOX", "FLC", "OXA",
+             "PEN", "CXM", "CZO", "KAN", "COL", "FOS",
+             "MUP", "TCY", "TEC", "IPM", "CHL", "FEP",
+             "MEM", "TZP", "GEN", "AMC", "AMX", "AMP")
+all_abs <- unique(all_abs)
+
+mic_vals <- c("0.25", "0.5", "1", "2", "4", "8", "16", "32")
+
+make_df <- function(n_rows, n_ab) {
   set.seed(42)
-  mics <- lapply(seq_len(n_ab), function(j) {
-    as.mic(sample(c("0.25", "0.5", "1", "2", "4", "8", "16", "32"), n, TRUE))
-  })
-  names(mics) <- c("AMC", "GEN", "CIP", "TZP", "IPM", "MEM")
+  ab_sel <- all_abs[seq_len(min(n_ab, length(all_abs)))]
+  mics <- lapply(ab_sel, function(a) as.mic(sample(mic_vals, n_rows, TRUE)))
+  names(mics) <- ab_sel
   data.frame(mo = "B_ESCHR_COLI", mics, stringsAsFactors = FALSE)
 }
 
-results <- do.call(rbind, lapply(sizes, function(n) {
-  df <- make_df(n)
-
+time_both <- function(n_rows, n_ab, label) {
+  df <- make_df(n_rows, n_ab)
   t_seq <- system.time(
     suppressMessages(as.sir(df, col_mo = "mo", info = FALSE, parallel = FALSE))
   )[["elapsed"]]
-
   t_par <- system.time(
     suppressMessages(as.sir(df, col_mo = "mo", info = FALSE, parallel = TRUE))
   )[["elapsed"]]
+  message(sprintf("%-28s  seq=%5.2fs  par=%5.2fs  speedup=%.1fx",
+                  label, t_seq, t_par, t_seq / t_par))
+  data.frame(group = label, mode = c("sequential", "parallel"),
+             seconds = c(t_seq, t_par), stringsAsFactors = FALSE)
+}
 
-  message(sprintf("n = %6d  seq = %.3fs  par = %.3fs  speedup = %.1fx",
-                  n, t_seq, t_par, t_seq / t_par))
+# ── warm-up (avoid first-call overhead biasing results) ───────────────────────
+message("Warming up cache ...")
+invisible(suppressMessages(as.sir(make_df(100, 6), col_mo = "mo", info = FALSE)))
+invisible(suppressMessages(as.sir(make_df(100, 6), col_mo = "mo", info = FALSE, parallel = TRUE)))
+sir_interpretation_history(clean = TRUE)
 
-  data.frame(n = n, mode = c("sequential", "parallel"),
-             seconds = c(t_seq, t_par))
+# ── panel 1: vary rows, fixed columns ─────────────────────────────────────────
+message(sprintf("\nPanel 1 – varying rows, %d fixed columns:", n_ab_fixed))
+res_rows <- do.call(rbind, lapply(row_sizes, function(n) {
+  time_both(n, n_ab_fixed, sprintf("rows=%d", n))
 }))
+res_rows$x     <- rep(row_sizes, each = 2)
+res_rows$panel <- "Vary rows (16 fixed AB columns)"
+
+# ── panel 2: vary columns, fixed rows ─────────────────────────────────────────
+message(sprintf("\nPanel 2 – varying columns, %d fixed rows:", n_rows_fixed))
+res_cols <- do.call(rbind, lapply(col_sizes, function(n_ab) {
+  time_both(n_rows_fixed, n_ab, sprintf("cols=%d", n_ab))
+}))
+res_cols$x     <- rep(col_sizes, each = 2)
+res_cols$panel <- sprintf("Vary columns (%d fixed rows)", n_rows_fixed)
+
+results <- rbind(res_rows, res_cols)
 
 if (requireNamespace("ggplot2", quietly = TRUE)) {
-  p <- ggplot2::ggplot(results, ggplot2::aes(x = n, y = seconds,
-                                             colour = mode, group = mode)) +
+  p <- ggplot2::ggplot(
+    results,
+    ggplot2::aes(x = x, y = seconds, colour = mode, group = mode)
+  ) +
     ggplot2::geom_line(linewidth = 1) +
-    ggplot2::geom_point(size = 3) +
-    ggplot2::scale_x_log10(
-      breaks = sizes,
-      labels = format(sizes, big.mark = ",", scientific = FALSE)
-    ) +
+    ggplot2::geom_point(size = 2.5) +
+    ggplot2::facet_wrap(~panel, scales = "free_x") +
     ggplot2::scale_colour_manual(
       values = c(sequential = "#E05C5C", parallel = "#2E86AB")
     ) +
     ggplot2::labs(
       title    = "as.sir() throughput: sequential vs parallel",
-      subtitle = sprintf("%d antibiotic columns, E. coli, EUCAST 2025", n_ab),
-      x        = "Number of rows (log scale)",
+      subtitle = sprintf("E. coli, EUCAST 2026, %d cores available", n_cores_avail),
+      x        = "Dataset dimension (rows  ·left·  or columns  ·right·)",
       y        = "Wall-clock time (seconds)",
       colour   = NULL
     ) +
-    ggplot2::theme_minimal(base_size = 13) +
+    ggplot2::theme_minimal(base_size = 12) +
     ggplot2::theme(legend.position = "top")
 
   out_file <- "tools/benchmark_parallel.png"
-  ggplot2::ggsave(out_file, p, width = 7, height = 5, dpi = 150)
-  message("Plot saved to ", out_file)
+  ggplot2::ggsave(out_file, p, width = 10, height = 5, dpi = 150)
+  message("\nPlot saved to ", out_file)
 } else {
   message("Install ggplot2 to get a plot; raw results:")
-  print(results)
+  print(results[, c("panel", "group", "mode", "seconds")])
 }