# =============================================================================
# test_doubleWishartHD_sweep.R
#
# FAST cross-parameter sweep harness for rootWishartHD upper-tail logSF.
#
# Purpose
#   Quickly probe correctness + timing across many (p, m_df, q_df) settings,
#   scaling up to p = 500, WITHOUT paying the Monte Carlo cost more than once.
#
# Design
#   * Monte Carlo theta draws are simulated ONCE per (p, m_df, q_df, nrep, seed)
#     and cached to disk (RDS). Re-running the sweep reuses the cached draws, so
#     only the (cheap) evaluation cost is repeated.
#   * For each setting we report three separated timings:
#       (A) MC simulation        -- paid once, then cached
#       (B) fast double CDF pass  -- whole vector, milliseconds
#       (C) exact multiprecision logSF      -- on a few upper-tail probe points (the slow
#                                    part); per-point timing is printed.
#   * n_exact is kept small so a full sweep up to p=500 stays tractable. The
#     exact cost grows ~ O(s^3) in the Pfaffian for large s, so the largest
#     settings use n_exact = 1.
#
# Run from the repository root:
#   Rscript test_doubleWishartHD_sweep.R              # full grid
#   Rscript test_doubleWishartHD_sweep.R p150 p500    # selected settings
#   N_EXACT=2 Rscript test_doubleWishartHD_sweep.R    # override probe count
#   RECOMPILE=1 Rscript test_doubleWishartHD_sweep.R  # force package rebuild
#
# Inside RStudio (no external Rscript needed):
#   source(system.file("validation", "test_doubleWishartHD_sweep.R", package = "rootWishartHD"))
#   res <- run_sweep(c("p40", "p150"))        # run selected settings
#   res <- run_sweep("p100", n_exact = 8)      # more probe points
#   res <- run_sweep("p150", plot = TRUE)      # also write CDF comparison PNGs
#   perf_table(res, "kable")                   # markdown performance table
#   cmp <- cdf_compare("p40"); plot_cdf_compare(cmp)   # MC vs analytic CDF
#   # regulate precision:
#   ex <- make_exact(start_digits10 = 150, stage1_max = 800, stage2_max = 20000)
#   run_sweep("p200", exact = ex, n_cores = 8)
# =============================================================================

suppressPackageStartupMessages({
  recompile <- identical(Sys.getenv("RECOMPILE"), "1")
  if (requireNamespace("rootWishartHD", quietly = TRUE) && !recompile) {
    library(rootWishartHD)
  } else if (requireNamespace("pkgload", quietly = TRUE)) {
    pkgload::load_all("rootWishartHD", quiet = TRUE, recompile = recompile,
                      compile = recompile)
  } else {
    stop("rootWishartHD is not installed and pkgload is unavailable.")
  }
  library(rWishart)
  library(corpcor)
})

options(rootWishartHD.force_multiprecision = TRUE)
Sys.setenv(OMP_NUM_THREADS = 1, OPENBLAS_NUM_THREADS = 1, MKL_NUM_THREADS = 1,
           BLIS_NUM_THREADS = 1, VECLIB_MAXIMUM_THREADS = 1)
# Pin BLAS/OMP threads at runtime too: env vars do not affect an already-loaded
# OpenBLAS, and oversubscribing 48 threads on the small (p x p) eigendecompositions
# in the MC step is much SLOWER than single-threaded. Workers set this themselves.
if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
  try(RhpcBLASctl::blas_set_num_threads(1L), silent = TRUE)
  try(RhpcBLASctl::omp_set_num_threads(1L), silent = TRUE)
}

# -----------------------------------------------------------------------------
# Configuration
# -----------------------------------------------------------------------------
CACHE_DIR <- file.path("diag_cache")
dir.create(CACHE_DIR, showWarnings = FALSE, recursive = TRUE)
SEED_SIM  <- 1L
FAST_SUB  <- 300L   # subsample size for the (diagnostic) fast double CDF pass

# Parallelization (same tools as test_doubleWishartHD_pit_fast_parallel.R):
#   BiocParallel::SnowParam (PSOCK/SOCK) if available, else base parallel PSOCK.
# PSOCK workers are fresh R processes, so rootWishartHD must be INSTALLED
# (R CMD INSTALL rootWishartHD), not merely load_all()'d.
N_CORES    <- {
  e <- Sys.getenv("N_CORES")
  if (nzchar(e)) as.integer(e) else min(16L, parallel::detectCores())
}
CHUNK_SIZE <- 1L    # probe points per worker task (1 = maximum parallelism)

# Precision knobs for the exact multiprecision logSF path, with a TWO-STAGE strategy:
#   stage 1 = moderate precision (start_digits10 -> stage1_max), used for ALL points.
#   stage 2 = high precision (-> stage2_max), used ONLY for points that stage 1
#             could not resolve (logSF returned -Inf / NaN because 1-F underflowed
#             at stage-1 precision). This keeps the common case fast while still
#             resolving deep upper-tail probes at large s.
# Tune these to trade speed vs. tail depth.
make_exact <- function(start_digits10 = 200L,
                       stage1_max     = 600L,
                       stage2_max     = 20000L,
                       tol1           = 1e-8,
                       tol2           = 1e-12,
                       pf_method      = "gauss",
                       direct_b       = FALSE,
                       scale_iter     = 12L) {
  list(start_digits10 = as.integer(start_digits10),
       stage1_max     = as.integer(stage1_max),
       stage2_max     = as.integer(stage2_max),
       tol1 = tol1, tol2 = tol2,
       pf_method = match.arg(pf_method, c("gauss", "lu", "auto", "svd", "schur")),
       direct_b = isTRUE(direct_b),
       scale_iter = as.integer(scale_iter))
}
# Default precision: moderate for most cases, but boosted for large-s settings (p200+)
# to avoid deep-tail underflows at the default 600-digit cap. The two-stage fallback
# will kick in on those points only if needed.
DEFAULT_EXACT <- make_exact(stage1_max = 800L)

# Parameter grid. `name` is the CLI selector. nrep is the MC sample size used
# for the cached draws; it is scaled DOWN as p grows because simulating large
# singular-Wishart matrices is itself expensive. n_exact = number of upper-tail
# probe points sent through the exact (slow) path; set 0 to skip exact entirely.
GRID <- list(
  list(name = "p40",  p = 40,  m_df = 28,  q_df = 20,  nrep = 3000, n_exact = 5),
  list(name = "p100", p = 100, m_df = 70,  q_df = 50,  nrep = 3000, n_exact = 5),
  list(name = "p150", p = 150, m_df = 100, q_df = 35,  nrep = 2000, n_exact = 5),
  list(name = "p200", p = 200, m_df = 140, q_df = 100, nrep = 1500, n_exact = 3),
  list(name = "p300", p = 300, m_df = 250, q_df = 150, nrep = 800,  n_exact = 2),
  list(name = "p500", p = 500, m_df = 499, q_df = 498, nrep = 300,  n_exact = 1)
)

# -----------------------------------------------------------------------------
# Helpers
# -----------------------------------------------------------------------------

# Identity-DSB reduction parameters (as used by the test scripts).
dsb_params <- function(p, m_df, q_df) {
  s   <- q_df
  mC  <- 0.5 * (abs(m_df - s) - 1)
  df2 <- p - m_df + q_df
  nC  <- 0.5 * (abs(df2 - s) - 1)
  list(s = s, m = mC, n = nC)
}

# One Monte Carlo draw of theta = lambda / (1 + lambda) for Roy's largest root.
one_draw_theta <- function(dfA, dfB, S) {
  A <- rWishart::rSingularWishart(1L, dfA, S)[, , 1]
  B <- rWishart::rSingularWishart(1L, dfB, S)[, , 1]

  eigA   <- eigen(A, symmetric = TRUE)
  rankVr <- min(dfA, nrow(A))
  V      <- eigA$vectors[, 1:rankVr, drop = FALSE]
  vals   <- eigA$values[1:rankVr]
  vmax   <- max(vals)
  vals   <- pmax(vals, 1e-12 * vmax)
  Xp     <- sweep(V, 2, sqrt(1 / vals), `*`)

  C    <- crossprod(Xp, B %*% Xp)
  sv   <- corpcor::fast.svd(C)
  W    <- Xp %*% sv$u
  lmax <- max(crossprod(W, B %*% W))
  lmax / (1 + lmax)
}

# Simulate (and cache) the MC theta vector for a setting. Returns list(theta, sim_sec, cached).
get_theta_cached <- function(p, m_df, q_df, nrep, seed = SEED_SIM) {
  key  <- sprintf("theta_p%d_m%d_q%d_n%d_s%d.rds", p, m_df, q_df, nrep, seed)
  path <- file.path(CACHE_DIR, key)
  if (file.exists(path)) {
    obj <- readRDS(path)
    obj$cached <- TRUE
    return(obj)
  }
  S <- diag(1, p)
  set.seed(seed)
  t0 <- proc.time()
  theta <- replicate(nrep, one_draw_theta(m_df, q_df, S))
  theta <- pmin(pmax(theta, 1e-12), 1 - 1e-12)
  sim_sec <- (proc.time() - t0)[["elapsed"]]
  obj <- list(theta = theta, sim_sec = sim_sec)
  saveRDS(obj, path)
  obj$cached <- FALSE
  obj
}

fmt_sec <- function(s) {
  if (is.na(s)) return("     n/a")
  if (s >= 90) sprintf("%6.2f min", s / 60) else sprintf("%7.3f s", s)
}

# Pick n upper-tail probe points spread across the empirical distribution.
pick_probe_points <- function(theta, n) {
  if (n <= 0) return(numeric(0))
  probs <- if (n == 1) 0.95 else seq(0.5, 0.999, length.out = n)
  pts <- as.numeric(quantile(theta, probs = probs, names = FALSE))
  pmin(pmax(pts, 1e-12), 1 - 1e-12)
}

exact_settings_label <- function(exact) {
  sprintf("adaptive %d->%d (stage2 ->%d), %s, direct_b=%s, scale_iter=%d",
          exact$start_digits10, exact$stage1_max, exact$stage2_max,
          exact$pf_method, exact$direct_b, exact$scale_iter)
}

# Look up a grid setting by name, or pass a config list directly.
get_cfg <- function(x) {
  if (is.list(x) && !is.null(x$p)) return(x)
  names_all <- vapply(GRID, `[[`, "", "name")
  i <- match(x, names_all)
  if (is.na(i)) stop("Unknown setting '", x, "'. Available: ",
                     paste(names_all, collapse = ", "))
  GRID[[i]]
}

# -----------------------------------------------------------------------------
# Parallel backend (PSOCK), mirroring test_doubleWishartHD_pit_fast_parallel.R
# -----------------------------------------------------------------------------

# Create a reusable PSOCK backend. Returns a handle used by parallel_logSF().
make_backend <- function(n_cores) {
  if (n_cores <= 1L) return(list(type = "seq"))
  if (requireNamespace("BiocParallel", quietly = TRUE)) {
    param <- BiocParallel::SnowParam(workers = n_cores, type = "SOCK",
                                     progressbar = FALSE)
    param <- BiocParallel::bpstart(param)   # keep workers alive across settings
    list(type = "bioc", param = param)
  } else {
    cl <- parallel::makeCluster(n_cores, type = "PSOCK")
    list(type = "base", cl = cl)
  }
}

stop_backend <- function(backend) {
  if (is.null(backend)) return(invisible())
  if (identical(backend$type, "base")) {
    try(parallel::stopCluster(backend$cl), silent = TRUE)
  } else if (identical(backend$type, "bioc")) {
    try(BiocParallel::bpstop(backend$param), silent = TRUE)
  }
}

# Build the per-setting worker function. Globals (params + precision settings)
# are captured into the closure so PSOCK serialization ships them to workers.
make_eval_chunk <- function(par, exact) {
  s0 <- par$s; m0 <- par$m; n0 <- par$n
  sd <- exact$start_digits10; s1 <- exact$stage1_max; s2 <- exact$stage2_max
  t1 <- exact$tol1; t2 <- exact$tol2
  pfm <- exact$pf_method; db <- exact$direct_b; si <- exact$scale_iter
  libs <- .libPaths()                          # ship master's library paths
  pkg_dir <- normalizePath("rootWishartHD")    # source dir for load_all fallback

  function(th_vec) {
    .libPaths(libs)
    if (!("rootWishartHD" %in% loadedNamespaces())) {
      suppressPackageStartupMessages({
        # Prefer the installed package (user's environment); fall back to loading
        # from source when it is not cleanly installed (this dev machine).
        if (requireNamespace("rootWishartHD", quietly = TRUE)) {
          library(rootWishartHD)
        } else {
          pkgload::load_all(pkg_dir, quiet = TRUE, compile = FALSE)
        }
      })
    }
    options(rootWishartHD.force_multiprecision = TRUE)
    Sys.setenv(OMP_NUM_THREADS = 1, OPENBLAS_NUM_THREADS = 1, MKL_NUM_THREADS = 1,
               BLIS_NUM_THREADS = 1, VECLIB_MAXIMUM_THREADS = 1)

    eval_at <- function(x, maxd, tol) {
      doubleWishart_log(
        x, s = s0, m = m0, n = n0,
        type = "arbitrary", tail = "upper", adaptive = TRUE,
        start_digits10 = sd, max_digits10 = maxd, tol = tol,
        pf_method = pfm, direct_b = db, scale_iter = si, verbose = FALSE)
    }

    th_vec <- as.numeric(th_vec)
    k <- length(th_vec)
    logSF <- rep(NA_real_, k); secs <- rep(NA_real_, k); stage <- rep(1L, k)
    for (i in seq_len(k)) {
      tt <- system.time({
        v <- eval_at(th_vec[i], s1, t1)          # stage 1 (moderate precision)
        if (!is.finite(v)) {                     # unresolved -> stage 2 (high precision)
          v <- eval_at(th_vec[i], s2, t2)
          stage[i] <- 2L
        }
      })[["elapsed"]]
      logSF[i] <- v; secs[i] <- tt
    }
    list(theta = th_vec, logSF = logSF, secs = secs, stage = stage)
  }
}

# Evaluate exact logSF on `pts` in parallel across PSOCK workers (chunked).
# Returns list(theta, logSF, secs) where secs is per-point COMPUTE time.
parallel_logSF <- function(pts, par, exact, backend, chunk_size = CHUNK_SIZE) {
  if (!length(pts)) return(list(theta = numeric(0), logSF = numeric(0), secs = numeric(0)))
  fn <- make_eval_chunk(par, exact)
  chunks <- split(pts, ceiling(seq_along(pts) / chunk_size))

  if (identical(backend$type, "seq")) {
    res <- lapply(chunks, fn)
  } else if (identical(backend$type, "bioc")) {
    res <- BiocParallel::bplapply(chunks, fn, BPPARAM = backend$param)
  } else {
    res <- parallel::parLapply(backend$cl, chunks, fn)
  }
  list(
    theta = unlist(lapply(res, `[[`, "theta"), use.names = FALSE),
    logSF = unlist(lapply(res, `[[`, "logSF"), use.names = FALSE),
    secs  = unlist(lapply(res, `[[`, "secs"),  use.names = FALSE),
    stage = unlist(lapply(res, `[[`, "stage"), use.names = FALSE)
  )
}

# -----------------------------------------------------------------------------
# Per-setting runner
# -----------------------------------------------------------------------------
run_one <- function(cfg, backend, exact = DEFAULT_EXACT, n_exact_override = NA,
                    n_cores_label = NA) {
  par <- dsb_params(cfg$p, cfg$m_df, cfg$q_df)
  n_exact <- if (!is.na(n_exact_override)) n_exact_override else cfg$n_exact

  cat(sprintf("\n=== %-5s p=%d m_df=%d q_df=%d | s=%d m=%.2f n=%.2f | nrep=%d n_exact=%d ===\n",
              cfg$name, cfg$p, cfg$m_df, cfg$q_df, par$s, par$m, par$n, cfg$nrep, n_exact))

  # (A) Monte Carlo (cached)
  sim <- get_theta_cached(cfg$p, cfg$m_df, cfg$q_df, cfg$nrep)
  theta <- sim$theta
  cat(sprintf("  [A] MC sim   : %s  (%s, nrep=%d)\n",
              fmt_sec(sim$sim_sec), if (sim$cached) "cached" else "fresh", length(theta)))

  # (B) fast double CDF pass on a subsample (single Pfaffian backend = gauss).
  # This only contextualizes how often the double path saturates; it is not the
  # quantity of interest. The Pfaffian is O(s^3), so we scale the subsample DOWN
  # with s to keep this diagnostic cheap (~a few seconds) even at s ~ 500.
  n_fast <- max(5L, min(FAST_SUB, length(theta), as.integer(round(3e8 / par$s^3))))
  sub <- theta[seq_len(n_fast)]
  tB <- system.time(
    Ffast <- doubleWishart_scaled(sub, s = par$s, m = par$m, n = par$n,
                                         type = "double", force_multiprecision = FALSE,
                                         pf_method = "gauss",
                                         adaptive = FALSE, verbose = FALSE))[["elapsed"]]
  sat <- mean(Ffast <= 1e-12 | Ffast >= 1 - 1e-12)
  cat(sprintf("  [B] fast CDF : %s  (n=%d subsample, saturated fraction = %.1f%%)\n",
              fmt_sec(tB), length(sub), 100 * sat))

  # (C) exact multiprecision logSF on probe points -- PARALLEL across PSOCK workers.
  pts <- pick_probe_points(theta, n_exact)
  per_pt <- numeric(0)
  wallC <- NA_real_
  n_stage2 <- 0L
  n_inf <- 0L
  if (length(pts)) {
    cat(sprintf("  [C] exact logSF (upper tail; %s) on %d pts, %d workers:\n",
                exact_settings_label(exact), length(pts),
                if (identical(backend$type, "seq")) 1L else n_cores_label))
    wallC <- system.time(out <- parallel_logSF(pts, par, exact, backend))[["elapsed"]]
    per_pt <- out$secs
    n_stage2 <- sum(out$stage == 2L, na.rm = TRUE)
    n_inf <- sum(!is.finite(out$logSF))
    ord <- order(out$theta)
    for (i in ord) {
      cat(sprintf("        theta=%.6f  logSF=%-14.6f  stage=%d  compute=%s\n",
                  out$theta[i], out$logSF[i], out$stage[i], fmt_sec(out$secs[i])))
    }
    cat(sprintf("        wall=%s  (sum compute=%s, speedup=%.1fx; stage2=%d, non-finite=%d)\n",
                fmt_sec(wallC), fmt_sec(sum(per_pt)),
                if (wallC > 0) sum(per_pt) / wallC else 1, n_stage2, n_inf))
  }

  data.frame(
    name = cfg$name, p = cfg$p, s = par$s, nrep = length(theta),
    mc_sec = round(sim$sim_sec, 2), fast_sec = round(tB, 3),
    exact_pts = length(pts),
    exact_per_pt_sec = round(if (length(per_pt)) mean(per_pt, na.rm = TRUE) else NA, 3),
    exact_wall_sec = round(wallC, 2),
    stage2 = n_stage2, non_finite = n_inf,
    stringsAsFactors = FALSE
  )
}

# -----------------------------------------------------------------------------
# CDF comparison: empirical (Monte Carlo) vs analytic (rootWishartHD)
# -----------------------------------------------------------------------------

# Compute the empirical CDF (from cached MC draws) and the analytic CDF from the
# package on a grid of theta values. Returns a data.frame(theta, F_emp, F_ana, abs_err).
# The analytic CDF uses the robust multiprecision log path: F = exp(logCDF). For large s this
# is expensive, so keep n_grid modest (analytic cost ~ n_grid * single-eval cost).
cdf_compare <- function(setting, n_grid = 25L, exact = DEFAULT_EXACT,
                        probs = NULL) {
  cfg <- get_cfg(setting)
  par <- dsb_params(cfg$p, cfg$m_df, cfg$q_df)
  sim <- get_theta_cached(cfg$p, cfg$m_df, cfg$q_df, cfg$nrep)
  theta <- sim$theta

  if (is.null(probs)) probs <- seq(0.02, 0.98, length.out = n_grid)
  grid <- as.numeric(quantile(theta, probs = probs, names = FALSE))
  grid <- sort(unique(pmin(pmax(grid, 1e-12), 1 - 1e-12)))

  logF <- doubleWishart_log(
    grid, s = par$s, m = par$m, n = par$n,
    type = "arbitrary", tail = "lower", adaptive = TRUE,
    start_digits10 = exact$start_digits10, max_digits10 = exact$stage1_max,
    tol = exact$tol1, pf_method = exact$pf_method, direct_b = exact$direct_b,
    scale_iter = exact$scale_iter, verbose = FALSE)
  F_ana <- exp(pmin(0, logF))
  F_emp <- stats::ecdf(theta)(grid)

  out <- data.frame(theta = grid, F_emp = F_emp, F_ana = F_ana,
                    abs_err = abs(F_emp - F_ana))
  attr(out, "setting") <- cfg$name
  attr(out, "s") <- par$s
  out
}

# Plot empirical vs analytic CDF. If `file` is given, writes a PNG.
plot_cdf_compare <- function(cmp, title = NULL, file = NULL) {
  if (is.null(title)) {
    title <- sprintf("%s (s=%s): empirical (MC) vs analytic CDF",
                     attr(cmp, "setting"), attr(cmp, "s"))
  }
  if (!is.null(file)) {
    dir.create(dirname(file), showWarnings = FALSE, recursive = TRUE)
    grDevices::png(file, width = 1000, height = 700, res = 130)
    on.exit(grDevices::dev.off(), add = TRUE)
  }
  op <- graphics::par(no.readonly = TRUE); on.exit(graphics::par(op), add = TRUE)
  plot(cmp$theta, cmp$F_emp, type = "s", col = "steelblue", lwd = 2,
       xlab = expression(theta), ylab = expression(F(theta)),
       main = title, ylim = c(0, 1))
  graphics::lines(cmp$theta, cmp$F_ana, col = "firebrick", lwd = 2, lty = 2)
  graphics::points(cmp$theta, cmp$F_ana, col = "firebrick", pch = 19, cex = 0.5)
  graphics::legend("bottomright", bty = "n", lwd = 2, lty = c(1, 2),
                   col = c("steelblue", "firebrick"),
                   legend = c("empirical (Monte Carlo)", "analytic (rootWishartHD)"))
  invisible(cmp)
}

# -----------------------------------------------------------------------------
# Performance table
# -----------------------------------------------------------------------------

# Format a run_sweep() summary as a performance table. format="kable" returns a
# markdown table (needs knitr); otherwise a tidy data.frame with friendly names.
perf_table <- function(summary, format = c("df", "kable")) {
  format <- match.arg(format)
  tab <- data.frame(
    Setting        = summary$name,
    p              = summary$p,
    s              = summary$s,
    `MC reps`      = summary$nrep,
    `MC sim (s)`   = summary$mc_sec,
    `fast CDF (s)` = summary$fast_sec,
    `exact pts`    = summary$exact_pts,
    `exact/pt (s)` = summary$exact_per_pt_sec,
    `exact wall (s)` = summary$exact_wall_sec,
    `stage2`       = summary$stage2,
    check.names = FALSE, stringsAsFactors = FALSE
  )
  if (format == "kable") {
    if (!requireNamespace("knitr", quietly = TRUE))
      stop("format='kable' requires the 'knitr' package.")
    return(knitr::kable(tab, format = "pipe", row.names = FALSE,
                        caption = "rootWishartHD performance across parameter settings"))
  }
  tab
}

# -----------------------------------------------------------------------------
# Main sweep
# -----------------------------------------------------------------------------

# which       : character vector of setting names (NULL = whole GRID), or config lists.
# n_exact     : override number of exact probe points per setting (NA = grid default).
# exact       : precision settings from make_exact().
# n_cores     : PSOCK workers (1 = sequential).
# plot        : if TRUE, also write empirical-vs-analytic CDF PNGs to plot_dir.
# plot_dir    : directory for CDF figures.
# cdf_grid    : grid size for the CDF comparison plots.
# cdf_max_s   : skip CDF plots when s exceeds this (analytic CDF gets too slow).
run_sweep <- function(which = NULL, n_exact = NA, exact = DEFAULT_EXACT,
                      n_cores = N_CORES, plot = FALSE, plot_dir = "diag_figs",
                      cdf_grid = 25L, cdf_max_s = 120L) {
  names_all <- vapply(GRID, `[[`, "", "name")
  if (is.null(which) || !length(which)) which <- names_all
  sel <- lapply(which, get_cfg)

  backend <- make_backend(n_cores)
  on.exit(stop_backend(backend), add = TRUE)
  cat(sprintf("Parallel backend: %s, workers=%d | precision: %s\n",
              backend$type, if (identical(backend$type, "seq")) 1L else n_cores,
              exact_settings_label(exact)))

  t0 <- proc.time()
  rows <- lapply(sel, run_one, backend = backend, exact = exact,
                 n_exact_override = n_exact, n_cores_label = n_cores)
  summary <- do.call(rbind, rows)

  cat("\n================ SWEEP SUMMARY ================\n")
  print(summary, row.names = FALSE)
  cat(sprintf("\nTotal wall time: %s\n", fmt_sec((proc.time() - t0)[["elapsed"]])))

  if (isTRUE(plot)) {
    dir.create(plot_dir, showWarnings = FALSE, recursive = TRUE)
    for (cfg in sel) {
      s <- dsb_params(cfg$p, cfg$m_df, cfg$q_df)$s
      if (s > cdf_max_s) {
        cat(sprintf("  [plot] skip %s (s=%d > cdf_max_s=%d; analytic CDF too slow)\n",
                    cfg$name, s, cdf_max_s)); next
      }
      cmp <- tryCatch(cdf_compare(cfg, n_grid = cdf_grid, exact = exact),
                      error = function(e) { cat("  [plot] error:", conditionMessage(e), "\n"); NULL })
      if (!is.null(cmp)) {
        f <- file.path(plot_dir, paste0("cdf_", cfg$name, ".png"))
        plot_cdf_compare(cmp, file = f)
        cat(sprintf("  [plot] %s  (max |F_emp - F_ana| = %.4f)\n", f, max(cmp$abs_err)))
      }
    }
  }
  invisible(summary)
}

if (sys.nframe() == 0L) {
  args <- commandArgs(trailingOnly = TRUE)
  nx_env <- Sys.getenv("N_EXACT")
  nx <- if (nzchar(nx_env)) as.integer(nx_env) else NA
  do_plot <- identical(Sys.getenv("PLOT"), "1")
  run_sweep(which = if (length(args)) args else NULL, n_exact = nx, plot = do_plot)
}
