## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.6,
  fig.align = "center"
)
has_mc <- requireNamespace("rWishart", quietly = TRUE) &&
          requireNamespace("corpcor", quietly = TRUE)

## ----load---------------------------------------------------------------------
library(rootWishartHD)

## ----dsb----------------------------------------------------------------------
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)
}

## ----cdf-basic----------------------------------------------------------------
par0 <- dsb_params(p = 20, m_df = 14, q_df = 10)   # s = 10
theta <- c(0.3, 0.5, 0.7, 0.9)
doubleWishart(theta, s = par0$s, m = par0$m, n = par0$n,
              type = "double", verbose = FALSE)

## ----multiprecision-option, eval=FALSE----------------------------------------
# options(rootWishartHD.force_multiprecision = TRUE)

## ----mpfr-enabled-------------------------------------------------------------
rootWishartHD_mpfr_enabled()

## ----logsf--------------------------------------------------------------------
# Upper-tail log-survival log(1 - F). 'lower' gives log F.
logSF <- doubleWishart_log(
  c(0.85, 0.92, 0.97), s = par0$s, m = par0$m, n = par0$n,
  type = "arbitrary", tail = "upper", verbose = FALSE)
data.frame(theta = c(0.85, 0.92, 0.97),
           logSF = logSF,
           SF    = exp(logSF))

## ----two-stage----------------------------------------------------------------
logsf_two_stage <- function(theta, s, m, n,
                            stage1_max = 600L, stage2_max = 20000L) {
  v <- doubleWishart_log(
    theta, s = s, m = m, n = n, type = "arbitrary", tail = "upper",
    adaptive = TRUE, start_digits10 = 200L, max_digits10 = stage1_max,
    tol = 1e-8, pf_method = "gauss", verbose = FALSE)
  unresolved <- !is.finite(v)
  if (any(unresolved)) {
    v[unresolved] <- doubleWishart_log(
      theta[unresolved], s = s, m = m, n = n, type = "arbitrary", tail = "upper",
      adaptive = TRUE, start_digits10 = 200L, max_digits10 = stage2_max,
      tol = 1e-12, pf_method = "gauss", verbose = FALSE)
  }
  v
}

logsf_two_stage(c(0.9, 0.97), s = par0$s, m = par0$m, n = par0$n)

## ----mc-helper, eval = has_mc-------------------------------------------------
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   <- pmax(eigA$values[1:rankVr], 1e-12 * max(eigA$values[1:rankVr]))
  Xp     <- sweep(V, 2, sqrt(1 / vals), `*`)
  C      <- crossprod(Xp, B %*% Xp)
  W      <- Xp %*% corpcor::fast.svd(C)$u
  lmax   <- max(crossprod(W, B %*% W))
  lmax / (1 + lmax)
}

## ----mc-figure, eval = has_mc, fig.cap = "Empirical (Monte Carlo) vs analytic CDF for a small double-Wishart case."----
set.seed(1)
p <- 20; m_df <- 14; q_df <- 10
par1 <- dsb_params(p, m_df, q_df)            # s = 10
S <- diag(1, p)
theta_mc <- replicate(100, one_draw_theta(m_df, q_df, S))
theta_mc <- pmin(pmax(theta_mc, 1e-12), 1 - 1e-12)

grid  <- as.numeric(quantile(theta_mc, probs = seq(0.02, 0.98, length.out = 30)))
F_ana <- doubleWishart(grid, s = par1$s, m = par1$m, n = par1$n,
                       type = "double", verbose = FALSE)
F_emp <- ecdf(theta_mc)(grid)

plot(grid, F_emp, type = "s", col = "steelblue", lwd = 2, ylim = c(0, 1),
     xlab = expression(theta), ylab = expression(F(theta)),
     main = sprintf("s = %d: empirical vs analytic CDF", par1$s))
lines(grid, F_ana, col = "firebrick", lwd = 2, lty = 2)
points(grid, F_ana, col = "firebrick", pch = 19, cex = 0.5)
legend("bottomright", bty = "n", lwd = 2, lty = c(1, 2),
       col = c("steelblue", "firebrick"),
       legend = c("empirical (Monte Carlo)", "analytic (rootWishartHD)"))

cat(sprintf("max |F_emp - F_ana| = %.4f\n", max(abs(F_emp - F_ana))))

## ----mc-skip, eval = !has_mc, echo = FALSE, results = "asis"------------------
# cat("> _The Monte Carlo validation figure is skipped because `rWishart` and/or",
#     "`corpcor` are not installed._\n")

## ----perf-table, echo = FALSE-------------------------------------------------
perf <- data.frame(
  setting   = c("p40", "p100", "p150", "p300", "p500"),
  p         = c(40, 100, 150, 300, 500),
  s         = c(20, 50, 35, 150, 498),
  `MC sim`  = c("4.1 s", "not run", "2.5 min", "2.5 min", "3.1 min"),
  `exact/pt`= c("2.5 s", "5.8 s", "4.1 s", "22 s", "3.7 min"),
  check.names = FALSE, stringsAsFactors = FALSE
)
knitr::kable(perf, caption = "Representative timings for one exact upper-tail logSF evaluation.")

## ----validation-sweep, eval=FALSE---------------------------------------------
# sweep_file <- system.file(
#   "validation", "test_doubleWishartHD_sweep.R",
#   package = "rootWishartHD"
# )
# 
# if (!nzchar(sweep_file)) {
#   stop("Validation script not found. Reinstall rootWishartHD with inst/validation included.")
# }
# 
# source(sweep_file)
# 
# # performance sweep over selected settings, 8 workers
# res <- run_sweep(c("p40", "p150"), n_exact = 6, n_cores = 8)
# perf_table(res, "kable")               # markdown performance table
# 
# # empirical-vs-analytic CDF comparison + figure
# cmp <- cdf_compare("p40", n_grid = 25)
# plot_cdf_compare(cmp)
# 
# # regulate precision: two-stage arbitrary-precision evaluation
# ex <- make_exact(start_digits10 = 150, stage1_max = 800, stage2_max = 20000)
# run_sweep("p200", exact = ex, plot = TRUE)   # writes CDF PNGs to diag_figs/