## =====================================================================
## benchmark_family.R
##
## A small, smart benchmark harness for the rootWishartHD *function
## family*: the CDF backends, the log-scale tail functions, and the
## applied p-value / critical-value wrappers, for both the single- and
## double-Wishart settings.
##
## It is a companion to test_doubleWishartHD_sweep.R (which focuses on
## scaling doubleWishart_log() to large p). This harness
## instead answers "which function should I call, and what does it
## cost?" across the whole exported API.
##
## Design goals
##   * sourceable from RStudio (no auto-run on source; see bottom guard)
##   * every call is timed and tagged with the method actually used
##     (fast double vs exact multiprecision) so the usability trade-off is visible
##   * exact (slow) evaluations are cached on disk and reused
##   * small default cases so the whole suite runs in seconds-to-minutes
##
## Quick start (RStudio):
##   source(system.file("validation", "benchmark_family.R", package = "rootWishartHD"))
##   res <- family_benchmark()      # curated small grid
##   family_report(res)             # compact usability table
##   bench_saturation(dsb(80,54,40))# show double-CDF saturation vs logSF
## =====================================================================

suppressMessages({
  if (!requireNamespace("rootWishartHD", quietly = TRUE)) {
    if (requireNamespace("pkgload", quietly = TRUE) &&
        dir.exists("rootWishartHD")) {
      pkgload::load_all("rootWishartHD", quiet = TRUE)
    } else {
      stop("Install rootWishartHD or run from a dir containing rootWishartHD/.")
    }
  } else {
    library(rootWishartHD)
  }
})
options(rootWishartHD.force_multiprecision = FALSE)

## ---------------------------------------------------------------------
## Parameter mapping and small helpers
## ---------------------------------------------------------------------

## MANOVA-style (p, m_df, q_df) -> Chiani (s, m, n) for the double Wishart.
dsb <- 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(p = p, m_df = m_df, q_df = q_df, s = s, m = mC, n = nC)
}

## time an expression; return value with elapsed seconds attached
.timed <- function(expr) {
  t <- system.time(v <- eval.parent(substitute(expr)))
  attr(v, "sec") <- unname(t[["elapsed"]])
  v
}

.cache_dir <- "bench_cache"
.cache_get <- function(key, compute) {
  if (!dir.exists(.cache_dir)) dir.create(.cache_dir, showWarnings = FALSE)
  f <- file.path(.cache_dir, paste0(key, ".rds"))
  if (file.exists(f)) {
    v <- readRDS(f); attr(v, "cached") <- TRUE; return(v)
  }
  v <- compute()
  saveRDS(v, f); attr(v, "cached") <- FALSE; v
}

## ---------------------------------------------------------------------
## 1. Double-Wishart CDF backends: agreement + speed
## ---------------------------------------------------------------------
## Compares the two double-precision CDF entry points
##   doubleWishart(), doubleWishart_scaled()
## on the same points. They should agree to ~1e-6 away from the deep
## lower tail; the scaled backend adds scaling/equilibration that
## matters only when the matrix is badly conditioned (large s).
bench_double_backends <- function(cfg, theta = seq(0.5, 0.97, length.out = 8)) {
  f <- function(fun) .timed(fun(theta, s = cfg$s, m = cfg$m, n = cfg$n,
                                type = "double", verbose = FALSE))
  v0 <- f(doubleWishart)
  v1 <- f(doubleWishart_scaled)
  data.frame(
    setting   = sprintf("p%d (s=%d)", cfg$p, cfg$s),
    backend   = c("doubleWishart", "scaled"),
    sec       = c(attr(v0, "sec"), attr(v1, "sec")),
    max_abs_vs_base = c(0, max(abs(v1 - v0))),
    stringsAsFactors = FALSE
  )
}

## ---------------------------------------------------------------------
## 2. Saturation: where the double CDF dies and the log tail survives
## ---------------------------------------------------------------------
## The plain double CDF saturates to exactly 1 in the upper tail (1 - F
## underflows). doubleWishart_log(tail="upper") returns a
## finite logSF instead. This is the single most important reason the
## log-scale functions exist.
bench_saturation <- function(cfg, theta = c(0.90, 0.95, 0.99)) {
  Fd  <- doubleWishart(theta, s = cfg$s, m = cfg$m, n = cfg$n,
                       type = "double", verbose = FALSE)
  key <- sprintf("sat_p%d_s%d_%s", cfg$p, cfg$s,
                 paste(formatC(theta, format = "g"), collapse = "_"))
  lse <- .cache_get(key, function()
    doubleWishart_log(theta, s = cfg$s, m = cfg$m, n = cfg$n,
                                   type = "arbitrary", tail = "upper",
                                   verbose = FALSE))
  data.frame(
    setting       = sprintf("p%d (s=%d)", cfg$p, cfg$s),
    theta         = theta,
    CDF_double    = Fd,
    one_minus_CDF = 1 - Fd,        # = 0/eps once saturated
    logSF_exact   = as.numeric(lse),
    SF_exact      = exp(as.numeric(lse)),
    saturated     = (1 - Fd) <= .Machine$double.eps,
    stringsAsFactors = FALSE
  )
}

## ---------------------------------------------------------------------
## 3. Applied wrappers: p-value + critical value, with method + timing
## ---------------------------------------------------------------------
## doubleWishart_pvalue() and doubleWishart_qalpha() pick a fast double
## pass when the tail is shallow and only fall back to the (slow) exact
## multiprecision path when the double CDF has saturated. We report which path was
## taken so the cost is transparent, and round-trip qalpha through pvalue.
bench_double_applied <- function(cfg, alphas = c(0.05, 0.01),
                                 probe = c(0.90, 0.999)) {
  pv <- .timed(doubleWishart_pvalue(probe, s = cfg$s, m = cfg$m, n = cfg$n,
                                    input = "theta", verbose = FALSE))
  rows <- vector("list", length(alphas))
  for (i in seq_along(alphas)) {
    a  <- alphas[i]
    qa <- .timed(doubleWishart_qalpha(a, s = cfg$s, m = cfg$m, n = cfg$n,
                                      return = "theta", verbose = FALSE))
    if (is.finite(as.numeric(qa))) {
      bk <- doubleWishart_pvalue(as.numeric(qa), s = cfg$s, m = cfg$m, n = cfg$n,
                                 input = "theta", verbose = FALSE)
      p_back <- as.numeric(bk)
      rel_err <- abs(p_back - a) / a
    } else {
      p_back <- NA_real_
      rel_err <- NA_real_
    }
    rows[[i]] <- data.frame(
      setting    = sprintf("p%d (s=%d)", cfg$p, cfg$s),
      alpha      = a,
      theta_crit = as.numeric(qa),
      p_back     = p_back,
      rel_err    = rel_err,
      qalpha_sec = attr(qa, "sec"),
      qalpha_method = paste(unique(attr(qa, "method")), collapse = ","),
      stringsAsFactors = FALSE
    )
  }
  list(
    pvalue = data.frame(
      setting = sprintf("p%d (s=%d)", cfg$p, cfg$s),
      theta   = probe,
      p       = as.numeric(pv),
      method  = attr(pv, "method"),
      sec     = attr(pv, "sec"),
      stringsAsFactors = FALSE),
    qalpha = do.call(rbind, rows)
  )
}

## ---------------------------------------------------------------------
## 4. Single-Wishart family: CDF / logSF / pvalue / qalpha
## ---------------------------------------------------------------------
## x is on the eigenvalue scale; CDF = P(lambda_max <= x). The exported
## single-Wishart functions take (x, n_min, n_max) directly (no z2u).
bench_single <- function(n_min, n_max, x = c(10, 20, 30, 40),
                         alphas = c(0.05, 0.01)) {
  F2  <- .timed(singleWishart_cdf(x, n_min = n_min, n_max = n_max,
                                            type = "double", verbose = FALSE))
  key <- sprintf("sW_log_%d_%d_%s", n_min, n_max,
                 paste(formatC(x, format = "g"), collapse = "_"))
  lse <- .cache_get(key, function()
    singleWishart_log(x, n_min = n_min, n_max = n_max,
                                   type = "arbitrary", tail = "upper",
                                   verbose = FALSE))
  pv  <- singleWishart_pvalue(x, n_min = n_min, n_max = n_max, verbose = FALSE)
  qa  <- vapply(alphas, function(a)
    as.numeric(singleWishart_qalpha(a, n_min = n_min, n_max = n_max,
                                    verbose = FALSE)), numeric(1))
  qa_back <- vapply(qa, function(z)
    as.numeric(singleWishart_pvalue(z, n_min = n_min, n_max = n_max,
                                    verbose = FALSE)), numeric(1))
  list(
    cdf = data.frame(
      pair = sprintf("(%d,%d)", n_min, n_max),
      x = x,
      CDF_double  = as.numeric(F2),
      CDF_from_SF = -expm1(as.numeric(lse)),
      logSF_exact = as.numeric(lse),
      p_value     = as.numeric(pv),
      method      = attr(pv, "method"),
      stringsAsFactors = FALSE),
    qalpha = data.frame(
      pair = sprintf("(%d,%d)", n_min, n_max),
      alpha = alphas, x_crit = qa, p_back = qa_back,
      rel_err = abs(qa_back - alphas) / alphas,
      stringsAsFactors = FALSE),
    cdf_sec = attr(F2, "sec")
  )
}

## ---------------------------------------------------------------------
## Curated grid + top-level driver
## ---------------------------------------------------------------------
## Small by default so the suite finishes quickly. Larger settings are
## included but exact evaluations are cached, so the second run is fast.
FAMILY_GRID <- list(
  small = dsb(20, 14, 10),   # s = 10
  mid   = dsb(40, 26, 20),   # s = 20
  large = dsb(80, 54, 40)    # s = 40 (exact path is the slow one)
)

## Run the whole family benchmark and return a list of tidy tables.
family_benchmark <- function(which = c("small", "mid"),
                             sat_theta = c(0.90, 0.95, 0.99),
                             single = list(c(5, 10), c(8, 16))) {
  which <- intersect(which, names(FAMILY_GRID))
  backends <- do.call(rbind, lapply(FAMILY_GRID[which], bench_double_backends))
  satur    <- do.call(rbind, lapply(FAMILY_GRID[which],
                                    function(cfg) bench_saturation(cfg, sat_theta)))
  applied  <- lapply(FAMILY_GRID[which], function(cfg) {
    tryCatch(bench_double_applied(cfg), error = function(e) {
      list(error = conditionMessage(e),
           pvalue = data.frame(setting = sprintf("p%d (s=%d)", cfg$p, cfg$s),
                               theta = NA_real_, p = NA_real_, method = "error",
                               sec = NA_real_, stringsAsFactors = FALSE),
           qalpha = data.frame(setting = sprintf("p%d (s=%d)", cfg$p, cfg$s),
                               alpha = NA_real_, theta_crit = NA_real_,
                               p_back = NA_real_, rel_err = NA_real_,
                               qalpha_sec = NA_real_, qalpha_method = "error",
                               stringsAsFactors = FALSE))
    })
  })
  singles  <- lapply(single, function(nn) bench_single(nn[1], nn[2]))
  rownames(backends) <- NULL; rownames(satur) <- NULL
  list(backends = backends, saturation = satur,
       applied = applied, single = singles)
}

## Compact usability summary printed to the console.
family_report <- function(res) {
  cat("================ rootWishartHD function-family benchmark ================\n\n")
  cat("[1] Double-Wishart CDF backends (double precision): speed + agreement\n")
  if (is.null(res) || is.null(res$backends)) {
    cat("  <no benchmark result available>
")
    return(invisible(res))
  }
  print(within(res$backends, {
    sec <- round(sec, 4); max_abs_vs_base <- signif(max_abs_vs_base, 3)
  }), row.names = FALSE)

  cat("\n[2] Saturation: plain double CDF vs exact log-survival (upper tail)\n")
  st <- res$saturation
  st$CDF_double <- round(st$CDF_double, 6)
  st$logSF_exact <- round(st$logSF_exact, 4)
  st$SF_exact <- signif(st$SF_exact, 4)
  print(st[, c("setting", "theta", "CDF_double", "saturated",
               "logSF_exact", "SF_exact")], row.names = FALSE)

  cat("\n[3] Applied wrappers: p-value method + qalpha round-trip\n")
  for (nm in names(res$applied)) {
    a <- res$applied[[nm]]
    cat(sprintf("  -- %s --\n", a$pvalue$setting[1]))
    pv <- a$pvalue; pv$p <- signif(pv$p, 4); pv$sec <- round(pv$sec, 3)
    print(pv[, c("theta", "p", "method", "sec")], row.names = FALSE)
    qa <- a$qalpha
    qa$theta_crit <- round(qa$theta_crit, 6)
    qa$p_back <- signif(qa$p_back, 4); qa$rel_err <- signif(qa$rel_err, 2)
    qa$qalpha_sec <- round(qa$qalpha_sec, 3)
    print(qa[, c("alpha", "theta_crit", "p_back", "rel_err",
                 "qalpha_sec", "qalpha_method")], row.names = FALSE)
  }

  cat("\n[4] Single-Wishart family (CDF / logSF / pvalue / qalpha)\n")
  for (s in res$single) {
    cat(sprintf("  -- pair %s --\n", s$cdf$pair[1]))
    cdf <- s$cdf
    cdf$CDF_double <- signif(cdf$CDF_double, 4)
    cdf$CDF_from_SF <- signif(cdf$CDF_from_SF, 4)
    cdf$logSF_exact <- round(cdf$logSF_exact, 4)
    cdf$p_value <- signif(cdf$p_value, 4)
    print(cdf[, c("x", "CDF_double", "CDF_from_SF", "logSF_exact",
                  "p_value", "method")], row.names = FALSE)
    qa <- s$qalpha; qa$x_crit <- round(qa$x_crit, 4)
    qa$p_back <- signif(qa$p_back, 4); qa$rel_err <- signif(qa$rel_err, 2)
    print(qa[, c("alpha", "x_crit", "p_back", "rel_err")], row.names = FALSE)
  }
  cat("\n========================================================================\n")
  invisible(res)
}

## Run only when executed via Rscript, not when sourced into RStudio.
if (sys.nframe() == 0L) {
  args <- commandArgs(trailingOnly = TRUE)
  which <- if (length(args)) args else c("small", "mid")
  res <- family_benchmark(which = which)
  family_report(res)
}
