#partial code for illustration of the SVEM loop
coef_matrix <- matrix(NA, nrow = nBoot, ncol = p + 1)
 for (i in 1:nBoot) {
      U <- runif(n)
      w_train <- -log(U)
      w_valid <- -log(1 - U)
      #match glmnet normalization of training weight vector
      w_train <- w_train * (n / sum(w_train))
      w_valid <- w_valid * (n / sum(w_valid))
      glmnet(
          X, y_numeric,
          alpha = alpha,
          weights = w_train,
          intercept = TRUE,
          standardize = standardize,
          maxit = 1e6,
          nlambda = 500
      )
      predict(fit, newx = X)
      val_errors <- colSums(w_valid * (y_numeric - pred_valid)^2)
      k_values <- fit$df
      n_obs <- length(y_numeric)
      aic_values <- n_obs * log(val_errors / n_obs) + 2 * k_values
         # Choose lambda
      if (objective == "wSSE") {
        idx_min <- which.min(val_errors)
        lambda_opt <- fit$lambda[idx_min]
        val_error <- val_errors[idx_min]
      } else if (objective == "wAIC") {
        idx_min <- which.min(aic_values)
        lambda_opt <- fit$lambda[idx_min]
        val_error <- aic_values[idx_min]
      }
      coef_matrix[i, ] <- as.vector(coef(fit, s = lambda_opt))
}

library(SVEMnet)

# Example data
data <- iris
svem_model <- SVEMnet(Sepal.Length ~ ., data = data, nBoot = 300)
coef(svem_model)

plot(svem_model)

predictions <- predict(svem_model, data)
print(predictions)

# Simulate data
set.seed(1)
n <- 25
X1 <- runif(n)
X2 <- runif(n)
X3 <- runif(n)
X4 <- runif(n)
X5 <- runif(n)

#y only depends on X1 and X2
y <- 1 + X1 +  X2 + X1 * X2 + X1^2 + rnorm(n)
data <- data.frame(y, X1, X2, X3, X4, X5)

# Perform the SVEM significance test
test_result <- svem_significance_test_parallel(
  y ~ (X1 + X2 + X3)^2 + I(X1^2) + I(X2^2) + I(X3^2),
  data = data

)

# View the p-value
print(test_result)
SVEM Significance Test p-value:
[1] 0.009399093


test_result2 <- svem_significance_test_parallel(
  y ~ (X1 + X2 )^2 + I(X1^2) + I(X2^2),
  data = data
)

# View the p-value
print(test_result2)
SVEM Significance Test p-value:
[1] 0.006475736

#note that the response does not depend on X4 or X5
test_result3 <- svem_significance_test_parallel(
  y ~ (X4 + X5)^2 + I(X4^2) + I(X5^2),
  data = data
)

# View the p-value
print(test_result3)
SVEM Significance Test p-value:
[1] 0.8968502

# Plot the Mahalanobis distances
plot(test_result,test_result2,test_result3)

## ======================================================================
## SVEMnet boundary sims: wAIC vs wSSE + debias(TRUE/FALSE)
## - standardize = TRUE
## - glmnet_alpha = c(1)
## - No leakage; test is disjoint from train; debias refit check
## - n_train ranges chosen to cross n_train_used vs p_full boundary
## - Outputs:
##   * res ........ per-fit records (you can model this)
##   * A tidy summary with 4 win-rates:
##       winrate_debias_TRUE, winrate_debias_FALSE,
##       winrate_wAIC, winrate_wSSE
##   * Two simple plots (comment out prints if running headless)
## ======================================================================

if (!requireNamespace("SVEMnet", quietly = TRUE)) install.packages("SVEMnet")
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
library(SVEMnet)
library(ggplot2)

## ---------- Logging / Progress ----------
VERBOSE     <- TRUE     # per-fit messages
ECHO_EVERY  <- 1        # print every k-th fit
PROGRESS    <- TRUE     # progress bar

stamp   <- function() format(Sys.time(), "%H:%M:%S")
log_msg <- function(...) if (VERBOSE) {
  cat(sprintf("[%s] %s\n", stamp(), sprintf(...))); flush.console()
}
options(warn = 1)       # surface warnings immediately

set.seed(20250814)

## ------------------ Controls ------------------
REPS      <- 10           # increase to 20+ for more power
NBOOT     <- 200
R2_GRID   <- c(0.30, 0.70, 0.90)
RHO_GRID  <- c(0.2, 0.6, 0.9)
P_GRID    <- c(4, 8, 12)
NTEST     <- 1000
MODELS    <- c( "main_plus_int", "full_quadratic")

## ------------------ Helpers -------------------
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))

safe_predict <- function(fit, newdata, debias=FALSE) {
  out <- try(predict(fit, newdata=newdata, debias=debias), silent=TRUE)
  if (inherits(out, "try-error")) return(rep(NA_real_, nrow(newdata)))
  as.numeric(out)
}

## Count p_full (columns in model.matrix minus intercept) WITHOUT response 'y'
count_p_full <- function(fm, p) {
  rhs_fml <- stats::reformulate(
    attr(stats::delete.response(stats::terms(fm)), "term.labels")
  )
  tmp <- as.data.frame(matrix(rnorm(4 * p), nrow = 4))
  names(tmp) <- paste0("X", seq_len(p))
  mm <- model.matrix(rhs_fml, data = tmp)
  sum(colnames(mm) != "(Intercept)")
}

verify_no_leakage <- function(fit, train_df_used, test_df) {
  ok_disjoint <- !any(rownames(train_df_used) %in% rownames(test_df))
  raw_pred_train <- safe_predict(fit, train_df_used, debias=FALSE)
  y_train <- train_df_used$y
  ref <- lm(y_train ~ raw_pred_train)
  co_ref <- coef(ref)
  if (is.null(fit$debias_fit)) {
    return(list(ok=FALSE, ok_disjoint=ok_disjoint, coef_close=FALSE))
  }
  co_pkg <- coef(fit$debias_fit)
  same_len <- length(co_pkg) == length(co_ref)
  coef_close <- same_len && all(abs(co_pkg - co_ref) < 1e-8)
  list(ok = (ok_disjoint && coef_close), ok_disjoint=ok_disjoint, coef_close=coef_close)
}

build_formulas <- function(p) {
  vars <- paste0("X", seq_len(p))
  main_terms <- paste(vars, collapse=" + ")
  int_terms  <- paste0("(", main_terms, ")^2")
  sq_terms   <- paste(sprintf("I(%s^2)", vars), collapse=" + ")
  list(
    main_effects   = as.formula(paste0("y ~ ", main_terms)),
    main_plus_int  = as.formula(paste0("y ~ ", int_terms)),
    full_quadratic = as.formula(paste0("y ~ ", int_terms, " + ", sq_terms))
  )
}

## Sparse truth
truth_for_p <- function(p) {
  list(
    active_main  = seq_len(min(3, p)),
    active_pairs = if (p >= 3) rbind(c(1,2), c(2,3)) else matrix(numeric(), ncol=2),
    active_sq    = if (p >= 2) 2L else integer(0)
  )
}

## Equal scales; AR-ish correlation
gen_X <- function(n, p, rho=0.5) {
  Z <- matrix(rnorm(n*p), n, p)
  if (p >= 2 && abs(rho) > 0) {
    for (j in 2:p) {
      Z[, j] <- rho * scale(Z[, j-1], TRUE, TRUE) + sqrt(1-rho^2) * scale(Z[, j], TRUE, TRUE)
    }
  }
  X <- Z
  colnames(X) <- paste0("X", seq_len(p))
  as.data.frame(X)
}

## Data with targeted in-sample R^2
make_sparse_data_R2 <- function(n, p, rho, target_R2,
                                active_main, active_pairs, active_sq,
                                coef_main=2, coef_int=1.5, coef_sq=1.0, seed=NULL) {
  if (!is.null(seed)) set.seed(seed)
  X <- gen_X(n, p, rho)
  mu <- 1
  y_sig <- rep(mu, n)
  for (j in active_main) y_sig <- y_sig + coef_main * X[[paste0("X", j)]]
  if (length(active_sq)) for (j in active_sq) y_sig <- y_sig + coef_sq * (X[[paste0("X", j)]]^2)
  if (nrow(active_pairs) > 0) {
    for (k in seq_len(nrow(active_pairs))) {
      i <- active_pairs[k,1]; j <- active_pairs[k,2]
      y_sig <- y_sig + coef_int * (X[[paste0("X", i)]] * X[[paste0("X", j)]])
    }
  }
  vs <- var(y_sig); sd_sig <- sqrt(max(vs, .Machine$double.eps))
  sd_eps <- sd_sig * sqrt((1 - target_R2) / target_R2)
  y <- y_sig + rnorm(n, 0, sd_eps)
  df <- cbind.data.frame(y=y, X)
  df$y_signal <- y_sig
  rownames(df) <- sprintf("row%06d", seq_len(n))
  df
}

## ------------------ RUN ------------------
rows <- list(); rid <- 1
leaks <- list()
TOTAL_FITS <- sum(sapply(P_GRID, function(p) length(MODELS)*length(RHO_GRID)*length(R2_GRID)*REPS))
fit_counter <- 0L
if (PROGRESS) pb <- txtProgressBar(min = 0, max = TOTAL_FITS, style = 3)

for (p in P_GRID) {
  fms <- build_formulas(p)
  thr <- truth_for_p(p)

  for (model in MODELS) {
    fm <- fms[[model]]
    p_full <- count_p_full(fm, p)

    ## choose n_train to cross the boundary; clamp [10, 80]
    n_grid_raw <- round(c(0.5, 0.75, 1.0, 1.25, 1.75, 2.5) * p_full)
    n_grid <- sort(unique(pmin(pmax(n_grid_raw, 10), 80)))

    for (rho in RHO_GRID) {
      for (R2_tgt in R2_GRID) {
        for (rep in seq_len(REPS)) {
          n_tr <- sample(n_grid, 1L)
          n_te <- NTEST

          df <- make_sparse_data_R2(
            n = n_tr + n_te, p = p, rho = rho, target_R2 = R2_tgt,
            active_main  = thr$active_main,
            active_pairs = thr$active_pairs,
            active_sq    = thr$active_sq,
            seed = 1000 + 7*rep + p + round(100*R2_tgt) + round(100*rho)
          )

          idx <- sample(seq_len(nrow(df)), size = n_tr)
          train_df <- df[idx, ]
          test_df  <- df[-idx, ]

          ## PRE-FILTER complete cases for THIS formula (avoid NA crash)
          mf_all <- model.frame(fm, data = train_df, na.action = na.pass)
          keep <- complete.cases(mf_all)
          if (sum(keep) < 2) {
            warning(sprintf("Skipping: n_train=%d but <2 complete cases for %s (p=%d)", n_tr, model, p))
            next
          }
          train_used <- train_df[keep, , drop = FALSE]
          n_used <- nrow(train_used)

          ## realized true R^2
          r2_true <- function(d) var(d$y_signal)/var(d$y)
          R2_train_true <- r2_true(train_used)
          R2_test_true  <- r2_true(test_df)

          fit_counter <- fit_counter + 1L
          if (PROGRESS) setTxtProgressBar(pb, fit_counter)
          if (fit_counter %% max(1L, ECHO_EVERY) == 0L) {
            log_msg("Start #%d | p=%d (p_full=%d) | model=%s | rho=%.2f | R2=%.2f | n_train=%d (used=%d) | ratio=%.2f",
                    fit_counter, p, p_full, model, rho, R2_tgt, n_tr, n_used, n_used/p_full)
          }

          ## ---- Fit wAIC ----
          t0 <- proc.time()[3]
          fit_waic <- SVEMnet(
            formula = fm, data = train_used,
            nBoot = NBOOT,
            glmnet_alpha = c(1),
            weight_scheme = "SVEM",
            objective = "wAIC",
            standardize = TRUE
          )
          t1 <- proc.time()[3]

          ## ---- Fit wSSE ----
          t2 <- proc.time()[3]
          fit_wsse <- SVEMnet(
            formula = fm, data = train_used,
            nBoot = NBOOT,
            glmnet_alpha = c(1),
            weight_scheme = "SVEM",
            objective = "wSSE",
            standardize = TRUE
          )
          t3 <- proc.time()[3]

          ## Predictions
          pr_tr_raw_waic <- safe_predict(fit_waic, train_used, debias = FALSE)
          pr_te_raw_waic <- safe_predict(fit_waic, test_df,   debias = FALSE)
          pr_tr_deb_waic <- safe_predict(fit_waic, train_used, debias = TRUE)
          pr_te_deb_waic <- safe_predict(fit_waic, test_df,    debias = TRUE)

          pr_tr_raw_wsse <- safe_predict(fit_wsse, train_used, debias = FALSE)
          pr_te_raw_wsse <- safe_predict(fit_wsse, test_df,   debias = FALSE)
          pr_tr_deb_wsse <- safe_predict(fit_wsse, train_used, debias = TRUE)
          pr_te_deb_wsse <- safe_predict(fit_wsse, test_df,    debias = TRUE)

          ## RMSEs
          rmse_tr_raw_waic <- rmse(train_used$y, pr_tr_raw_waic)
          rmse_te_raw_waic <- rmse(test_df$y,    pr_te_raw_waic)
          rmse_tr_deb_waic <- rmse(train_used$y, pr_tr_deb_waic)
          rmse_te_deb_waic <- rmse(test_df$y,    pr_te_deb_waic)

          rmse_tr_raw_wsse <- rmse(train_used$y, pr_tr_raw_wsse)
          rmse_te_raw_wsse <- rmse(test_df$y,    pr_te_raw_wsse)
          rmse_tr_deb_wsse <- rmse(train_used$y, pr_tr_deb_wsse)
          rmse_te_deb_wsse <- rmse(test_df$y,    pr_te_deb_wsse)

          log_msg("Done  #%d | AIC raw=%.3f deb=%.3f | SSE raw=%.3f deb=%.3f | Δ_deb(raw): AIC=%.3f SSE=%.3f | %.1fs+%.1fs",
                  fit_counter, rmse_te_raw_waic, rmse_te_deb_waic,
                  rmse_te_raw_wsse, rmse_te_deb_wsse,
                  rmse_te_deb_waic - rmse_te_raw_waic,
                  rmse_te_deb_wsse - rmse_te_raw_wsse,
                  (t1 - t0), (t3 - t2))

          ## Save row
          rows[[rid]] <- data.frame(
            p = p, p_full = p_full, model = model, rho = rho, R2_target = R2_tgt,
            n_train = n_tr, n_train_used = n_used, n_test = n_te,
            ratio_n_over_p = n_used / p_full,
            above_boundary = as.integer(n_used > p_full),
            R2_true_train = R2_train_true, R2_true_test = R2_test_true,
            rmse_test_raw_waic = rmse_te_raw_waic,
            rmse_test_deb_waic = rmse_te_deb_waic,
            rmse_test_raw_wsse = rmse_te_raw_wsse,
            rmse_test_deb_wsse = rmse_te_deb_wsse,
            rmse_train_raw_waic = rmse_tr_raw_waic,
            rmse_train_deb_waic = rmse_tr_deb_waic,
            rmse_train_raw_wsse = rmse_tr_raw_wsse,
            rmse_train_deb_wsse = rmse_tr_deb_wsse,
            time_sec_waic = round(t1 - t0, 2),
            time_sec_wsse = round(t3 - t2, 2),
            stringsAsFactors = FALSE
          )
          rid <- rid + 1

          ## light leakage check
          if (rep <= 2) {
            chk <- verify_no_leakage(fit_waic, train_used, test_df)
            leaks[[length(leaks) + 1]] <- data.frame(
              p = p, model = model, rho = rho, R2_target = R2_tgt, n_train = n_tr, n_train_used = n_used,
              ok = chk$ok, ok_disjoint = chk$ok_disjoint, coef_close = chk$coef_close
            )
          }
        }
      }
    }
  }
}

if (!length(rows)) stop("No successful fits were recorded. Check logs/warnings above.")
res <- do.call(rbind, rows)
leak_df <- if (length(leaks)) do.call(rbind, leaks) else NULL

## Ratio bins for readability on plots/tables
res$ratio_bin <- cut(
  res$ratio_n_over_p,
  breaks = c(-Inf, 0.75, 1.0, 1.25, 1.75, Inf),
  labels = c("<0.75", "0.75–1.0", "1.0–1.25", "1.25–1.75", ">1.75"),
  right = TRUE
)

## ------------------ SUMMARY (with 4 win-rates) ------------------

## CI helper
ci95 <- function(x) {
  x <- x[is.finite(x)]
  n <- length(x); if (!n) return(c(mean=NA, lo=NA, hi=NA, n=0))
  m <- mean(x); se <- sd(x)/sqrt(n)
  c(mean=m, lo=m-1.96*se, hi=m+1.96*se, n=n)
}

## Grouping key
grp_vars <- c("model","R2_target","ratio_bin","above_boundary")

## Build summary per group using by() to avoid unpack glitches
summ <- do.call(rbind, by(res, INDICES = res[grp_vars], FUN = function(d){
  ## win indicators (row-wise)
  deb_win_waic <- as.numeric(d$rmse_test_deb_waic < d$rmse_test_raw_waic)
  deb_win_wsse <- as.numeric(d$rmse_test_deb_wsse < d$rmse_test_raw_wsse)
  deb_lose_waic <- as.numeric(d$rmse_test_deb_waic > d$rmse_test_raw_waic)
  deb_lose_wsse <- as.numeric(d$rmse_test_deb_wsse > d$rmse_test_raw_wsse)

  waic_win_raw <- as.numeric(d$rmse_test_raw_waic < d$rmse_test_raw_wsse)
  waic_win_deb <- as.numeric(d$rmse_test_deb_waic < d$rmse_test_deb_wsse)
  wsse_win_raw <- as.numeric(d$rmse_test_raw_wsse < d$rmse_test_raw_waic)
  wsse_win_deb <- as.numeric(d$rmse_test_deb_wsse < d$rmse_test_deb_waic)

  ## pooled winrates (redundant on purpose)
  winrate_debias_TRUE  <- mean(c(deb_win_waic,  deb_win_wsse),  na.rm = TRUE)
  winrate_debias_FALSE <- mean(c(deb_lose_waic, deb_lose_wsse), na.rm = TRUE)
  winrate_wAIC         <- mean(c(waic_win_raw,  waic_win_deb),  na.rm = TRUE)
  winrate_wSSE         <- mean(c(wsse_win_raw,  wsse_win_deb),  na.rm = TRUE)

  ## RMSE means + CI
  a_raw <- ci95(d$rmse_test_raw_waic)
  a_deb <- ci95(d$rmse_test_deb_waic)
  s_raw <- ci95(d$rmse_test_raw_wsse)
  s_deb <- ci95(d$rmse_test_deb_wsse)

  ## debias deltas (for context)
  dA <- ci95(d$rmse_test_deb_waic - d$rmse_test_raw_waic)
  dS <- ci95(d$rmse_test_deb_wsse - d$rmse_test_raw_wsse)

  ## realized true R2
  rtr <- ci95(d$R2_true_train); rte <- ci95(d$R2_true_test)

  ## timings
  tt  <- ci95(rowMeans(cbind(d$time_sec_waic, d$time_sec_wsse), na.rm = TRUE))

  data.frame(
    model = d$model[1],
    R2_target = d$R2_target[1],
    ratio_bin = d$ratio_bin[1],
    above_boundary = d$above_boundary[1],

    R2_true_train_mean = rtr["mean"], R2_true_test_mean = rte["mean"],

    rmse_raw_waic_mean = a_raw["mean"],  rmse_raw_waic_ci95_lo = a_raw["lo"],  rmse_raw_waic_ci95_hi = a_raw["hi"],
    rmse_deb_waic_mean = a_deb["mean"],  rmse_deb_waic_ci95_lo = a_deb["lo"],  rmse_deb_waic_ci95_hi = a_deb["hi"],
    rmse_raw_wsse_mean = s_raw["mean"],  rmse_raw_wsse_ci95_lo = s_raw["lo"],  rmse_raw_wsse_ci95_hi = s_raw["hi"],
    rmse_deb_wsse_mean = s_deb["mean"],  rmse_deb_wsse_ci95_lo = s_deb["lo"],  rmse_deb_wsse_ci95_hi = s_deb["hi"],

    ## the four requested winrates
    winrate_debias_TRUE  = winrate_debias_TRUE,
    winrate_debias_FALSE = winrate_debias_FALSE,
    winrate_wAIC         = winrate_wAIC,
    winrate_wSSE         = winrate_wSSE,

    ## debias deltas by objective (context)
    deb_raw_diff_AIC_mean = dA["mean"], deb_raw_diff_AIC_ci95_lo = dA["lo"], deb_raw_diff_AIC_ci95_hi = dA["hi"],
    deb_raw_diff_SSE_mean = dS["mean"], deb_raw_diff_SSE_ci95_lo = dS["lo"], deb_raw_diff_SSE_ci95_hi = dS["hi"],

    mean_time = tt["mean"],
    row.names = NULL,
    check.names = FALSE
  )
}))
rownames(summ) <- NULL
summ <- summ[order(summ$model, summ$R2_target, summ$ratio_bin, summ$above_boundary), ]

## Leakage summary
cat("\n---- Leakage summary (should be ~100%) ----\n")
if (!is.null(leak_df)) {
  print(aggregate(cbind(ok, ok_disjoint, coef_close) ~ model, data = leak_df,
                  FUN = function(x) mean(x, na.rm=TRUE)), row.names = FALSE)
} else {
  cat("No leakage samples collected.\n")
}

## Print compact summary columns
ord <- c("model","R2_target","ratio_bin","above_boundary",
         "R2_true_train_mean","R2_true_test_mean",
         "rmse_raw_waic_mean","rmse_deb_waic_mean",
         "rmse_raw_wsse_mean","rmse_deb_wsse_mean",
         "winrate_debias_TRUE","winrate_debias_FALSE",
         "winrate_wAIC","winrate_wSSE",
         "deb_raw_diff_AIC_mean","deb_raw_diff_AIC_ci95_lo","deb_raw_diff_AIC_ci95_hi",
         "deb_raw_diff_SSE_mean","deb_raw_diff_SSE_ci95_lo","deb_raw_diff_SSE_ci95_hi",
         "mean_time")
cat("\n================ SUMMARY (with 4 win-rates) ================\n")
print(summ[, ord], row.names = FALSE, digits = 4)

## ------------------ Plots (optional) ------------------

## Debias − Raw (negative helps) for each objective
plot_df <- within(res, {
  deb_raw_AIC  <- rmse_test_deb_waic - rmse_test_raw_waic
  deb_raw_SSE  <- rmse_test_deb_wsse - rmse_test_raw_wsse
})
ggA <- ggplot(plot_df, aes(ratio_n_over_p, deb_raw_AIC)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_point(alpha = 0.2) +
  geom_smooth(se = TRUE, method = "loess", span = 0.7) +
  facet_grid(R2_target ~ model) +
  labs(x = expression(paste(n[train[used]], " / ", p[full])),
       y = "Debias − Raw (TEST RMSE) under wAIC",
       title = "When does debiasing help? (wAIC)") +
  theme_bw()

ggS <- ggplot(plot_df, aes(ratio_n_over_p, deb_raw_SSE)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_point(alpha = 0.2) +
  geom_smooth(se = TRUE, method = "loess", span = 0.7) +
  facet_grid(R2_target ~ model) +
  labs(x = expression(paste(n[train[used]], " / ", p[full])),
       y = "Debias − Raw (TEST RMSE) under wSSE",
       title = "When does debiasing help? (wSSE)") +
  theme_bw()

print(ggA)
print(ggS)

## ------------------ Export (for your own modeling) ------------------
keep_cols <- c("p","p_full","model","rho","R2_target",
               "n_train","n_train_used","n_test",
               "ratio_n_over_p","above_boundary",
               "R2_true_train","R2_true_test",
               "rmse_test_raw_waic","rmse_test_deb_waic",
               "rmse_test_raw_wsse","rmse_test_deb_wsse",
               "time_sec_waic","time_sec_wsse")
res_export <- res[, keep_cols]
cat("\nFirst few rows of res_export:\n")
print(utils::head(res_export), row.names = FALSE)

## Write to CSV if you like:
## write.csv(res_export, "res_export.csv", row.names = FALSE)

## ---- Oneway-style plots: block on Row and show ANOM-like needles --------
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)

# If your results frame is called res_export, just do:
# res <- res_export

# 1) Long format with a Row (block) ID and log(RMSE)
df <- res_export %>%
  mutate(Row = row_number()) %>%                           # block label = run id
  pivot_longer(
    cols = c(rmse_test_deb_waic, rmse_test_deb_wsse,
             rmse_test_raw_waic, rmse_test_raw_wsse),
    names_to  = "settings",
    values_to = "rmse"
  ) %>%
  mutate(
    settings = factor(settings,
                      levels = c("rmse_test_deb_waic","rmse_test_deb_wsse",
                                 "rmse_test_raw_waic","rmse_test_raw_wsse")),
    lmrse    = log(pmax(rmse, .Machine$double.eps))
  )

# 2) Block-centering (remove Row mean, then add back overall mean so the scale is familiar)
overall_mean <- mean(df$lmrse, na.rm = TRUE)
df <- df %>%
  group_by(Row) %>%
  mutate(lmrse_block = lmrse - mean(lmrse, na.rm = TRUE) + overall_mean) %>%
  ungroup()

# 3) Fit two-way fixed effects model to get a pooled MSE (settings + block)
fit <- lm(lmrse_block ~ settings + factor(Row), data = df)
MSE <- summary(fit)$sigma^2
df_resid <- df.residual(fit)

# Effective sample size per setting (usually each setting appears once per Row)
n_i <- as.numeric(table(df$settings))
SE_i <- sqrt(MSE / n_i)            # SE for each setting mean (RCBD-style)
names(SE_i) <- names(table(df$settings))

# 4) Means per setting and an ANOM-like CI
alpha <- 0.05
tcrit <- qt(1 - alpha/2, df = df_resid)

means <- df %>%
  group_by(settings) %>%
  summarise(mean_lmrse = mean(lmrse_block, na.rm = TRUE), .groups = "drop") %>%
  mutate(
    SE   = SE_i[as.character(settings)],
    lo   = mean_lmrse - tcrit * SE,
    hi   = mean_lmrse + tcrit * SE
  )

# Grand mean and a constant decision band
grand <- mean(df$lmrse_block, na.rm = TRUE)
SE_bar <- mean(SE_i)                      # one band for the panel
UDL <- grand + tcrit * SE_bar
LDL <- grand - tcrit * SE_bar
band_df <- data.frame(LDL = LDL, UDL = UDL)

# 5) Panel 1: “Block Centered lmrse” (violin + dots)
p1 <- ggplot(df, aes(settings, lmrse_block)) +
  geom_violin(fill = "grey90", color = "grey40", scale = "width") +
  geom_jitter(width = 0.10, height = 0, alpha = 0.25, size = 1) +
  geom_hline(yintercept = grand, linetype = 2, color = "grey35") +
  labs(title = "Oneway Analysis of lmrse by settings",
       subtitle = "Block: Row (run ID) — block-centered values",
       y = "Block-Centered lmrse", x = "settings") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 6) Panel 2: “Analysis of Means” (lollipops + CI + decision band)
p2 <- ggplot(means, aes(settings, mean_lmrse)) +
  # decision band
  geom_rect(data = band_df,
            aes(xmin = -Inf, xmax = Inf, ymin = LDL, ymax = UDL),
            inherit.aes = FALSE, alpha = 0.12, fill = "skyblue") +
  geom_hline(yintercept = grand, color = "grey25") +
  # “needles”
  geom_segment(aes(x = settings, xend = settings, y = grand, yend = mean_lmrse),
               linewidth = 0.6, color = "black") +
  # point + CI
  geom_point(size = 3, color = "firebrick") +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.15) +
  annotate("text", x = Inf, y = UDL, hjust = 1.05, vjust = -0.4,
           label = sprintf("UDL=%.4f", UDL), size = 3.3) +
  annotate("text", x = Inf, y = grand, hjust = 1.05, vjust = -0.4,
           label = sprintf("Avg = %.4f", grand), size = 3.3) +
  annotate("text", x = Inf, y = LDL, hjust = 1.05, vjust = 1.2,
           label = sprintf("LDL=%.4f", LDL), size = 3.3) +
  labs(title = "Analysis of Means (ANOM-style)",
       subtitle = expression(paste("95% t-band; pooled MSE from ",
                                   "lm(lmrse_block ~ settings + block)")),
       y = "lmrse", x = "settings") +
  coord_cartesian(ylim = range(c(means$lo, means$hi, LDL, UDL))) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 7)
(p1 / p2) + plot_annotation(
  title = "Blocked Oneway (Row as Block) for lmrse vs. settings",
  theme = theme(plot.title = element_text(face = "bold"))
)

## ============================================================
## SVEM (wAIC / wAICc / wSSE, no debias) vs glmnet_with_cv (defaults)
## Focus: small training sizes; compare test RMSE head-to-head
## Parallelized with future.apply (5 workers)
## ============================================================

## ----- Packages -----
if (!requireNamespace("SVEMnet", quietly = TRUE)) install.packages("SVEMnet")
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("dplyr", quietly = TRUE)) install.packages("dplyr")
if (!requireNamespace("tidyr", quietly = TRUE)) install.packages("tidyr")
if (!requireNamespace("patchwork", quietly = TRUE)) install.packages("patchwork")
if (!requireNamespace("future.apply", quietly = TRUE)) install.packages("future.apply")
# Optional progress bar:
# if (!requireNamespace("progressr", quietly = TRUE)) install.packages("progressr")

library(SVEMnet)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(future.apply)
 library(progressr)  # uncomment if you want progress bars

## ----- Controls -----
set.seed(425)
REPS         <- 20
NBOOT        <- 200
R2_GRID      <- c(0.30, 0.50, 0.70, 0.90)
RHO_GRID     <- c(0, -0.5, 0.5, -0.9, 0.9)  # correlation among factors
P_GRID       <- c(3, 5, 7)
MODELS       <- c("main_effects","main_plus_int", "full_quadratic")
NTEST        <- 1000
DENSITY_GRID <- c(.05, .10, .20, .30, .50,.75)  # fraction of non-intercept terms active
N_WORKERS    <- 7

## ----- Helpers -----
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))

safe_predict_svem <- function(fit, newdata, agg = "mean") {
  out <- try(predict(fit, newdata = newdata, debias = FALSE, agg = agg), silent = TRUE)
  if (inherits(out, "try-error")) return(rep(NA_real_, nrow(newdata)))
  as.numeric(out)
}
safe_predict_cv <- function(fitcv, newdata) {
  out <- try(predict_cv(fitcv, newdata = newdata, debias = FALSE), silent = TRUE)
  if (inherits(out, "try-error")) return(rep(NA_real_, nrow(newdata)))
  as.numeric(out)
}

## Count p_full (columns in model matrix – intercept)
count_p_full <- function(fm, p) {
  rhs <- stats::reformulate(attr(stats::delete.response(stats::terms(fm)), "term.labels"))
  tmp <- as.data.frame(matrix(rnorm(4 * p), nrow = 4))
  names(tmp) <- paste0("X", seq_len(p))
  mm <- model.matrix(rhs, data = tmp)
  sum(colnames(mm) != "(Intercept)")
}

build_formulas <- function(p) {
  vars <- paste0("X", seq_len(p))
  main_terms <- paste(vars, collapse = " + ")
  int_terms  <- paste0("(", main_terms, ")^2")
  sq_terms   <- paste(sprintf("I(%s^2)", vars), collapse = " + ")
  list(
    main_effects   = as.formula(paste0("y ~ ", main_terms)),
    main_plus_int  = as.formula(paste0("y ~ ", int_terms)),
    full_quadratic = as.formula(paste0("y ~ ", int_terms, " + ", sq_terms))
  )
}

gen_X <- function(n, p, rho = 0.5) {
  Z <- matrix(rnorm(n * p), n, p)
  if (p >= 2 && abs(rho) > 0) {
    for (j in 2:p) {
      Z[, j] <- rho * scale(Z[, j - 1], TRUE, TRUE) +
        sqrt(1 - rho^2) * scale(Z[, j], TRUE, TRUE)
    }
  }
  X <- Z
  colnames(X) <- paste0("X", seq_len(p))
  as.data.frame(X)
}

## Density-driven truth generator using the selected formula fm
make_sparse_data_R2 <- function(n, p, rho, target_R2, fm, density = 0.2, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  Xdf <- gen_X(n, p, rho); names(Xdf) <- paste0("X", seq_len(p))
  terms_obj <- stats::delete.response(stats::terms(fm, data = transform(Xdf, y = 0)))
  mf  <- model.frame(terms_obj, data = Xdf)
  MM  <- model.matrix(terms_obj, mf)   # n x (M+1) with intercept at column 1

  M <- ncol(MM) - 1L
  if (M < 1L) stop("Formula produces no predictors beyond intercept.")
  n_active   <- max(1L, floor(density * M))
  active_idx <- sample.int(M, n_active, replace = FALSE) + 1L  # skip intercept

  beta <- numeric(M + 1L)
  beta[active_idx] <- rexp(n_active) - rexp(n_active)
  y_signal <- drop(MM %*% beta)

  sd_sig <- sqrt(max(var(y_signal), .Machine$double.eps))
  sd_eps <- sd_sig * sqrt((1 - target_R2) / target_R2)
  y <- y_signal + rnorm(n, 0, sd_eps)

  out <- cbind.data.frame(y = y, Xdf)
  out$y_signal <- y_signal
  rownames(out) <- sprintf("row%06d", seq_len(n))
  out
}

## Convenience wrapper to fit SVEM for a given objective (namespaced; robust)
fit_and_predict_svem <- function(fm, train_used, test_df, objective, nBoot = NBOOT, agg = "mean") {
  t0 <- proc.time()[3]
  mod <- try(SVEMnet::SVEMnet(
    formula = fm, data = train_used,
    nBoot = nBoot,
    glmnet_alpha = c(0,.5,1),
    weight_scheme = "SVEM",
    objective = objective,
    standardize = TRUE
  ), silent = TRUE)
  elapsed <- round(proc.time()[3] - t0, 2)
  if (inherits(mod, "try-error")) {
    preds <- rep(NA_real_, nrow(test_df))
  } else {
    preds <- safe_predict_svem(mod, test_df, agg = agg)
  }
  list(pred = preds, time = elapsed)
}

## ----- Build scenario grid -----
grid <- expand.grid(
  p       = P_GRID,
  model   = MODELS,
  rho     = RHO_GRID,
  R2_tgt  = R2_GRID,
  dens    = DENSITY_GRID,
  rep     = seq_len(REPS),
  stringsAsFactors = FALSE
)

## ----- Parallel plan (5 workers) -----
old_plan <- future::plan()
on.exit(future::plan(old_plan), add = TRUE)
future::plan(future::multisession, workers = N_WORKERS)

## Optional progress bar:
# handlers(global = TRUE)
# with_progress({
#   pbar <- progressor(along = seq_len(nrow(grid)))
#   ... call future_lapply with on.exit(pbar()) ...
# })

## ----- One scenario run -----
run_one <- function(row) {
  p      <- row$p
  model  <- row$model
  rho    <- row$rho
  R2_tgt <- row$R2_tgt
  dens   <- row$dens
  repi   <- row$rep

  fms    <- build_formulas(p)
  fm     <- fms[[model]]
  p_full <- count_p_full(fm, p)

  # n around the boundary; clamp to [12, 40]
  n_grid <- sort(unique(pmin(pmax(p_full + c(-12,-8,-6,-4,-2,0,2,4,6,8,12), 12), 40)))

  # per-row deterministic seed (so sampling is reproducible too)
  seed_i <- 3000 +
    19 * repi + 100 * p +
    1000 * match(model, MODELS) +
    round(1000 * R2_tgt) +
    round(1000 * rho) +
    round(1000 * dens)
  set.seed(seed_i)

  n_tr <- sample(n_grid, 1L)
  n_te <- NTEST

  df <- make_sparse_data_R2(
    n = n_tr + n_te, p = p, rho = rho, target_R2 = R2_tgt,
    fm = fm, density = dens, seed = seed_i + 1L
  )

  set.seed(seed_i + 2L)
  idx      <- sample(seq_len(nrow(df)), size = n_tr)
  train_df <- df[idx, ]
  test_df  <- df[-idx, ]

  keep <- complete.cases(model.frame(fm, data = train_df, na.action = stats::na.pass))
  if (sum(keep) < 2) return(NULL)
  train_used <- train_df[keep, , drop = FALSE]
  n_used     <- nrow(train_used)

  r2_true <- function(d) var(d$y_signal) / var(d$y)
  R2_train_true <- r2_true(train_used)
  R2_test_true  <- r2_true(test_df)

  # SVEM variants (alpha=1; no debias)
  sv_waic  <- fit_and_predict_svem(fm, train_used, test_df, "wAIC",  nBoot = NBOOT)
  sv_waicc <- fit_and_predict_svem(fm, train_used, test_df, "wAICc", nBoot = NBOOT)
  sv_wsse  <- fit_and_predict_svem(fm, train_used, test_df, "wSSE",  nBoot = NBOOT)

  # glmnet_with_cv (alpha=1; no debias)
  t1 <- proc.time()[3]
  fit_cv <- glmnet_with_cv(
    formula = fm, data = train_used,
    glmnet_alpha = c(0,.5,1),
    standardize = TRUE
  )
  pr_te_cv <- safe_predict_cv(fit_cv, test_df)
  time_cv  <- round(proc.time()[3] - t1, 2)

  data.frame(
    p = p, model = model, rho = rho, R2_target = R2_tgt, density = dens,
    p_full = p_full, n_train = n_tr, n_train_used = n_used, n_test = n_te,
    n_train_minus_p_full = n_tr - p_full,
    ratio_n_over_p = n_used / p_full,
    above_boundary = as.integer(n_used > p_full),
    R2_true_train = R2_train_true, R2_true_test = R2_test_true,
    rmse_test_svem_waic  = rmse(test_df$y, sv_waic$pred),
    rmse_test_svem_waicc = rmse(test_df$y, sv_waicc$pred),
    rmse_test_svem_wsse  = rmse(test_df$y, sv_wsse$pred),
    rmse_test_glmnet_cv  = rmse(test_df$y, pr_te_cv),
    time_sec_svem_waic   = sv_waic$time,
    time_sec_svem_waicc  = sv_waicc$time,
    time_sec_svem_wsse   = sv_wsse$time,
    time_sec_cv          = time_cv,
    stringsAsFactors = FALSE
  )
}

## ----- Run in parallel -----
res_list <- future.apply::future_lapply(
  seq_len(nrow(grid)),
  function(i) {
    # if using progressr:
    # on.exit(pbar(), add = TRUE)
    run_one(grid[i, , drop = FALSE])
  },
  future.seed = TRUE
)

res <- dplyr::bind_rows(Filter(Negate(is.null), res_list))
stopifnot(nrow(res) > 0)

## ---------- Summary & win-rates (robust, single block) ----------
stopifnot(all(c("rmse_test_svem_waic","rmse_test_svem_waicc","rmse_test_svem_wsse","rmse_test_glmnet_cv") %in% names(res)))

# Long format for plotting and summaries
res_long <- res %>%
  mutate(
    Row           = dplyr::row_number(),
    lmrse_waic    = log(pmax(rmse_test_svem_waic,  .Machine$double.eps)),
    lmrse_waicc   = log(pmax(rmse_test_svem_waicc, .Machine$double.eps)),
    lmrse_wsse    = log(pmax(rmse_test_svem_wsse,  .Machine$double.eps)),
    lmrse_cv      = log(pmax(rmse_test_glmnet_cv,  .Machine$double.eps))
  ) %>%
  tidyr::pivot_longer(
    cols = c(lmrse_waic, lmrse_waicc, lmrse_wsse, lmrse_cv),
    names_to = "metric", values_to = "lmrse"
  ) %>%
  mutate(
    settings = dplyr::recode(metric,
                             lmrse_waic  = "SVEM_wAIC",
                             lmrse_waicc = "SVEM_wAICc",
                             lmrse_wsse  = "SVEM_wSSE",
                             lmrse_cv    = "glmnet_cv"),
    settings = factor(settings, levels = c("SVEM_wAIC","SVEM_wAICc","SVEM_wSSE","glmnet_cv"))
  )

# Head-to-head win-rates vs glmnet_cv (paired by row)
win_rates_vs_cv <- res %>%
  summarise(
    SVEM_wAIC  = mean(rmse_test_svem_waic  < rmse_test_glmnet_cv, na.rm = TRUE),
    SVEM_wAICc = mean(rmse_test_svem_waicc < rmse_test_glmnet_cv, na.rm = TRUE),
    SVEM_wSSE  = mean(rmse_test_svem_wsse  < rmse_test_glmnet_cv, na.rm = TRUE)
  ) %>%
  tidyr::pivot_longer(everything(), names_to = "settings", values_to = "winrate_vs_cv")

# Overall summary (log RMSE)
summ <- res_long %>%
  group_by(settings) %>%
  summarise(
    mean_lmrse = mean(lmrse, na.rm = TRUE),
    sd_lmrse   = sd(lmrse,   na.rm = TRUE),
    n          = dplyr::n(),
    se         = sd_lmrse / sqrt(pmax(n, 1)),
    ci_lo      = mean_lmrse - 1.96 * se,
    ci_hi      = mean_lmrse + 1.96 * se,
    .groups    = "drop"
  ) %>%
  left_join(win_rates_vs_cv, by = "settings")

cat("\n================ SUMMARY (log RMSE) ================\n")
print(summ, row.names = FALSE, digits = 4)

## ---------- By-boundary summaries (attach boundary to long) ----------
if ("above_boundary" %in% names(res)) {
  bound_key <- res %>%
    mutate(Row = dplyr::row_number()) %>%
    select(Row, above_boundary)
} else if (all(c("n_train_used","p_full") %in% names(res))) {
  bound_key <- res %>%
    mutate(Row = dplyr::row_number(),
           above_boundary = as.integer(n_train_used > p_full)) %>%
    select(Row, above_boundary)
} else if ("ratio_n_over_p" %in% names(res)) {
  bound_key <- res %>%
    mutate(Row = dplyr::row_number(),
           above_boundary = as.integer(ratio_n_over_p > 1)) %>%
    select(Row, above_boundary)
} else {
  bound_key <- res %>%
    mutate(Row = dplyr::row_number(), above_boundary = NA_integer_) %>%
    select(Row, above_boundary)
}

res_long <- res_long %>% left_join(bound_key, by = "Row")

if ("above_boundary.x" %in% names(res_long) || "above_boundary.y" %in% names(res_long)) {
  res_long <- res_long %>%
    dplyr::mutate(
      above_boundary = dplyr::coalesce(
        if ("above_boundary.y" %in% names(.)) .data$above_boundary.y else NA_integer_,
        if ("above_boundary.x" %in% names(.)) .data$above_boundary.x else NA_integer_
      )
    ) %>%
    dplyr::select(-dplyr::any_of(c("above_boundary.x","above_boundary.y")))
}

# Means of log RMSE by boundary & setting
summ_by_boundary <- res_long %>%
  group_by(above_boundary, settings) %>%
  summarise(
    mean_lmrse = mean(lmrse, na.rm = TRUE),
    sd_lmrse   = sd(lmrse,   na.rm = TRUE),
    n          = dplyr::n(),
    .groups    = "drop"
  )

cat("\n---- Means by boundary (n_train_used > p_full) ----\n")
print(summ_by_boundary, row.names = FALSE, digits = 4)

# Boundary-specific win-rate for each SVEM variant vs glmnet
win_by_boundary <- res %>%
  mutate(
    above_boundary = if ("above_boundary" %in% names(.)) above_boundary
    else if (all(c("n_train_used","p_full") %in% names(.))) as.integer(n_train_used > p_full)
    else if ("ratio_n_over_p" %in% names(.)) as.integer(ratio_n_over_p > 1)
    else NA_integer_,
    win_waic  = rmse_test_svem_waic  < rmse_test_glmnet_cv,
    win_waicc = rmse_test_svem_waicc < rmse_test_glmnet_cv,
    win_wsse  = rmse_test_svem_wsse  < rmse_test_glmnet_cv
  ) %>%
  group_by(above_boundary) %>%
  summarise(
    winrate_waic  = mean(win_waic,  na.rm = TRUE),
    winrate_waicc = mean(win_waicc, na.rm = TRUE),
    winrate_wsse  = mean(win_wsse,  na.rm = TRUE),
    n             = sum(!is.na(win_waic) | !is.na(win_waicc) | !is.na(win_wsse)),
    .groups       = "drop"
  )

cat("\n---- Win-rate vs glmnet (by boundary) ----\n")
print(win_by_boundary, row.names = FALSE, digits = 4)

## ========== ANOM-style, row-blocked plot ==========
res_long_anom <- res_long %>%
  dplyr::mutate(Row = as.factor(Row)) %>%
  dplyr::filter(is.finite(lmrse))

fit_anom <- lm(lmrse ~ settings + Row, data = res_long_anom)
aov_tbl  <- anova(fit_anom)

MSE    <- aov_tbl[["Mean Sq"]][nrow(aov_tbl)]
df_res <- aov_tbl[["Df"]][nrow(aov_tbl)]

t_methods <- nlevels(res_long_anom$settings)
b_blocks  <- nlevels(res_long_anom$Row)
grand_mu  <- mean(res_long_anom$lmrse, na.rm = TRUE)

# Var(ȳ_i. − ȳ..) = σ^2 * (t − 1) / (t * b) under RBD
se_group <- sqrt(MSE * (t_methods - 1) / (t_methods * b_blocks))

alpha <- 0.05
crit  <- qt(1 - alpha / (2 * t_methods), df = df_res)
UCL   <- grand_mu + crit * se_group
LCL   <- grand_mu - crit * se_group

means_df <- res_long_anom %>%
  dplyr::group_by(settings) %>%
  dplyr::summarise(mean_lmrse = mean(lmrse, na.rm = TRUE), .groups = "drop") %>%
  dplyr::mutate(flag = dplyr::case_when(
    mean_lmrse > UCL ~ "Above UCL",
    mean_lmrse < LCL ~ "Below LCL",
    TRUE             ~ "Within Limits"
  ))

p_anom <- ggplot(means_df, aes(x = settings, y = mean_lmrse)) +
  geom_hline(yintercept = grand_mu, linetype = 2) +
  geom_hline(yintercept = UCL,      linetype = 3) +
  geom_hline(yintercept = LCL,      linetype = 3) +
  geom_segment(aes(xend = settings, y = grand_mu, yend = mean_lmrse), linewidth = 1) +
  geom_point(aes(color = flag), size = 3) +
  scale_color_manual(values = c("Within Limits" = "black",
                                "Above UCL"     = "red",
                                "Below LCL"     = "red")) +
  labs(
    title = "Blocked ANOM-style plot for lmrse (Row = block)",
    subtitle = sprintf("Grand mean = %.3f | Limits = [%.3f, %.3f] | df = %d",
                       grand_mu, LCL, UCL, df_res),
    x = NULL, y = "Mean lmrse", color = NULL
  ) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"))

p_pairs <- ggplot(res_long_anom, aes(x = settings, y = lmrse, group = Row)) +
  geom_line(alpha = 0.25) +
  geom_point(alpha = 0.6, position = position_dodge(width = 0.05)) +
  stat_summary(fun = mean, geom = "point", size = 3, color = "black") +
  labs(title = "Paired runs by Row", x = NULL, y = "lmrse") +
  theme_bw()

(p_anom / p_pairs)

## ---------- Geometric-mean RMSE ratio (SVEM wAIC vs CV) ----------
eps <- .Machine$double.eps
rmse_ratio <- with(res, pmax(rmse_test_svem_waic, eps) / pmax(rmse_test_glmnet_cv, eps))
logdiff    <- log(rmse_ratio)

gm_ratio <- exp(mean(logdiff, na.rm = TRUE))           # primary
q75      <- quantile(rmse_ratio, 0.75, na.rm = TRUE)   # tail
q90      <- quantile(rmse_ratio, 0.90, na.rm = TRUE)
winrate  <- mean(rmse_ratio < 1, na.rm = TRUE)

# simple bootstrap CI for GM ratio
set.seed(123)
B <- 2000
idx <- replicate(B, sample.int(length(logdiff), replace = TRUE))
gm_boot <- apply(idx, 2, function(ii) exp(mean(logdiff[ii], na.rm = TRUE)))
gm_ci   <- quantile(gm_boot, c(0.025, 0.975), na.rm = TRUE)

cat(sprintf(
  "\nPaired geometric-mean RMSE ratio (SVEM/CV): %.3f  (95%% CI %.3f – %.3f)\n",
  gm_ratio, gm_ci[1], gm_ci[2]
))
cat(sprintf("Win rate (SVEM better): %.1f%%\n", 100*winrate))
cat(sprintf("Tail ratios: 75th=%.3f, 90th=%.3f\n", q75, q90))

# Trim top/bottom 5% of logdiff before averaging
trim <- function(x, p=0.05) x[x >= quantile(x, p, na.rm=TRUE) & x <= quantile(x, 1-p, na.rm=TRUE)]
gm_ratio_trim <- exp(mean(trim(logdiff), na.rm = TRUE))
cat(sprintf("Trimmed GM ratio (5%%/5%%): %.3f\n", gm_ratio_trim))