# Combined tests from unit/
# Auto-generated by workspace refactor script

# ===== BEGIN unit/test-dispatch-separation.R =====
test_that("scalar vs vector dispatch separation is enforced", {

  spec <- list(
    meta = list(backend = "sb", kernel = "gamma", GPD = FALSE),
    dispatch = list(backend = "sb", GPD = FALSE),
    kernel = list(key = "gamma")
  )

  scalar <- .get_dispatch_scalar(spec)
  expect_false(isTRUE(attr(scalar$d, "vectorized_wrapper")))
  expect_false(isTRUE(attr(scalar$p, "vectorized_wrapper")))
  expect_false(isTRUE(attr(scalar$q, "vectorized_wrapper")))
  expect_false(isTRUE(attr(scalar$r, "vectorized_wrapper")))

  vec <- .get_dispatch(spec)
  expect_true(isTRUE(attr(vec$d, "vectorized_wrapper")))
  expect_true(isTRUE(attr(vec$p, "vectorized_wrapper")))
  expect_true(isTRUE(attr(vec$q, "vectorized_wrapper")))
  expect_true(isTRUE(attr(vec$r, "vectorized_wrapper")))
})

# ===== END unit/test-dispatch-separation.R =====

# ===== BEGIN unit/test-distributions.R =====
# test-distributions.R
# Consolidated per-kernel distribution integration tests (Tier B)
# Merged from: test-amoroso.R, test-cauchy.R, test-gamma.R, test-invgauss.R,
#             test-laplace.R, test-lognormal.R, test-normal.R

if (!exists(".run_predict_case")) {
  helper_path <- file.path("tests", "testthat", "helper-03-predict-helpers.R")
  if (file.exists(helper_path)) source(helper_path)
}

# =============================================================================
# Amoroso kernel
# =============================================================================
test_that("SB amoroso: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB amoroso: no GPD, no X", "amoroso", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB amoroso: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB amoroso: GPD, no X", "amoroso", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB amoroso: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB amoroso: no GPD, with X", "amoroso", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB amoroso: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB amoroso: GPD, with X", "amoroso", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP amoroso: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP amoroso: no GPD, no X", "amoroso", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP amoroso: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP amoroso: GPD, no X", "amoroso", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP amoroso: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP amoroso: no GPD, with X", "amoroso", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP amoroso: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP amoroso: GPD, with X", "amoroso", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# =============================================================================
# Cauchy kernel
# =============================================================================
test_that("SB cauchy: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB cauchy: no GPD, no X", "cauchy", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB cauchy: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB cauchy: GPD, no X", "cauchy", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB cauchy: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB cauchy: no GPD, with X", "cauchy", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB cauchy: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB cauchy: GPD, with X", "cauchy", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP cauchy: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP cauchy: no GPD, no X", "cauchy", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP cauchy: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP cauchy: GPD, no X", "cauchy", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP cauchy: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP cauchy: no GPD, with X", "cauchy", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP cauchy: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP cauchy: GPD, with X", "cauchy", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# =============================================================================
# Gamma kernel
# =============================================================================
test_that("SB gamma: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB gamma: no GPD, no X", "gamma", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB gamma: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB gamma: GPD, no X", "gamma", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB gamma: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB gamma: no GPD, with X", "gamma", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB gamma: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB gamma: GPD, with X", "gamma", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP gamma: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP gamma: no GPD, no X", "gamma", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP gamma: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP gamma: GPD, no X", "gamma", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP gamma: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP gamma: no GPD, with X", "gamma", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP gamma: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP gamma: GPD, with X", "gamma", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# =============================================================================
# InvGauss kernel
# =============================================================================
test_that("SB invgauss: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB invgauss: no GPD, no X", "invgauss", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB invgauss: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB invgauss: GPD, no X", "invgauss", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB invgauss: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB invgauss: no GPD, with X", "invgauss", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB invgauss: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB invgauss: GPD, with X", "invgauss", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP invgauss: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP invgauss: no GPD, no X", "invgauss", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP invgauss: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP invgauss: GPD, no X", "invgauss", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP invgauss: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP invgauss: no GPD, with X", "invgauss", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP invgauss: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP invgauss: GPD, with X", "invgauss", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# =============================================================================
# Laplace kernel
# =============================================================================
test_that("SB laplace: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB laplace: no GPD, no X", "laplace", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB laplace: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB laplace: GPD, no X", "laplace", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB laplace: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB laplace: no GPD, with X", "laplace", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB laplace: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB laplace: GPD, with X", "laplace", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP laplace: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP laplace: no GPD, no X", "laplace", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP laplace: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP laplace: GPD, no X", "laplace", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP laplace: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP laplace: no GPD, with X", "laplace", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP laplace: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP laplace: GPD, with X", "laplace", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# =============================================================================
# Lognormal kernel
# =============================================================================
test_that("SB lognormal: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB lognormal: no GPD, no X", "lognormal", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB lognormal: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB lognormal: GPD, no X", "lognormal", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB lognormal: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB lognormal: no GPD, with X", "lognormal", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB lognormal: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB lognormal: GPD, with X", "lognormal", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP lognormal: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP lognormal: no GPD, no X", "lognormal", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP lognormal: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP lognormal: GPD, no X", "lognormal", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP lognormal: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP lognormal: no GPD, with X", "lognormal", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP lognormal: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP lognormal: GPD, with X", "lognormal", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# =============================================================================
# Normal kernel
# =============================================================================
test_that("SB normal: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB normal: no GPD, no X", "normal", "sb", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB normal: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB normal: GPD, no X", "normal", "sb", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB normal: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB normal: no GPD, with X", "normal", "sb", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("SB normal: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("SB normal: GPD, with X", "normal", "sb", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP normal: no GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP normal: no GPD, no X", "normal", "crp", FALSE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP normal: GPD, no X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP normal: GPD, no X", "normal", "crp", TRUE, FALSE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP normal: no GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP normal: no GPD, with X", "normal", "crp", FALSE, TRUE)
  expect_true(res$ok, info = res$msg)
})

test_that("CRP normal: GPD, with X", {
  skip_if_not_test_level("ci")
  res <- .run_predict_case("CRP normal: GPD, with X", "normal", "crp", TRUE, TRUE)
  expect_true(res$ok, info = res$msg)
})

# ===== END unit/test-distributions.R =====

# ===== BEGIN unit/test-global-contracts.R =====
# Tests for global contract functions (00-global-contracts.R)

# ============================================================================
# is_allowed_kernel() tests
# ============================================================================

test_that("is_allowed_kernel() returns TRUE for valid kernels", {
  is_allowed_kernel <- getFromNamespace("is_allowed_kernel", "CausalMixGPD")

  expect_true(is_allowed_kernel("normal"))
  expect_true(is_allowed_kernel("gamma"))
  expect_true(is_allowed_kernel("lognormal"))
  expect_true(is_allowed_kernel("invgauss"))
  expect_true(is_allowed_kernel("laplace"))
  expect_true(is_allowed_kernel("cauchy"))
  expect_true(is_allowed_kernel("amoroso"))
})

test_that("is_allowed_kernel() returns FALSE for invalid kernels", {
  is_allowed_kernel <- getFromNamespace("is_allowed_kernel", "CausalMixGPD")

  expect_false(is_allowed_kernel("invalid"))
  expect_false(is_allowed_kernel("weibull"))
  expect_false(is_allowed_kernel("exponential"))
  expect_false(is_allowed_kernel(""))
  expect_false(is_allowed_kernel("Normal"))  # Case sensitive
})

# ============================================================================
# check_gpd_contract() tests
# ============================================================================

test_that("check_gpd_contract() errors when GPD=TRUE with cauchy kernel", {
  check_gpd_contract <- getFromNamespace("check_gpd_contract", "CausalMixGPD")

  expect_error(
    check_gpd_contract(GPD = TRUE, kernel = "cauchy"),
    "Cauchy kernels are never paired with GPD"
  )
})

test_that("check_gpd_contract() returns invisibly for valid combinations", {
  check_gpd_contract <- getFromNamespace("check_gpd_contract", "CausalMixGPD")

  # GPD = TRUE with non-cauchy kernels should work
  expect_invisible(check_gpd_contract(GPD = TRUE, kernel = "normal"))
  expect_invisible(check_gpd_contract(GPD = TRUE, kernel = "gamma"))
  expect_invisible(check_gpd_contract(GPD = TRUE, kernel = "lognormal"))
  expect_invisible(check_gpd_contract(GPD = TRUE, kernel = "invgauss"))
  expect_invisible(check_gpd_contract(GPD = TRUE, kernel = "laplace"))
  expect_invisible(check_gpd_contract(GPD = TRUE, kernel = "amoroso"))

  # GPD = FALSE should always work
  expect_invisible(check_gpd_contract(GPD = FALSE, kernel = "cauchy"))
  expect_invisible(check_gpd_contract(GPD = FALSE, kernel = "normal"))
})

test_that("check_gpd_contract() returns NULL invisibly", {
  check_gpd_contract <- getFromNamespace("check_gpd_contract", "CausalMixGPD")

  result <- check_gpd_contract(GPD = TRUE, kernel = "normal")
  expect_null(result)
})

# ============================================================================
# allowed_* constant tests
# ============================================================================

test_that("allowed_backends constant is correct", {
  allowed_backends <- getFromNamespace("allowed_backends", "CausalMixGPD")

  expect_equal(allowed_backends, c("crp", "sb", "spliced"))
})

test_that("allowed_kernels constant contains all 7 kernels", {
  allowed_kernels <- getFromNamespace("allowed_kernels", "CausalMixGPD")

  expected <- c("gamma", "lognormal", "invgauss", "normal", "laplace", "cauchy", "amoroso")
  expect_equal(sort(allowed_kernels), sort(expected))
})

test_that("positive_support_kernels constant is correct", {
  positive_support_kernels <- getFromNamespace("positive_support_kernels", "CausalMixGPD")

  expect_true("gamma" %in% positive_support_kernels)
  expect_true("lognormal" %in% positive_support_kernels)
  expect_true("invgauss" %in% positive_support_kernels)
  expect_true("amoroso" %in% positive_support_kernels)
  expect_false("normal" %in% positive_support_kernels)
  expect_false("laplace" %in% positive_support_kernels)
  expect_false("cauchy" %in% positive_support_kernels)
})

test_that("real_support_kernels constant is correct", {
  real_support_kernels <- getFromNamespace("real_support_kernels", "CausalMixGPD")

  expect_true("normal" %in% real_support_kernels)
  expect_true("laplace" %in% real_support_kernels)
  expect_true("cauchy" %in% real_support_kernels)
  expect_false("gamma" %in% real_support_kernels)
  expect_false("lognormal" %in% real_support_kernels)
  expect_false("invgauss" %in% real_support_kernels)
  expect_false("amoroso" %in% real_support_kernels)
})

# ===== END unit/test-global-contracts.R =====

# ===== BEGIN unit/test-glue-validity.R =====
# test-glue-validity.R
# CI-safe smoke test for glue diagnostics.

test_that("check_glue_validity passes on fixture (short grid, few draws)", {
  testthat::skip_on_cran()
  skip_if_not_test_level("ci")

  fit <- .load_fixture("fit_gpd_small.rds")
  small_grid <- seq(0.1, 5, length.out = 20)

  out <- check_glue_validity(fit, n_draws = 5L, grid = small_grid, check_continuity = FALSE)
  expect_true(is.list(out))
  expect_true(is.list(out$pass))
  expect_true(isTRUE(out$pass$cdf_range))
  expect_true(isTRUE(out$pass$cdf_monotone))
  expect_true(isTRUE(out$pass$density_nonneg))
})

# ===== END unit/test-glue-validity.R =====

# ===== BEGIN unit/test-hpd-intervals.R =====
# Test HPD interval functionality across all S3 functions
# Tests that interval = "hpd" works correctly and produces valid intervals

# ==============================================================================
# Test .compute_interval() helper function
# ==============================================================================

test_that(".compute_interval computes equal-tailed intervals correctly", {
  set.seed(42)
  draws <- rnorm(1000)

  iv <- .compute_interval(draws, level = 0.95, type = "credible")

  expect_named(iv, c("lower", "upper"))
  expect_true(iv["lower"] < iv["upper"])
  expect_true(iv["lower"] < 0)  # For standard normal
  expect_true(iv["upper"] > 0)

  # Check approximate quantiles
  q <- quantile(draws, probs = c(0.025, 0.975))
  expect_equal(unname(iv["lower"]), unname(q[1]))
  expect_equal(unname(iv["upper"]), unname(q[2]))
})


test_that(".compute_interval computes HPD intervals correctly", {
  set.seed(42)
  draws <- rnorm(1000)

  iv <- .compute_interval(draws, level = 0.95, type = "hpd")

  expect_named(iv, c("lower", "upper"))
  expect_true(iv["lower"] < iv["upper"])
  expect_true(iv["lower"] < 0)  # For standard normal
  expect_true(iv["upper"] > 0)

  # HPD should approximately match coda::HPDinterval
  hpd_coda <- coda::HPDinterval(coda::as.mcmc(draws), prob = 0.95)
  expect_equal(unname(iv["lower"]), unname(hpd_coda[1, "lower"]), tolerance = 0.001)
  expect_equal(unname(iv["upper"]), unname(hpd_coda[1, "upper"]), tolerance = 0.001)
})


test_that(".compute_interval handles skewed distributions differently", {
  set.seed(42)
  # Skewed distribution - HPD should give shorter interval
  draws <- rexp(1000)

  iv_credible <- .compute_interval(draws, level = 0.95, type = "credible")
  iv_hpd <- .compute_interval(draws, level = 0.95, type = "hpd")

  width_credible <- iv_credible["upper"] - iv_credible["lower"]
  width_hpd <- iv_hpd["upper"] - iv_hpd["lower"]

  # HPD interval should be narrower for skewed distributions
  expect_true(width_hpd <= width_credible + 0.01)
})


test_that(".compute_interval handles edge cases", {
  # Empty or NA draws
  iv <- .compute_interval(c(NA, NA, NA), level = 0.95, type = "credible")
  expect_true(is.na(iv["lower"]))
  expect_true(is.na(iv["upper"]))

  # Single value
  iv <- .compute_interval(c(1), level = 0.95, type = "credible")
  expect_true(is.na(iv["lower"]))
  expect_true(is.na(iv["upper"]))
})


# ==============================================================================
# Test predict.mixgpd_fit with interval parameter
# ==============================================================================

test_that("predict.mixgpd_fit supports unified interval parameter", {
  skip_if_not_test_level("ci")

  set.seed(42)
  y <- abs(rnorm(40)) + 0.2

  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "gamma",
    GPD = FALSE,
    components = 6,
    mcmc = list(niter = 200, nburnin = 50, nchains = 1)
  )

  nullfile <- if (.Platform$OS.type == "windows") "NUL" else "/dev/null"
  utils::capture.output(
    fit <- run_mcmc_bundle_manual(bundle, show_progress = FALSE),
    file = nullfile
  )

  # Test quantile prediction with both interval types
  pred_credible <- predict(fit, type = "quantile", index = c(0.5, 0.9),
                           interval = "credible")
  pred_hpd <- predict(fit, type = "quantile", index = c(0.5, 0.9),
                      interval = "hpd")
  pred_none <- predict(fit, type = "quantile", index = c(0.5, 0.9),
                       interval = NULL)

  expect_s3_class(pred_credible, "mixgpd_predict")
  expect_s3_class(pred_hpd, "mixgpd_predict")
  expect_s3_class(pred_none, "mixgpd_predict")

  # Both credible and hpd should have valid intervals
  expect_true(all(pred_credible$fit$lower <= pred_credible$fit$estimate))
  expect_true(all(pred_credible$fit$estimate <= pred_credible$fit$upper))
  expect_true(all(pred_hpd$fit$lower <= pred_hpd$fit$estimate))
  expect_true(all(pred_hpd$fit$estimate <= pred_hpd$fit$upper))

  # NULL interval should have NA for lower/upper
  expect_true(all(is.na(pred_none$fit$lower)))
  expect_true(all(is.na(pred_none$fit$upper)))

  # Default should be credible
  pred_default <- predict(fit, type = "quantile", index = c(0.5, 0.9))
  expect_equal(pred_default$fit$lower, pred_credible$fit$lower)
  expect_equal(pred_default$fit$upper, pred_credible$fit$upper)
})


test_that("predict.mixgpd_fit mean prediction works with HPD", {
  skip_if_not_test_level("ci")

  set.seed(42)
  y <- abs(rnorm(40)) + 0.2

  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "gamma",
    GPD = FALSE,
    components = 6,
    mcmc = list(niter = 200, nburnin = 50, nchains = 1)
  )

  nullfile <- if (.Platform$OS.type == "windows") "NUL" else "/dev/null"
  utils::capture.output(
    fit <- run_mcmc_bundle_manual(bundle, show_progress = FALSE),
    file = nullfile
  )

  pred_credible <- predict(fit, type = "mean", interval = "credible", nsim_mean = 100)
  pred_hpd <- predict(fit, type = "mean", interval = "hpd", nsim_mean = 100)

  expect_s3_class(pred_credible, "mixgpd_predict")
  expect_s3_class(pred_hpd, "mixgpd_predict")

  # Both should have valid intervals
  expect_true(pred_credible$fit$lower[1] <= pred_credible$fit$estimate[1])
  expect_true(pred_hpd$fit$lower[1] <= pred_hpd$fit$estimate[1])
})

test_that("bulk-only dpmix mean is analytical across kernels and cauchy rejects it", {
  init_kernel_registry()

  testthat::local_mocked_bindings(
    .validate_fit = function(object) TRUE,
    .cmgpd_progress_start = function(...) list(),
    .cmgpd_progress_step = function(...) invisible(NULL),
    .cmgpd_progress_done = function(...) invisible(NULL),
    .extract_draws_matrix = function(object) object$draw_mat,
    .package = "CausalMixGPD"
  )

  assemble_draw_mat <- function(weights, params = list(), betas = list(), scalar_draws = list()) {
    n_draws <- nrow(weights)
    cols <- list()
    for (k in seq_len(ncol(weights))) cols[[paste0("w[", k, "]")]] <- weights[, k]
    for (nm in names(params)) {
      mat <- params[[nm]]
      for (k in seq_len(ncol(mat))) cols[[paste0(nm, "[", k, "]")]] <- mat[, k]
    }
    for (nm in names(scalar_draws)) cols[[nm]] <- as.numeric(scalar_draws[[nm]])
    for (nm in names(betas)) {
      arr <- betas[[nm]]
      for (k in seq_len(dim(arr)[2])) {
        for (j in seq_len(dim(arr)[3])) {
          cols[[paste0("beta_", nm, "[", k, ",", j, "]")]] <- arr[, k, j]
        }
      }
    }
    out <- matrix(NA_real_, nrow = n_draws, ncol = length(cols))
    colnames(out) <- names(cols)
    for (j in seq_along(cols)) out[, j] <- cols[[j]]
    out
  }

  make_fit <- function(kernel, draw_mat, X = NULL, link_params = list(), GPD = FALSE, backend = "sb", gpd = list()) {
    structure(
      list(
        spec = list(
          meta = list(backend = backend, kernel = kernel, GPD = GPD, has_X = !is.null(X)),
          dispatch = list(backend = backend, GPD = GPD, link_params = link_params, gpd = gpd)
        ),
        data = list(
          X = X,
          y = rep(1, if (is.null(X)) 1L else nrow(X))
        ),
        draw_mat = draw_mat,
        epsilon = 1e-8
      ),
      class = "mixgpd_fit"
    )
  }

  weighted_mean <- function(w, component_mean) {
    sum((w / sum(w)) * component_mean)
  }

  kernel_component_mean <- function(kernel, pars) {
    switch(kernel,
      normal = pars$mean,
      lognormal = exp(pars$meanlog + 0.5 * pars$sdlog^2),
      gamma = pars$shape * pars$scale,
      laplace = pars$location,
      invgauss = pars$mean,
      amoroso = pars$loc + pars$scale * exp(lgamma(pars$shape1 + 1 / pars$shape2) - lgamma(pars$shape1))
    )
  }

  kernel_component_trunc_mean <- function(kernel, pars, threshold) {
    switch(kernel,
      normal = {
        a <- (threshold - pars$mean) / pars$sd
        pars$mean * pnorm(a) - pars$sd * dnorm(a)
      },
      lognormal = {
        if (threshold <= 0) {
          rep(0, length(pars$meanlog))
        } else {
          exp(pars$meanlog + 0.5 * pars$sdlog^2) *
            pnorm((log(threshold) - pars$meanlog - pars$sdlog^2) / pars$sdlog)
        }
      },
      gamma = {
        if (threshold <= 0) {
          rep(0, length(pars$shape))
        } else {
          pars$shape * pars$scale * pgamma(threshold, shape = pars$shape + 1, scale = pars$scale)
        }
      },
      laplace = ifelse(
        threshold < pars$location,
        0.5 * (threshold - pars$scale) * exp((threshold - pars$location) / pars$scale),
        pars$location - 0.5 * (threshold + pars$scale) * exp(-(threshold - pars$location) / pars$scale)
      ),
      invgauss = {
        if (threshold <= 0) {
          rep(0, length(pars$mean))
        } else {
          sqrt_term <- sqrt(pars$shape / threshold)
          a <- sqrt_term * (threshold / pars$mean - 1)
          b <- -sqrt_term * (threshold / pars$mean + 1)
          pars$mean * (pnorm(a) - exp(2 * pars$shape / pars$mean) * pnorm(b))
        }
      },
      amoroso = {
        z <- pmax((threshold - pars$loc) / pars$scale, 0)^pars$shape2
        pars$loc * pgamma(z, shape = pars$shape1, scale = 1) +
          pars$scale * exp(lgamma(pars$shape1 + 1 / pars$shape2) - lgamma(pars$shape1)) *
          pgamma(z, shape = pars$shape1 + 1 / pars$shape2, scale = 1)
      }
    )
  }

  kernel_spliced_mean <- function(kernel, w, pars, threshold, tail_scale, tail_shape) {
    trunc_mean <- kernel_component_trunc_mean(kernel, pars, threshold)
    Fu <- switch(kernel,
      normal = pnormmix(threshold, w = w, mean = pars$mean, sd = pars$sd),
      lognormal = plognormalmix(threshold, w = w, meanlog = pars$meanlog, sdlog = pars$sdlog),
      gamma = pgammamix(threshold, w = w, shape = pars$shape, scale = pars$scale),
      laplace = plaplacemix(threshold, w = w, location = pars$location, scale = pars$scale),
      invgauss = pinvgaussmix(threshold, w = w, mean = pars$mean, shape = pars$shape),
      amoroso = pamorosomix(threshold, w = w, loc = pars$loc, scale = pars$scale, shape1 = pars$shape1, shape2 = pars$shape2)
    )
    weighted_mean(w, trunc_mean) + (1 - Fu) * (threshold + tail_scale / (1 - tail_shape))
  }

  unconditional_cases <- list(
    normal = list(
      weights = rbind(c(0.7, 0.3), c(0.2, 0.8)),
      params = list(mean = rbind(c(1, 3), c(2, 4)), sd = rbind(c(1.1, 1.4), c(1.2, 1.5)))
    ),
    lognormal = list(
      weights = rbind(c(0.6, 0.4), c(0.3, 0.7)),
      params = list(meanlog = rbind(c(0.1, 0.4), c(0.2, 0.5)), sdlog = rbind(c(0.3, 0.2), c(0.25, 0.35)))
    ),
    gamma = list(
      weights = rbind(c(0.5, 0.5), c(0.25, 0.75)),
      params = list(shape = rbind(c(2, 4), c(3, 5)), scale = rbind(c(1.0, 0.8), c(0.9, 1.1)))
    ),
    laplace = list(
      weights = rbind(c(0.65, 0.35), c(0.45, 0.55)),
      params = list(location = rbind(c(-1, 2), c(0.5, 1.5)), scale = rbind(c(1.2, 0.7), c(1.1, 0.9)))
    ),
    invgauss = list(
      weights = rbind(c(0.55, 0.45), c(0.4, 0.6)),
      params = list(mean = rbind(c(1.5, 2.5), c(2.0, 3.5)), shape = rbind(c(4, 6), c(5, 7)))
    ),
    amoroso = list(
      weights = rbind(c(0.6, 0.4), c(0.35, 0.65)),
      params = list(
        loc = rbind(c(0.2, 0.5), c(0.3, 0.6)),
        scale = rbind(c(1.1, 0.9), c(1.2, 1.0)),
        shape1 = rbind(c(2.5, 3.0), c(2.2, 2.8)),
        shape2 = rbind(c(1.4, 1.6), c(1.5, 1.8))
      )
    )
  )

  for (kernel in names(unconditional_cases)) {
    case <- unconditional_cases[[kernel]]
    draw_mat <- assemble_draw_mat(case$weights, params = case$params)
    fit <- make_fit(kernel, draw_mat = draw_mat)
    pred <- .predict_mixgpd(fit, type = "mean", store_draws = TRUE, show_progress = FALSE)
    expected_draws <- vapply(seq_len(nrow(case$weights)), function(s) {
      pars <- lapply(case$params, function(mat) mat[s, ])
      weighted_mean(case$weights[s, ], kernel_component_mean(kernel, pars))
    }, numeric(1))
    expect_equal(as.numeric(pred$draws), expected_draws, tolerance = 1e-10, info = kernel)
    expect_equal(pred$fit$estimate, mean(expected_draws), tolerance = 1e-10, info = kernel)
    expect_equal(pred$diagnostics$mean_method, "analytic", info = kernel)
    expect_null(pred$diagnostics$nsim_mean, info = kernel)
  }

  Xpred <- matrix(c(1, 2), ncol = 1, dimnames = list(NULL, "x1"))
  conditional_cases <- list(
    normal = list(
      weights = rbind(c(0.7, 0.3), c(0.2, 0.8)),
      params = list(sd = rbind(c(1.1, 1.4), c(1.2, 1.5))),
      betas = list(mean = array(c(1, 2, 0, 3), dim = c(2, 2, 1))),
      link_params = list(mean = list(mode = "link", link = "identity"))
    ),
    lognormal = list(
      weights = rbind(c(0.6, 0.4), c(0.3, 0.7)),
      params = list(sdlog = rbind(c(0.3, 0.2), c(0.25, 0.35))),
      betas = list(meanlog = array(c(0.1, 0.3, 0.2, 0.4), dim = c(2, 2, 1))),
      link_params = list(meanlog = list(mode = "link", link = "identity"))
    ),
    gamma = list(
      weights = rbind(c(0.5, 0.5), c(0.25, 0.75)),
      params = list(shape = rbind(c(2, 4), c(3, 5))),
      betas = list(scale = array(c(log(1.0), log(0.6), log(0.8), log(0.5)), dim = c(2, 2, 1))),
      link_params = list(scale = list(mode = "link", link = "exp"))
    ),
    laplace = list(
      weights = rbind(c(0.65, 0.35), c(0.45, 0.55)),
      params = list(scale = rbind(c(1.2, 0.7), c(1.1, 0.9))),
      betas = list(location = array(c(-0.5, 0.8, 0.25, 1.0), dim = c(2, 2, 1))),
      link_params = list(location = list(mode = "link", link = "identity"))
    ),
    invgauss = list(
      weights = rbind(c(0.55, 0.45), c(0.4, 0.6)),
      params = list(shape = rbind(c(4, 6), c(5, 7))),
      betas = list(mean = array(c(log(1.2), log(2.0), log(1.5), log(2.2)), dim = c(2, 2, 1))),
      link_params = list(mean = list(mode = "link", link = "exp"))
    ),
    amoroso = list(
      weights = rbind(c(0.6, 0.4), c(0.35, 0.65)),
      params = list(
        scale = rbind(c(1.1, 0.9), c(1.2, 1.0)),
        shape1 = rbind(c(2.5, 3.0), c(2.2, 2.8)),
        shape2 = rbind(c(1.4, 1.6), c(1.5, 1.8))
      ),
      betas = list(loc = array(c(0.2, 0.4, 0.3, 0.5), dim = c(2, 2, 1))),
      link_params = list(loc = list(mode = "link", link = "identity"))
    )
  )

  apply_link <- function(x, link) {
    if (identical(link, "exp")) return(exp(x))
    x
  }

  for (kernel in names(conditional_cases)) {
    case <- conditional_cases[[kernel]]
    draw_mat <- assemble_draw_mat(case$weights, params = case$params, betas = case$betas)
    fit <- make_fit(kernel, draw_mat = draw_mat, X = Xpred, link_params = case$link_params)
    pred <- .predict_mixgpd(fit, x = Xpred, type = "mean", store_draws = TRUE, show_progress = FALSE)
    expected_draws <- matrix(NA_real_, nrow = nrow(case$weights), ncol = nrow(Xpred))
    link_name <- names(case$link_params)[1]
    link_fun <- case$link_params[[1]]$link
    beta_arr <- case$betas[[link_name]]
    for (s in seq_len(nrow(case$weights))) {
      for (i in seq_len(nrow(Xpred))) {
        pars <- lapply(case$params, function(mat) mat[s, ])
        pars[[link_name]] <- apply_link(Xpred[i, 1] * beta_arr[s, , 1], link_fun)
        expected_draws[s, i] <- weighted_mean(case$weights[s, ], kernel_component_mean(kernel, pars))
      }
    }
    expect_equal(pred$draws, expected_draws, tolerance = 1e-10, info = kernel)
    expect_equal(pred$fit$estimate, colMeans(expected_draws), tolerance = 1e-10, info = kernel)
    expect_equal(pred$diagnostics$mean_method, "analytic", info = kernel)
  }

  gpd_cases <- list(
    normal = list(
      weights = rbind(c(0.7, 0.3), c(0.2, 0.8)),
      params = list(mean = rbind(c(1, 3), c(2, 4)), sd = rbind(c(1.1, 1.4), c(1.2, 1.5))),
      threshold = c(2.5, 3.0),
      tail_scale = c(0.8, 0.6),
      tail_shape = c(0.2, 0.4)
    ),
    lognormal = list(
      weights = rbind(c(0.6, 0.4), c(0.3, 0.7)),
      params = list(meanlog = rbind(c(0.1, 0.4), c(0.2, 0.5)), sdlog = rbind(c(0.3, 0.2), c(0.25, 0.35))),
      threshold = c(1.8, 2.4),
      tail_scale = c(0.7, 0.5),
      tail_shape = c(0.25, 0.3)
    ),
    gamma = list(
      weights = rbind(c(0.5, 0.5), c(0.25, 0.75)),
      params = list(shape = rbind(c(2, 4), c(3, 5)), scale = rbind(c(1.0, 0.8), c(0.9, 1.1))),
      threshold = c(3.0, 4.2),
      tail_scale = c(0.9, 0.7),
      tail_shape = c(0.15, 0.35)
    ),
    laplace = list(
      weights = rbind(c(0.65, 0.35), c(0.45, 0.55)),
      params = list(location = rbind(c(-1, 2), c(0.5, 1.5)), scale = rbind(c(1.2, 0.7), c(1.1, 0.9))),
      threshold = c(1.2, 2.1),
      tail_scale = c(0.75, 0.55),
      tail_shape = c(0.1, 0.2)
    ),
    invgauss = list(
      weights = rbind(c(0.55, 0.45), c(0.4, 0.6)),
      params = list(mean = rbind(c(1.5, 2.5), c(2.0, 3.5)), shape = rbind(c(4, 6), c(5, 7))),
      threshold = c(2.5, 3.2),
      tail_scale = c(0.65, 0.8),
      tail_shape = c(0.2, 0.45)
    ),
    amoroso = list(
      weights = rbind(c(0.6, 0.4), c(0.35, 0.65)),
      params = list(
        loc = rbind(c(0.2, 0.5), c(0.3, 0.6)),
        scale = rbind(c(1.1, 0.9), c(1.2, 1.0)),
        shape1 = rbind(c(2.5, 3.0), c(2.2, 2.8)),
        shape2 = rbind(c(1.4, 1.6), c(1.5, 1.8))
      ),
      threshold = c(2.1, 2.4),
      tail_scale = c(0.7, 0.6),
      tail_shape = c(0.2, 0.3)
    )
  )

  for (kernel in names(gpd_cases)) {
    case <- gpd_cases[[kernel]]
    draw_mat <- assemble_draw_mat(
      case$weights,
      params = case$params,
      scalar_draws = list(
        threshold = case$threshold,
        tail_scale = case$tail_scale,
        tail_shape = case$tail_shape
      )
    )
    fit <- make_fit(
      kernel,
      draw_mat = draw_mat,
      GPD = TRUE,
      gpd = list(
        threshold = list(mode = "constant"),
        tail_scale = list(mode = "constant"),
        tail_shape = list(mode = "constant")
      )
    )
    pred <- .predict_mixgpd(fit, type = "mean", store_draws = TRUE, show_progress = FALSE)
    expected_draws <- vapply(seq_len(nrow(case$weights)), function(s) {
      pars <- lapply(case$params, function(mat) mat[s, ])
      kernel_spliced_mean(
        kernel,
        w = case$weights[s, ],
        pars = pars,
        threshold = case$threshold[s],
        tail_scale = case$tail_scale[s],
        tail_shape = case$tail_shape[s]
      )
    }, numeric(1))
    expect_equal(as.numeric(pred$draws), expected_draws, tolerance = 1e-8, info = kernel)
    expect_equal(pred$fit$estimate, mean(expected_draws), tolerance = 1e-8, info = kernel)
    expect_equal(pred$diagnostics$mean_method, "analytic", info = kernel)
    expect_null(pred$diagnostics$nsim_mean, info = kernel)
  }

  gpd_cond_draw_mat <- assemble_draw_mat(
    weights = rbind(c(0.7, 0.3), c(0.2, 0.8)),
    params = list(sd = rbind(c(1.1, 1.4), c(1.2, 1.5))),
    betas = list(mean = array(c(1, 2, 0, 3), dim = c(2, 2, 1))),
    scalar_draws = list(
      threshold = c(2.5, 3.0),
      tail_scale = c(0.8, 0.6),
      tail_shape = c(0.2, 0.4)
    )
  )
  gpd_cond_fit <- make_fit(
    "normal",
    draw_mat = gpd_cond_draw_mat,
    X = Xpred,
    link_params = list(mean = list(mode = "link", link = "identity")),
    GPD = TRUE,
    gpd = list(
      threshold = list(mode = "constant"),
      tail_scale = list(mode = "constant"),
      tail_shape = list(mode = "constant")
    )
  )
  gpd_cond_pred <- .predict_mixgpd(gpd_cond_fit, x = Xpred, type = "mean", store_draws = TRUE, show_progress = FALSE)
  gpd_cond_expected <- matrix(NA_real_, nrow = 2, ncol = nrow(Xpred))
  beta_arr <- array(c(1, 2, 0, 3), dim = c(2, 2, 1))
  for (s in 1:2) {
    for (i in seq_len(nrow(Xpred))) {
      pars <- list(
        mean = Xpred[i, 1] * beta_arr[s, , 1],
        sd = c(1.1, 1.4)
      )
      if (s == 2) pars$sd <- c(1.2, 1.5)
      gpd_cond_expected[s, i] <- kernel_spliced_mean(
        "normal",
        w = if (s == 1) c(0.7, 0.3) else c(0.2, 0.8),
        pars = pars,
        threshold = if (s == 1) 2.5 else 3.0,
        tail_scale = if (s == 1) 0.8 else 0.6,
        tail_shape = if (s == 1) 0.2 else 0.4
      )
    }
  }
  expect_equal(gpd_cond_pred$draws, gpd_cond_expected, tolerance = 1e-8)
  expect_equal(gpd_cond_pred$fit$estimate, colMeans(gpd_cond_expected), tolerance = 1e-8)
  expect_equal(gpd_cond_pred$diagnostics$mean_method, "analytic")

  xi_bad_draw_mat <- assemble_draw_mat(
    weights = rbind(c(0.7, 0.3), c(0.2, 0.8)),
    params = list(mean = rbind(c(1, 3), c(2, 4)), sd = rbind(c(1.1, 1.4), c(1.2, 1.5))),
    scalar_draws = list(threshold = c(2.5, 3.0), tail_scale = c(0.8, 0.6), tail_shape = c(1.1, 0.4))
  )
  xi_bad_fit <- make_fit(
    "normal",
    draw_mat = xi_bad_draw_mat,
    GPD = TRUE,
    gpd = list(
      threshold = list(mode = "constant"),
      tail_scale = list(mode = "constant"),
      tail_shape = list(mode = "constant")
    )
  )
  expect_error(
    .predict_mixgpd(xi_bad_fit, type = "mean", show_progress = FALSE),
    "tail_shape \\(xi\\) >= 1; use type='rmean'"
  )

  cauchy_draw_mat <- assemble_draw_mat(
    weights = rbind(c(0.6, 0.4), c(0.5, 0.5)),
    params = list(location = rbind(c(0, 1), c(-1, 2)), scale = rbind(c(1, 2), c(1.5, 0.75)))
  )
  cauchy_fit <- make_fit("cauchy", draw_mat = cauchy_draw_mat)
  expect_error(
    .predict_mixgpd(cauchy_fit, type = "mean", show_progress = FALSE),
    "Mean is not supported for the Cauchy kernel; use type='rmean'"
  )
})


# ==============================================================================
# Test fitted.mixgpd_fit with interval parameter
# ==============================================================================

test_that("fitted.mixgpd_fit supports unified interval parameter", {
  skip_if_not_test_level("ci")

  set.seed(42)
  n <- 30
  y <- abs(rnorm(n)) + 0.2
  X <- data.frame(x1 = rnorm(n), x2 = runif(n))

  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "sb",
    kernel = "gamma",
    GPD = FALSE,
    components = 6,
    mcmc = list(niter = 200, nburnin = 50, nchains = 1)
  )

  nullfile <- if (.Platform$OS.type == "windows") "NUL" else "/dev/null"
  utils::capture.output(
    fit <- run_mcmc_bundle_manual(bundle, show_progress = FALSE),
    file = nullfile
  )

  fitted_credible <- fitted(fit, type = "mean", interval = "credible")
  fitted_hpd <- fitted(fit, type = "mean", interval = "hpd")
  fitted_none <- fitted(fit, type = "mean", interval = NULL)

  expect_s3_class(fitted_credible, "mixgpd_fitted")
  expect_s3_class(fitted_hpd, "mixgpd_fitted")
  expect_s3_class(fitted_none, "mixgpd_fitted")

  # Both should have valid structure
  expect_true(all(c("fit", "lower", "upper", "residuals") %in% names(fitted_credible)))
  expect_true(all(c("fit", "lower", "upper", "residuals") %in% names(fitted_hpd)))
  expect_true(all(c("fit", "lower", "upper", "residuals") %in% names(fitted_none)))

  # Check interval is stored
  expect_equal(attr(fitted_credible, "interval"), "credible")
  expect_equal(attr(fitted_hpd, "interval"), "hpd")
  expect_null(attr(fitted_none, "interval"))
})


# ==============================================================================
# Test .posterior_summarize with interval parameter
# ==============================================================================

test_that(".posterior_summarize uses interval correctly", {
  set.seed(42)
  draws <- matrix(rexp(500), nrow = 5, ncol = 100)

  summ_credible <- .posterior_summarize(draws, interval = "credible")
  summ_hpd <- .posterior_summarize(draws, interval = "hpd")
  summ_none <- .posterior_summarize(draws, interval = NULL)

  # Check structure
  expect_named(summ_credible, c("estimate", "lower", "upper", "q"))
  expect_named(summ_hpd, c("estimate", "lower", "upper", "q"))
  expect_named(summ_none, c("estimate", "lower", "upper", "q"))

  # Estimates should be the same (all use mean)
  expect_equal(summ_credible$estimate, summ_hpd$estimate)
  expect_equal(summ_credible$estimate, summ_none$estimate)

  # Intervals should be different for skewed distribution (credible vs hpd)
  expect_false(all(summ_credible$lower == summ_hpd$lower))
  expect_false(all(summ_credible$upper == summ_hpd$upper))

  # NULL interval should give NA for lower/upper
  expect_true(all(is.na(summ_none$lower)))
  expect_true(all(is.na(summ_none$upper)))
})

# ===== END unit/test-hpd-intervals.R =====

# ===== BEGIN unit/test-internal-utils.R =====
# tests/testthat/test-internal-utils.R
# Unit tests for internal utility helpers (02-utilities-internal.R)

# ======================================================================
# .validate_nimble_reserved_names
# ======================================================================

test_that(".validate_nimble_reserved_names accepts valid names", {
  expect_invisible(.validate_nimble_reserved_names(c("alpha", "beta", "sigma")))
  expect_invisible(.validate_nimble_reserved_names(character(0)))
  expect_invisible(.validate_nimble_reserved_names(NULL))
  expect_invisible(.validate_nimble_reserved_names(c("x_if", "myfor")))
})

test_that(".validate_nimble_reserved_names rejects reserved keywords", {
  expect_error(
    .validate_nimble_reserved_names(c("if", "alpha"), context = "test"),
    regexp = "reserved NIMBLE keywords"
  )
  expect_error(
    .validate_nimble_reserved_names(c("FOR", "WHILE")),
    regexp = "reserved NIMBLE keywords"
  )
  expect_error(
    .validate_nimble_reserved_names(c("T", "F")),
    regexp = "reserved NIMBLE keywords"
  )
  expect_error(
    .validate_nimble_reserved_names(c("na", "nan", "inf")),
    regexp = "reserved NIMBLE keywords"
  )
})

test_that(".validate_nimble_reserved_names handles edge cases", {
  expect_invisible(.validate_nimble_reserved_names(c("", NA_character_)))
  expect_invisible(.validate_nimble_reserved_names(c("iffy", "truly", "nullify")))
})

# ======================================================================
# Plot helpers / code wrappers
# ======================================================================

test_that(".plot_palette returns requested length", {
  expect_length(.plot_palette(3), 3)
  expect_length(.plot_palette(8), 8)
  expect_length(.plot_palette(12), 12)
})

test_that(".extract_nimble_code and .wrap_nimble_code handle wrapped/list code", {
  code_obj <- quote({ a <- 1 })
  wrapped <- .wrap_nimble_code(code_obj)
  expect_true(is.list(wrapped))
  expect_true(!is.null(wrapped$nimble))

  # extract returns $nimble when stored in a list wrapper
  extracted <- .extract_nimble_code(wrapped)
  expect_equal(deparse(extracted), deparse(code_obj))

  # if already a list, wrapper should pass-through
  pass <- .wrap_nimble_code(list(nimble = code_obj))
  expect_true(is.list(pass))
  expect_true(!is.null(pass$nimble))
})

test_that(".wrap_plotly returns plotly when available, otherwise passthrough", {
  p <- ggplot2::ggplot() + ggplot2::geom_point(ggplot2::aes(1, 1))

  out <- .wrap_plotly(p)
  if (requireNamespace("plotly", quietly = TRUE)) {
    old_opt <- getOption("CausalMixGPD.plotly")
    options(CausalMixGPD.plotly = TRUE)
    on.exit(options(CausalMixGPD.plotly = old_opt), add = TRUE)
    out <- .wrap_plotly(p)
    expect_true(is.list(out))
    expect_true(inherits(out, "plotly") || inherits(out, "htmlwidget"))
  } else {
    expect_true(inherits(out, "ggplot"))
  }

  plots <- list(a = p, b = p)
  class(plots) <- c("mixgpd_fit_plots", "list")
  out2 <- .wrap_plotly(plots)
  expect_true(inherits(out2, "mixgpd_fit_plots"))
  expect_true(is.list(out2))
})

# ======================================================================
# .coerce_fit_df
# ======================================================================

test_that(".coerce_fit_df handles data.frame input", {
  df <- data.frame(estimate = 1:3, lower = 0:2, upper = 2:4)
  result <- .coerce_fit_df(df)
  expect_s3_class(result, "data.frame")
  expect_true(all(c("estimate", "lower", "upper", "id") %in% names(result)))
  expect_equal(result$estimate, 1:3)
})

test_that(".coerce_fit_df handles data.frame missing estimate", {
  df <- data.frame(fit = c(1, 2), other = c(3, 4))
  result <- .coerce_fit_df(df)
  expect_equal(result$estimate, c(1, 2))
  expect_true("id" %in% names(result))
})

test_that(".coerce_fit_df handles data.frame missing fit/estimate column", {
  df <- data.frame(a = 1:3, b = 4:6)
  result <- .coerce_fit_df(df)
  expect_equal(result$estimate, 1:3)  # uses first numeric column
})

test_that(".coerce_fit_df handles numeric matrix with named columns", {
  mat <- matrix(1:9, nrow = 3, ncol = 3)
  colnames(mat) <- c("estimate", "lower", "upper")
  result <- .coerce_fit_df(mat)
  expect_s3_class(result, "data.frame")
  expect_true("id" %in% names(result))
  expect_equal(nrow(result), 3)
})

test_that(".coerce_fit_df handles matrix with probs expansion", {
  mat <- matrix(1:6, nrow = 2, ncol = 3)
  probs <- c(0.25, 0.5, 0.75)
  result <- .coerce_fit_df(mat, n_pred = 2, probs = probs)
  expect_s3_class(result, "data.frame")
  expect_true(all(c("id", "index", "estimate") %in% names(result)))
  expect_equal(nrow(result), 6)  # 2 rows x 3 probs
})

test_that(".coerce_fit_df handles matrix default (first column)", {

  mat <- matrix(1:6, nrow = 3, ncol = 2)
  result <- .coerce_fit_df(mat)
  expect_equal(result$estimate, 1:3)
  expect_equal(result$lower, 4:6)
})

test_that(".coerce_fit_df handles numeric vector", {
  vec <- c(1.0, 2.0, 3.0)
  result <- .coerce_fit_df(vec)
  expect_s3_class(result, "data.frame")
  expect_equal(result$estimate, vec)
  expect_equal(result$id, 1:3)
})

test_that(".coerce_fit_df errors on unsupported type", {
  expect_error(.coerce_fit_df(list(a = 1)), regexp = "Cannot coerce")
})

# ======================================================================
# .compute_interval
# ======================================================================

test_that(".compute_interval credible returns expected quantiles", {
  draws <- rnorm(1000, mean = 5, sd = 1)
  iv <- .compute_interval(draws, level = 0.95, type = "credible")
  expect_named(iv, c("lower", "upper"))
  expect_true(iv["lower"] < 5)

  expect_true(iv["upper"] > 5)
  expect_true(iv["lower"] < iv["upper"])
})

test_that(".compute_interval hpd returns valid interval", {
  skip_if_not_installed("coda")
  draws <- rnorm(1000, mean = 0, sd = 1)
  iv <- .compute_interval(draws, level = 0.95, type = "hpd")
  expect_named(iv, c("lower", "upper"))
  expect_true(iv["lower"] < iv["upper"])
})
test_that(".compute_interval handles insufficient draws", {
  expect_equal(.compute_interval(numeric(0)), c(lower = NA_real_, upper = NA_real_))
  expect_equal(.compute_interval(c(1)), c(lower = NA_real_, upper = NA_real_))
})

test_that(".compute_interval filters non-finite", {
  draws <- c(1, 2, 3, NA, Inf, -Inf, 4, 5)
  iv <- .compute_interval(draws, level = 0.95, type = "credible")
  expect_true(is.finite(iv["lower"]))
  expect_true(is.finite(iv["upper"]))
})

# ======================================================================
# .posterior_summarize
# ======================================================================

test_that(".posterior_summarize works on numeric vector", {
  draws <- rnorm(100, mean = 10, sd = 2)
  result <- .posterior_summarize(draws, probs = c(0.025, 0.5, 0.975), interval = "credible")
  expect_named(result, c("estimate", "lower", "upper", "q"))
  expect_length(result$estimate, 1)
  expect_true(result$lower < result$upper)
})

test_that(".posterior_summarize works on 2D matrix", {
  draws <- matrix(rnorm(300), nrow = 3, ncol = 100)
  result <- .posterior_summarize(draws, probs = c(0.025, 0.975), interval = "credible")
  expect_length(result$estimate, 3)
  expect_length(result$lower, 3)
  expect_length(result$upper, 3)
})

test_that(".posterior_summarize works on 3D array", {
  draws <- array(rnorm(600), dim = c(2, 3, 100))
  result <- .posterior_summarize(draws, probs = c(0.025, 0.975), interval = "credible")
  expect_equal(dim(result$estimate), c(2, 3))
  expect_equal(dim(result$lower), c(2, 3))
  expect_equal(dim(result$upper), c(2, 3))
})

test_that(".posterior_summarize handles interval = NULL", {
  draws <- rnorm(100)
  result <- .posterior_summarize(draws, probs = c(0.025, 0.975), interval = NULL)
  expect_true(is.na(result$lower))
  expect_true(is.na(result$upper))
})

test_that(".posterior_summarize hpd mode", {
  skip_if_not_installed("coda")
  draws <- rnorm(200)
  result <- .posterior_summarize(draws, probs = c(0.025, 0.975), interval = "hpd")
  expect_true(!is.na(result$lower))
  expect_true(!is.na(result$upper))
})

# ======================================================================
# .truncate_components_one_draw
# ======================================================================

test_that(".truncate_components_one_draw sorts by weight", {
  w <- c(0.1, 0.6, 0.3)
  params <- list(mu = c(1, 2, 3), sigma = c(0.5, 1, 1.5))
  result <- .truncate_components_one_draw(w, params, epsilon = 0.01)

  expect_equal(result$ord, c(2, 3, 1))  # sorted descending by weight
  expect_true(sum(result$weights) > 0.99)
})

test_that(".truncate_components_one_draw truncates small components", {
  w <- c(0.05, 0.9, 0.05)
  params <- list(mu = c(1, 2, 3))
  result <- .truncate_components_one_draw(w, params, epsilon = 0.1)

  expect_true(result$k <= 3)
  expect_true(length(result$weights) <= 3)
})

test_that(".truncate_components_one_draw errors on invalid epsilon", {
  w <- c(0.5, 0.5)
  params <- list(mu = c(1, 2))

  expect_error(.truncate_components_one_draw(w, params, epsilon = -0.1))
  expect_error(.truncate_components_one_draw(w, params, epsilon = 1.0))
  expect_error(.truncate_components_one_draw(w, params, epsilon = NA))
  expect_error(.truncate_components_one_draw(w, params, epsilon = c(0.1, 0.2)))
})

test_that(".truncate_components_one_draw handles empty params", {
  w <- c(0.3, 0.7)
  result <- .truncate_components_one_draw(w, params = list(), epsilon = 0.01)
  expect_true(result$k >= 1)
})

test_that(".truncate_components_one_draw validates params length", {
  w <- c(0.5, 0.5)
  params <- list(mu = c(1, 2, 3))  # wrong length

  expect_error(.truncate_components_one_draw(w, params, epsilon = 0.1))
})

test_that(".truncate_components_one_draw adjusts weights to sum to 1", {
  w <- c(0.6, 0.3, 0.1)
  params <- list(mu = c(1, 2, 3))
  result <- .truncate_components_one_draw(w, params, epsilon = 0.05)

  expect_equal(sum(result$weights), 1, tolerance = 1e-10)
})

# ======================================================================
# Additional internal utility helpers
# ======================================================================

# ======================================================================
# .extract_nimble_code tests
# ======================================================================

test_that(".extract_nimble_code extracts from list with nimble", {
  code_obj <- list(nimble = quote(y ~ dnorm(0, 1)))
  result <- .extract_nimble_code(code_obj)
  expect_equal(result, quote(y ~ dnorm(0, 1)))
})

test_that(".extract_nimble_code extracts from list with code", {
  code_obj <- list(code = quote(y ~ dnorm(0, 1)))
  result <- .extract_nimble_code(code_obj)
  expect_equal(result, quote(y ~ dnorm(0, 1)))
})

test_that(".extract_nimble_code returns nimbleCode as-is", {
  code <- quote(y ~ dnorm(0, 1))
  class(code) <- c("nimbleCode", class(code))
  result <- .extract_nimble_code(code)
  expect_equal(result, code)
})

test_that(".extract_nimble_code returns non-list as-is", {
  result <- .extract_nimble_code("some_code")
  expect_equal(result, "some_code")
})

# ======================================================================
# .wrap_nimble_code tests
# ======================================================================

test_that(".wrap_nimble_code wraps non-list code", {
  code <- quote(y ~ dnorm(0, 1))
  result <- .wrap_nimble_code(code)
  expect_true(is.list(result))
  expect_equal(result$nimble, code)
})

test_that(".wrap_nimble_code returns list as-is", {
  code_list <- list(nimble = quote(y ~ dnorm(0, 1)))
  result <- .wrap_nimble_code(code_list)
  expect_equal(result, code_list)
})

# ======================================================================
# .backend_label tests
# ======================================================================

test_that(".backend_label returns correct labels", {
  expect_equal(.backend_label("sb"), "Stick-Breaking Process")
  expect_equal(.backend_label("crp"), "Chinese Restaurant Process")
  expect_equal(.backend_label("unknown"), "unknown")
})

# ======================================================================
# .kernel_label tests
# ======================================================================

test_that(".kernel_label returns correct labels", {
  expect_equal(.kernel_label("normal"), "Normal Distribution")
  expect_equal(.kernel_label("gamma"), "Gamma Distribution")
  expect_equal(.kernel_label("lognormal"), "Lognormal Distribution")
  expect_equal(.kernel_label("laplace"), "Laplace Distribution")
  expect_equal(.kernel_label("invgauss"), "Inverse Gaussian Distribution")
  expect_equal(.kernel_label("amoroso"), "Amoroso Distribution")
  expect_equal(.kernel_label("cauchy"), "Cauchy Distribution")
  expect_equal(.kernel_label("unknown"), "unknown")
})

# ======================================================================
# .plot_palette tests
# ======================================================================

test_that(".plot_palette returns valid hex colors", {
  pal <- .plot_palette(4)
  expect_length(pal, 4)
  expect_true(all(grepl("^#[0-9A-Fa-f]{6}$", pal)))
})

test_that(".plot_palette recycles for n > 8", {
  pal <- .plot_palette(10)
  expect_length(pal, 10)
})

test_that(".plot_palette handles NULL", {
  pal <- .plot_palette(NULL)
  expect_length(pal, 8)
})

# ======================================================================
# .get_epsilon tests
# ======================================================================

test_that(".get_epsilon returns provided epsilon", {
  obj <- list(spec = list(meta = list(epsilon = 0.05)))
  result <- .get_epsilon(obj, epsilon = 0.1)
  expect_equal(result, 0.1)
})

test_that(".get_epsilon extracts from spec$meta", {
  obj <- list(spec = list(meta = list(epsilon = 0.05)))
  result <- .get_epsilon(obj, epsilon = NULL)
  expect_equal(result, 0.05)
})

test_that(".get_epsilon extracts from object$epsilon", {
  obj <- list(epsilon = 0.03)
  result <- .get_epsilon(obj, epsilon = NULL)
  expect_equal(result, 0.03)
})

test_that(".get_epsilon returns default 0.025", {
  obj <- list()
  result <- .get_epsilon(obj, epsilon = NULL)
  expect_equal(result, 0.025)
})

# ======================================================================
# .get_nobs tests
# ======================================================================

test_that(".get_nobs extracts from data$y", {
  obj <- list(data = list(y = 1:100))
  result <- .get_nobs(obj)
  expect_equal(result, 100)
})

test_that(".get_nobs extracts from y directly", {
  obj <- list(y = 1:50)
  result <- .get_nobs(obj)
  expect_equal(result, 50)
})

test_that(".get_nobs returns NA when not found", {
  obj <- list()
  result <- .get_nobs(obj)
  expect_true(is.na(result))
})

# ===== END unit/test-internal-utils.R =====

# ===== BEGIN unit/test-kernels.R =====
# test-kernels.R
# Consolidated kernel tests: registry, support table, distribution roundtrips (Tier A/B)
# Merged from: test-kernel-registry.R, test-kernel-support.R, test-kernels.R

# =============================================================================
# Kernel registry (from test-kernel-registry.R)
# =============================================================================
test_that("init_kernel_registry() returns TRUE", {
  result <- init_kernel_registry()
  expect_true(result)
  result2 <- init_kernel_registry()
  expect_true(result2)
})

test_that("get_kernel_registry() returns a complete registry", {
  registry <- get_kernel_registry()
  expect_true(is.list(registry))
  expected_kernels <- c("normal", "lognormal", "invgauss", "gamma", "laplace", "amoroso", "cauchy")
  expect_true(all(expected_kernels %in% names(registry)))
  expect_equal(length(registry), 7L)
})

test_that("Each kernel has required fields", {
  registry <- get_kernel_registry()
  for (kernel_name in names(registry)) {
    kernel <- registry[[kernel_name]]
    expect_true("key" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'key'"))
    expect_true("bulk_params" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'bulk_params'"))
    expect_true("bulk_support" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'bulk_support'"))
    expect_true("param_types" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'param_types'"))
    expect_true("allow_gpd" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'allow_gpd'"))
    expect_true("defaults_X" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'defaults_X'"))
    expect_true("sb" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'sb'"))
    expect_true("crp" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'crp'"))
    expect_true("signatures" %in% names(kernel), info = paste0("kernel '", kernel_name, "' missing 'signatures'"))
    expect_equal(kernel$key, kernel_name, info = paste0("kernel '", kernel_name, "' key mismatch"))
  }
})

test_that("Kernel bulk_params are character vectors", {
  registry <- get_kernel_registry()
  for (kernel_name in names(registry)) {
    kernel <- registry[[kernel_name]]
    expect_true(is.character(kernel$bulk_params), info = paste0("kernel '", kernel_name, "' bulk_params not character"))
    expect_true(length(kernel$bulk_params) >= 2L, info = paste0("kernel '", kernel_name, "' should have at least 2 bulk params"))
  }
})

test_that("Kernel SB signatures are properly structured", {
  registry <- get_kernel_registry()
  for (kernel_name in names(registry)) {
    kernel <- registry[[kernel_name]]
    sb <- kernel$sb
    expect_true(is.list(sb), info = paste0("kernel '", kernel_name, "' sb not a list"))
    expect_true("d" %in% names(sb), info = paste0("kernel '", kernel_name, "' sb missing 'd'"))
    expect_true("args" %in% names(sb), info = paste0("kernel '", kernel_name, "' sb missing 'args'"))
  }
})

test_that("Kernel CRP signatures are properly structured", {
  registry <- get_kernel_registry()
  for (kernel_name in names(registry)) {
    kernel <- registry[[kernel_name]]
    crp <- kernel$crp
    expect_true(is.list(crp), info = paste0("kernel '", kernel_name, "' crp not a list"))
    expect_true("d_base" %in% names(crp), info = paste0("kernel '", kernel_name, "' crp missing 'd_base'"))
  }
})

test_that("get_tail_registry() returns tail metadata", {
  tail_reg <- get_tail_registry()
  expect_true(is.list(tail_reg))
  expect_true("params" %in% names(tail_reg))
  expect_true("support" %in% names(tail_reg))
  expect_true("indexed_by_cluster_in_crp" %in% names(tail_reg))
  expect_equal(tail_reg$params, c("threshold", "tail_scale", "tail_shape"))
  expect_true(is.character(tail_reg$support))
  expect_equal(names(tail_reg$support), c("threshold", "tail_scale", "tail_shape"))
  expect_true(is.logical(tail_reg$indexed_by_cluster_in_crp))
})

test_that("Cauchy kernel does not allow GPD", {
  registry <- get_kernel_registry()
  expect_false(registry$cauchy$allow_gpd)
})

test_that("Most kernels allow GPD", {
  registry <- get_kernel_registry()
  gpd_kernels <- c("normal", "lognormal", "invgauss", "gamma", "laplace", "amoroso")
  for (k in gpd_kernels) {
    expect_true(registry[[k]]$allow_gpd, info = paste0("kernel '", k, "' should allow GPD"))
  }
})

# =============================================================================
# Kernel support table (from test-kernel-support.R)
# =============================================================================
test_that("kernel_support_table() returns a data frame", {
  result <- kernel_support_table()
  expect_true(is.data.frame(result))
})

test_that("kernel_support_table() has correct columns", {
  result <- kernel_support_table()
  expected_cols <- c("kernel", "gpd", "covariates", "sb", "crp")
  expect_equal(names(result), expected_cols)
})

test_that("kernel_support_table() includes all 7 kernels", {
  result <- kernel_support_table()
  expected_kernels <- c("normal", "lognormal", "invgauss", "gamma", "laplace", "amoroso", "cauchy")
  expect_equal(nrow(result), 7L)
  expect_true(all(expected_kernels %in% result$kernel))
})

test_that("kernel_support_table(round=TRUE) produces checkmark/cross symbols", {
  result <- kernel_support_table(round = TRUE)
  checkmark <- "\u2714"
  cross <- "\u274C"
  for (col in c("gpd", "covariates", "sb", "crp")) {
    expect_true(all(result[[col]] %in% c(checkmark, cross)), info = paste0("column '", col, "' should contain only checkmark/cross"))
  }
})

test_that("kernel_support_table(round=FALSE) produces logical values", {
  result <- kernel_support_table(round = FALSE)
  for (col in c("gpd", "covariates", "sb", "crp")) {
    expect_true(is.logical(result[[col]]), info = paste0("column '", col, "' should be logical"))
  }
})

test_that("Cauchy kernel does not support GPD in support table", {
  result <- kernel_support_table(round = FALSE)
  cauchy_row <- result[result$kernel == "cauchy", ]
  expect_false(cauchy_row$gpd)
})

test_that("All kernels support both sb and crp backends in support table", {
  result <- kernel_support_table(round = FALSE)
  expect_true(all(result$sb))
  expect_true(all(result$crp))
})

test_that("All kernels support covariates in support table", {
  result <- kernel_support_table(round = FALSE)
  expect_true(all(result$covariates))
})

test_that("Most kernels support GPD in support table (except Cauchy)", {
  result <- kernel_support_table(round = FALSE)
  gpd_kernels <- result[result$kernel != "cauchy", ]
  expect_true(all(gpd_kernels$gpd))
})

test_that("kernel_support_table() kernel column is character", {
  result <- kernel_support_table()
  expect_true(is.character(result$kernel))
})

# =============================================================================
# Kernel distribution roundtrips (from test-kernels.R)
# =============================================================================
test_that(
  "Kernels: mixture families (mix / mix+GPD / scalar+GPD) semantic + roundtrip checks",
  {
    # This test runs kernel math checks - fast, no MCMC compilation

    # ---------- helpers ----------
    .norm_w <- function(w) {
      w <- as.numeric(w)
      s <- sum(w)
      if (!is.finite(s) || s <= 0) rep(1 / length(w), length(w)) else w / s
    }

    .expect_density_ok <- function(dval, label = "") {
      expect_true(is.numeric(dval), info = label)
      expect_true(length(dval) >= 1, info = label)
      expect_true(all(is.finite(dval)), info = label)
      expect_true(all(dval >= 0), info = label)
    }

    .expect_cdf_ok <- function(pval, label = "") {
      expect_true(is.numeric(pval), info = label)
      expect_true(length(pval) >= 1, info = label)
      expect_true(all(is.finite(pval)), info = label)
      expect_true(all(pval >= 0 & pval <= 1), info = label)
    }

    .expect_rng_ok <- function(x, support = c(-Inf, Inf), label = "") {
      expect_true(is.numeric(x), info = label)
      expect_true(length(x) == 1, info = label)
      expect_true(is.finite(x), info = label)
      expect_true(x >= support[1] && x <= support[2], info = label)
    }

    # Roundtrip q(p) -> p(q) with diagnostics (realistic tolerances for numeric inversion).
    .roundtrip_qp <- function(qFun, pFun, probs, args, tol, label) {
      qv <- do.call(qFun, c(list(p = probs), args))
      expect_true(is.numeric(qv), info = paste0(label, " : qFun returned non-numeric"))
      expect_true(length(qv) == length(probs), info = paste0(label, " : qFun length mismatch"))

      okq <- is.finite(qv) & probs > 0 & probs < 1
      if (!any(okq)) {
        succeed(paste0(label, " : roundtrip skipped (no finite q)"))
        return(invisible(TRUE))
      }

      # IMPORTANT: for nimbleFunctions, pass integer flags (1/0), not TRUE/FALSE.
      pb <- vapply(
        qv[okq],
        function(qq) do.call(pFun, c(list(q = qq), args, list(lower.tail = 1L, log.p = 0L))),
        numeric(1)
      )

      if (any(!is.finite(pb))) {
        bad <- which(!is.finite(pb))[1]
        msg <- paste0(
          label, " : pFun returned non-finite at index ", bad,
          " (p=", probs[okq][bad], ", q=", qv[okq][bad], ")."
        )
        fail(msg)
      }

      .expect_cdf_ok(pb, label = paste0(label, " : pFun range"))

      err <- max(abs(pb - probs[okq]))
      expect_true(
        err < tol,
        info = paste0(label, " : max|p_back - p| = ", signif(err, 6), " (tol=", tol, ")")
      )
      invisible(TRUE)
    }

    .check_mix_family <- function(
    name,
    support = c(-Inf, Inf),
    dMix, pMix, rMix, qMix,
    args_mix,
    grid = NULL,
    dMixGpd = NULL, pMixGpd = NULL, rMixGpd = NULL, qMixGpd = NULL,
    args_mixgpd = NULL,
    dScalarGpd = NULL, pScalarGpd = NULL, rScalarGpd = NULL, qScalarGpd = NULL,
    args_scalargpd = NULL,
    tol_mix = 1e-4,
    tol_splice = 5e-4
    ) {
      # normalize weights for stability
      if ("w" %in% names(args_mix)) args_mix$w <- .norm_w(args_mix$w)
      if (!is.null(args_mixgpd) && "w" %in% names(args_mixgpd)) args_mixgpd$w <- .norm_w(args_mixgpd$w)

      if (is.null(grid)) {
        grid <- if (is.finite(support[1])) support[1] + c(0.01, 0.2, 1, 2) else c(-2, -0.5, 0, 0.5, 2)
        grid <- grid[grid >= support[1]]
        if (length(grid) == 0) grid <- support[1] + 0.5
      }

      # --- MIX checks ---
      dvals <- vapply(grid, function(xx) do.call(dMix, c(list(x = xx), args_mix, list(log = 0L))), numeric(1))
      pvals <- vapply(grid, function(xx) do.call(pMix, c(list(q = xx), args_mix, list(lower.tail = 1L, log.p = 0L))), numeric(1))

      .expect_density_ok(dvals, label = paste0(name, " : dMix"))
      .expect_cdf_ok(pvals, label = paste0(name, " : pMix"))

      r1 <- do.call(rMix, c(list(n = 1L), args_mix))
      .expect_rng_ok(r1, support = support, label = paste0(name, " : rMix"))

      probs <- c(0.05, 0.2, 0.5, 0.8, 0.95)
      .roundtrip_qp(qMix, pMix, probs, args_mix, tol = tol_mix, label = paste0(name, " : mix roundtrip"))

      succeed(paste0("mix checks ok: ", name))

      # --- MIX+GPD checks ---
      if (!is.null(dMixGpd) && !is.null(pMixGpd) && !is.null(rMixGpd) && !is.null(qMixGpd)) {
        expect_true(!is.null(args_mixgpd), info = paste0(name, " : args_mixgpd missing"))

        u <- args_mixgpd$threshold

        # continuity at threshold: F_mix(u) == F_mixgpd(u)
        Fu0 <- do.call(pMix,    c(list(q = u), args_mix,    list(lower.tail = 1L, log.p = 0L)))
        Fu1 <- do.call(pMixGpd, c(list(q = u), args_mixgpd, list(lower.tail = 1L, log.p = 0L)))

        .expect_cdf_ok(Fu0, label = paste0(name, " : Fu0"))
        .expect_cdf_ok(Fu1, label = paste0(name, " : Fu1"))

        expect_equal(
          as.numeric(Fu1), as.numeric(Fu0),
          tolerance = 1e-6,
          info = paste0(name, " : continuity at threshold")
        )

        grid2 <- sort(unique(c(u - 0.1, u, u + 0.1, u + 1)))
        d2 <- vapply(grid2, function(xx) do.call(dMixGpd, c(list(x = xx), args_mixgpd, list(log = 0L))), numeric(1))
        p2 <- vapply(grid2, function(xx) do.call(pMixGpd, c(list(q = xx), args_mixgpd, list(lower.tail = 1L, log.p = 0L))), numeric(1))

        if (any(!is.finite(p2))) {
          bad <- which(!is.finite(p2))[1]
          fail(paste0(
            name, " : pMixGpd produced non-finite at q=", grid2[bad],
            " (threshold=", u, ")."
          ))
        }

        .expect_density_ok(d2, label = paste0(name, " : dMixGpd"))
        .expect_cdf_ok(p2, label = paste0(name, " : pMixGpd"))

        r2 <- do.call(rMixGpd, c(list(n = 1L), args_mixgpd))
        .expect_rng_ok(r2, support = support, label = paste0(name, " : rMixGpd"))

        probs2 <- c(0.05, 0.2, 0.5, 0.9, 0.98)
        .roundtrip_qp(qMixGpd, pMixGpd, probs2, args_mixgpd, tol = tol_splice, label = paste0(name, " : mixGpd roundtrip"))

        succeed(paste0("mix+GPD checks ok: ", name))
      }

      # --- Scalar+GPD checks (if present) ---
      if (!is.null(dScalarGpd) && !is.null(pScalarGpd) && !is.null(rScalarGpd) && !is.null(qScalarGpd)) {
        expect_true(!is.null(args_scalargpd), info = paste0(name, " : args_scalargpd missing"))

        u <- args_scalargpd$threshold
        grid3 <- sort(unique(c(u - 0.1, u, u + 0.1, u + 1)))

        d3 <- vapply(grid3, function(xx) do.call(dScalarGpd, c(list(x = xx), args_scalargpd, list(log = 0L))), numeric(1))
        p3 <- vapply(grid3, function(xx) do.call(pScalarGpd, c(list(q = xx), args_scalargpd, list(lower.tail = 1L, log.p = 0L))), numeric(1))

        if (any(!is.finite(p3))) {
          bad <- which(!is.finite(p3))[1]
          fail(paste0(
            name, " : pScalarGpd produced non-finite at q=", grid3[bad],
            " (threshold=", u, ")."
          ))
        }

        .expect_density_ok(d3, label = paste0(name, " : dScalarGpd"))
        .expect_cdf_ok(p3, label = paste0(name, " : pScalarGpd"))

        r3 <- do.call(rScalarGpd, c(list(n = 1L), args_scalargpd))
        .expect_rng_ok(r3, support = support, label = paste0(name, " : rScalarGpd"))

        probs3 <- c(0.05, 0.2, 0.5, 0.9, 0.98)
        .roundtrip_qp(qScalarGpd, pScalarGpd, probs3, args_scalargpd, tol = tol_splice, label = paste0(name, " : scalarGpd roundtrip"))

        succeed(paste0("scalar+GPD checks ok: ", name))
      }

      invisible(TRUE)
    }

    # ---------- shared ----------
    # 3-component mixtures everywhere (per your rule)
    w3 <- c(0.50, 0.30, 0.20)

    gpd_tail <- list(
      threshold  = 2.0,
      tail_scale = 1.0,
      tail_shape = 0.2
    )

    # Gamma
    .check_mix_family(
      name = "Gamma",
      support = c(0, Inf),
      dMix = dGammaMix, pMix = pGammaMix, rMix = rGammaMix, qMix = qGammaMix,
      dMixGpd = dGammaMixGpd, pMixGpd = pGammaMixGpd, rMixGpd = rGammaMixGpd, qMixGpd = qGammaMixGpd,
      dScalarGpd = dGammaGpd, pScalarGpd = pGammaGpd, rScalarGpd = rGammaGpd, qScalarGpd = qGammaGpd,
      args_mix = list(w = w3, shape = c(2, 4, 6), scale = c(1, 1.5, 2.0)),
      args_mixgpd = c(list(w = w3, shape = c(2, 4, 6), scale = c(1, 1.5, 2.0)), gpd_tail),
      args_scalargpd = c(list(shape = 3, scale = 1), gpd_tail),
      grid = c(0.1, 0.5, 1, 2, 4, 6),
      tol_mix = 5e-4, tol_splice = 1e-3
    )

    # Normal
    .check_mix_family(
      name = "Normal",
      support = c(-Inf, Inf),
      dMix = dNormMix, pMix = pNormMix, rMix = rNormMix, qMix = qNormMix,
      dMixGpd = dNormMixGpd, pMixGpd = pNormMixGpd, rMixGpd = rNormMixGpd, qMixGpd = qNormMixGpd,
      dScalarGpd = dNormGpd, pScalarGpd = pNormGpd, rScalarGpd = rNormGpd, qScalarGpd = qNormGpd,
      args_mix = list(w = w3, mean = c(-1, 1, 3), sd = c(1, 0.8, 0.6)),
      args_mixgpd = c(list(w = w3, mean = c(-1, 1, 3), sd = c(1, 0.8, 0.6)), gpd_tail),
      args_scalargpd = c(list(mean = 0, sd = 1), gpd_tail),
      grid = c(-3, -1, 0, 1, 3, 5),
      tol_mix = 5e-4, tol_splice = 1e-3
    )

    # Lognormal
    .check_mix_family(
      name = "Lognormal",
      support = c(0, Inf),
      dMix = dLognormalMix, pMix = pLognormalMix, rMix = rLognormalMix, qMix = qLognormalMix,
      dMixGpd = dLognormalMixGpd, pMixGpd = pLognormalMixGpd, rMixGpd = rLognormalMixGpd, qMixGpd = qLognormalMixGpd,
      dScalarGpd = dLognormalGpd, pScalarGpd = pLognormalGpd, rScalarGpd = rLognormalGpd, qScalarGpd = qLognormalGpd,
      args_mix = list(w = w3, meanlog = c(-0.2, 0.6, 1.1), sdlog = c(0.7, 0.5, 0.4)),
      args_mixgpd = c(list(w = w3, meanlog = c(-0.2, 0.6, 1.1), sdlog = c(0.7, 0.5, 0.4)), gpd_tail),
      args_scalargpd = c(list(meanlog = 0, sdlog = 0.6), gpd_tail),
      grid = c(0.1, 0.5, 1, 2, 5, 10),
      tol_mix = 5e-4, tol_splice = 1e-3
    )

    # Laplace
    .check_mix_family(
      name = "Laplace",
      support = c(-Inf, Inf),
      dMix = dLaplaceMix, pMix = pLaplaceMix, rMix = rLaplaceMix, qMix = qLaplaceMix,
      dMixGpd = dLaplaceMixGpd, pMixGpd = pLaplaceMixGpd, rMixGpd = rLaplaceMixGpd, qMixGpd = qLaplaceMixGpd,
      dScalarGpd = dLaplaceGpd, pScalarGpd = pLaplaceGpd, rScalarGpd = rLaplaceGpd, qScalarGpd = qLaplaceGpd,
      args_mix = list(w = w3, location = c(-1, 1, 3), scale = c(1.2, 0.8, 0.6)),
      args_mixgpd = c(list(w = w3, location = c(-1, 1, 3), scale = c(1.2, 0.8, 0.6)), gpd_tail),
      args_scalargpd = c(list(location = 0, scale = 1), gpd_tail),
      grid = c(-3, -1, 0, 1, 3, 5),
      tol_mix = 5e-4, tol_splice = 1e-3
    )

    # Inverse Gaussian
    .check_mix_family(
      name = "InvGauss",
      support = c(0, Inf),
      dMix = dInvGaussMix, pMix = pInvGaussMix, rMix = rInvGaussMix, qMix = qInvGaussMix,
      dMixGpd = dInvGaussMixGpd, pMixGpd = pInvGaussMixGpd, rMixGpd = rInvGaussMixGpd, qMixGpd = qInvGaussMixGpd,
      dScalarGpd = dInvGaussGpd, pScalarGpd = pInvGaussGpd, rScalarGpd = rInvGaussGpd, qScalarGpd = qInvGaussGpd,
      args_mix = list(w = w3, mean = c(0.8, 1.6, 2.4), shape = c(1.5, 2.5, 3.5)),
      args_mixgpd = c(list(w = w3, mean = c(0.8, 1.6, 2.4), shape = c(1.5, 2.5, 3.5)), gpd_tail),
      args_scalargpd = c(list(mean = 1.5, shape = 2.0), gpd_tail),
      grid = c(0.1, 0.5, 1, 2, 4, 8),
      tol_mix = 1e-3, tol_splice = 2e-3
    )

    # Amoroso
    .check_mix_family(
      name = "Amoroso",
      support = c(0, Inf),
      dMix = dAmorosoMix, pMix = pAmorosoMix, rMix = rAmorosoMix, qMix = qAmorosoMix,
      dMixGpd = dAmorosoMixGpd, pMixGpd = pAmorosoMixGpd, rMixGpd = rAmorosoMixGpd, qMixGpd = qAmorosoMixGpd,
      dScalarGpd = dAmorosoGpd, pScalarGpd = pAmorosoGpd, rScalarGpd = rAmorosoGpd, qScalarGpd = qAmorosoGpd,
      args_mix = list(
        w = w3,
        loc = c(0, 0.3, 0.6),
        scale = c(1, 1, 1),
        shape1 = c(2, 3, 5),
        shape2 = c(1.0, 1.3, 1.6)
      ),
      args_mixgpd = c(list(
        w = w3,
        loc = c(0, 0.3, 0.6),
        scale = c(1, 1, 1),
        shape1 = c(2, 3, 5),
        shape2 = c(1.0, 1.3, 1.6)
      ), gpd_tail),
      args_scalargpd = c(list(loc = 0, scale = 1, shape1 = 2, shape2 = 1.3), gpd_tail),
      grid = c(0.01, 0.2, 0.8, 1.5, 3, 6),
      tol_mix = 1e-3, tol_splice = 2e-3
    )

    # Cauchy (mix only)
    .check_mix_family(
      name = "Cauchy",
      support = c(-Inf, Inf),
      dMix = dCauchyMix, pMix = pCauchyMix, rMix = rCauchyMix, qMix = qCauchyMix,
      args_mix = list(w = w3, location = c(-1, 1, 3), scale = c(1.2, 0.8, 0.6)),
      grid = c(-5, -2, -1, 0, 1, 3, 6),
      tol_mix = 1e-3
    )
  }
)
test_that("All nimble::nimbleFunction-defined objects in R/ compile", {
  # Keep this out of Tier A (cran) because it's expensive; allow it for CI + coverage.
  skip_if_not_test_level("ci")

  r_dir <- test_path("..", "..", "R")
  expect_true(dir.exists(r_dir), info = paste("Missing R dir:", r_dir))

  r_files <- list.files(r_dir, pattern = "\\.[Rr]$", full.names = TRUE)
  expect_true(length(r_files) > 0, info = "No R files found under R/")

  # Find:  NAME <- nimble::nimbleFunction(   OR   NAME <- nimbleFunction(
  find_defs <- function(files) {
    rx <- "^\\s*([.A-Za-z][.A-Za-z0-9_]*)\\s*<-\\s*(nimble::)?nimbleFunction\\s*\\("
    out <- list()
    for (ff in files) {
      lines <- readLines(ff, warn = FALSE)
      hit <- grepl(rx, lines)
      if (any(hit)) {
        idx <- which(hit)
        nm  <- sub(rx, "\\1", lines[idx], perl = TRUE)
        out[[ff]] <- data.frame(
          name = nm,
          file = basename(ff),
          line = idx,
          stringsAsFactors = FALSE
        )
      }
    }
    if (!length(out)) {
      return(data.frame(name = character(), file = character(), line = integer()))
    }
    do.call(rbind, out)
  }

  defs <- find_defs(r_files)
  expect_true(nrow(defs) > 0, info = "No nimbleFunction definitions found in R/")

  # Use the package namespace so dependencies among nimbleFunctions are resolved.
  ns <- tryCatch(asNamespace("CausalMixGPD"), error = function(e) NULL)
  if (is.null(ns)) {
    skip("CausalMixGPD namespace not available for compile checks.")
  }

  # Pretty status symbols (colored if crayon is available)
  tick  <- if (requireNamespace("crayon", quietly = TRUE))
    crayon::green("\u2714") else "\u2714"
  cross <- if (requireNamespace("crayon", quietly = TRUE))
    crayon::red("\u2718") else "\u2718"

  failures <- character(0)
  results <- data.frame(
    status = character(),
    name   = character(),
    where  = character(),
    detail = character(),
    stringsAsFactors = FALSE
  )

  for (i in seq_len(nrow(defs))) {
    nm <- defs$name[i]
    where <- paste0(defs$file[i], ":", defs$line[i])

    if (!exists(nm, envir = ns, inherits = FALSE)) {
      msg <- paste0("object not found in package namespace")
      results <- rbind(results, data.frame(
        status = cross, name = nm, where = where, detail = msg,
        stringsAsFactors = FALSE
      ))
      failures <- c(failures, paste0("- ", nm, " @ ", where, " :: ", msg))
      next
    }

    obj <- get(nm, envir = ns, inherits = FALSE)

    # If a nimbleFunction was wrapped for vectorization, compile its *_nf alias.
    if (!inherits(obj, "nimbleFunction")) {
      nf_name <- paste0(nm, "_nf")
      if (exists(nf_name, envir = ns, inherits = FALSE)) {
        obj <- get(nf_name, envir = ns, inherits = FALSE)
      }
    }

    err <- tryCatch({
      suppressWarnings(nimble::compileNimble(obj))
      NULL
    }, error = function(e) e)

    if (is.null(err)) {
      # Per-object PASS entry in testthat output
      expect_true(TRUE, info = paste0("compileNimble succeeded for ", nm, " @ ", where))

      results <- rbind(results, data.frame(
        status = tick, name = nm, where = where, detail = "",
        stringsAsFactors = FALSE
      ))
    } else {
      e1 <- strsplit(conditionMessage(err), "\n", fixed = TRUE)[[1]][1]
      cls <- paste(class(obj), collapse = "/")
      msg <- paste0("class=", cls, " :: ", e1)

      results <- rbind(results, data.frame(
        status = cross, name = nm, where = where, detail = msg,
        stringsAsFactors = FALSE
      ))
      failures <- c(failures, paste0("- ", nm, " @ ", where, " :: ", msg))
    }
  }

  # Print a compact summary table
  cat("\nCompilation summary (nimbleFunction -> compileNimble):\n")
  print(results[, c("status", "name", "where")], row.names = FALSE)

  # If there are failures, print details and fail the test
  if (length(failures) > 0) {
    cat("\nFailure details:\n")
    cat(paste0(failures, collapse = "\n"), "\n")
    fail(paste(
      "The following nimbleFunction-defined objects failed to compile:",
      paste(failures, collapse = "\n"),
      sep = "\n"
    ))
  }
})

# ===== END unit/test-kernels.R =====

# ===== BEGIN unit/test-s3-methods.R =====
# tests/testthat/test-s3-methods.R
# Unit tests for S3 print/summary/plot methods (04-S3-Methods.R)
# Uses minimal mock objects to avoid heavy MCMC runs

pkg <- "CausalMixGPD"
print_mixgpd_params <- getFromNamespace("print.mixgpd_params", pkg)
print_mixgpd_params_pair <- getFromNamespace("print.mixgpd_params_pair", pkg)
print_causalmixgpd_qte <- getFromNamespace("print.causalmixgpd_qte", pkg)
print_causalmixgpd_ate <- getFromNamespace("print.causalmixgpd_ate", pkg)
summary_causalmixgpd_qte <- getFromNamespace("summary.causalmixgpd_qte", pkg)
summary_causalmixgpd_ate <- getFromNamespace("summary.causalmixgpd_ate", pkg)
print_summary_causalmixgpd_qte <- getFromNamespace("print.summary.causalmixgpd_qte", pkg)
print_summary_causalmixgpd_ate <- getFromNamespace("print.summary.causalmixgpd_ate", pkg)
print_summary_causalmixgpd_ps_fit <- getFromNamespace("print.summary.causalmixgpd_ps_fit", pkg)
print_summary_causalmixgpd_causal_fit <- getFromNamespace("print.summary.causalmixgpd_causal_fit", pkg)
plot_mixgpd_predict <- getFromNamespace("plot.mixgpd_predict", pkg)
print_causalmixgpd_causal_predict_plots <- getFromNamespace("print.causalmixgpd_causal_predict_plots", pkg)
print_mixgpd_predict_plots <- getFromNamespace("print.mixgpd_predict_plots", pkg)
print_mixgpd_fit_plots <- getFromNamespace("print.mixgpd_fit_plots", pkg)
print_causalmixgpd_causal_fit_plots <- getFromNamespace("print.causalmixgpd_causal_fit_plots", pkg)
plot_mixgpd_fitted <- getFromNamespace("plot.mixgpd_fitted", pkg)
print_mixgpd_fitted_plots <- getFromNamespace("print.mixgpd_fitted_plots", pkg)
plot_causalmixgpd_qte <- getFromNamespace("plot.causalmixgpd_qte", pkg)
plot_causalmixgpd_ate <- getFromNamespace("plot.causalmixgpd_ate", pkg)

# ======================================================================
# Mock object constructors
# ======================================================================

make_mock_mixgpd_params <- function(empty = FALSE) {
  if (empty) {
    out <- list()
  } else {
    out <- list(
      alpha = 1.5,
      w = c(0.6, 0.3, 0.1),
      mu = c(1, 2, 3),
      sigma = c(0.5, 1, 1.5),
      beta_mu = matrix(c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6), nrow = 3, ncol = 2,
                       dimnames = list(paste0("comp", 1:3), c("x1", "x2")))
    )
  }
  class(out) <- "mixgpd_params"
  out
}

make_mock_mixgpd_params_pair <- function() {
  out <- list(
    treated = make_mock_mixgpd_params(),
    control = make_mock_mixgpd_params()
  )
  class(out) <- "mixgpd_params_pair"
  out
}

make_mock_qte <- function() {
  out <- list(
    probs = c(0.25, 0.5, 0.75),
    grid = c(0.25, 0.5, 0.75),
    n_pred = 5,
    level = 0.95,
    interval = "credible",
    meta = list(backend = list(trt = "sb", con = "sb"),
                kernel = list(trt = "normal", con = "normal"),
                GPD = list(trt = FALSE, con = FALSE)),
    ps = runif(5),
    x = matrix(rnorm(10), nrow = 5, ncol = 2),
    qte = list(
      fit = data.frame(
        id = rep(1:5, each = 3),
        index = rep(c(0.25, 0.5, 0.75), times = 5),
        estimate = rnorm(15),
        lower = rnorm(15) - 0.5,
        upper = rnorm(15) + 0.5
      )
    ),
    fit = matrix(rnorm(15), nrow = 5, ncol = 3),
    lower = matrix(rnorm(15) - 0.5, nrow = 5, ncol = 3),
    upper = matrix(rnorm(15) + 0.5, nrow = 5, ncol = 3),
    trt = list(fit = data.frame(estimate = rnorm(15), id = rep(1:5, 3), index = rep(c(0.25, 0.5, 0.75), each = 5))),
    con = list(fit = data.frame(estimate = rnorm(15), id = rep(1:5, 3), index = rep(c(0.25, 0.5, 0.75), each = 5)))
  )
  class(out) <- "causalmixgpd_qte"
  out
}

make_mock_ate <- function() {
  out <- list(
    n_pred = 5,
    level = 0.95,
    interval = "credible",
    nsim_mean = 200,
    meta = list(backend = list(trt = "sb", con = "sb"),
                kernel = list(trt = "normal", con = "normal"),
                GPD = list(trt = FALSE, con = FALSE)),
    ps = runif(5),
    x = matrix(rnorm(10), nrow = 5, ncol = 2),
    ate = list(
      fit = data.frame(
        id = 1:5,
        estimate = rnorm(5),
        lower = rnorm(5) - 0.5,
        upper = rnorm(5) + 0.5
      )
    ),
    fit = rnorm(5),
    lower = rnorm(5) - 0.5,
    upper = rnorm(5) + 0.5,
    trt = list(fit = data.frame(estimate = rnorm(5), id = 1:5)),
    con = list(fit = data.frame(estimate = rnorm(5), id = 1:5))
  )
  class(out) <- "causalmixgpd_ate"
  out
}

make_mock_mixgpd_predict <- function(type = "quantile") {
  if (type == "quantile") {
    out <- list(
      fit = data.frame(
        estimate = c(1, 2, 3),
        index = c(0.25, 0.5, 0.75),
        lower = c(0.5, 1.5, 2.5),
        upper = c(1.5, 2.5, 3.5)
      ),
      type = "quantile",
      grid = c(0.25, 0.5, 0.75)
    )
  } else if (type == "sample") {
    out <- list(
      fit = rnorm(100),
      type = "sample",
      grid = NULL
    )
  } else if (type == "mean") {
    out <- list(
      fit = data.frame(estimate = 5.2, lower = 4.8, upper = 5.6),
      type = "mean",
      grid = NULL
    )
  } else if (type == "density") {
    y_grid <- seq(0, 5, length.out = 20)
    out <- list(
      fit = data.frame(
        id = rep(1, 20),
        y = y_grid,
        density = dnorm(y_grid, mean = 2.5, sd = 1),
        lower = dnorm(y_grid, mean = 2.5, sd = 1) - 0.05,
        upper = dnorm(y_grid, mean = 2.5, sd = 1) + 0.05
      ),
      type = "density",
      grid = y_grid
    )
  } else if (type == "survival") {
    y_grid <- seq(0, 5, length.out = 20)
    out <- list(
      fit = data.frame(
        id = rep(1, 20),
        y = y_grid,
        survival = 1 - pnorm(y_grid, mean = 2.5, sd = 1),
        lower = pmax(0, 1 - pnorm(y_grid, mean = 2.5, sd = 1) - 0.05),
        upper = pmin(1, 1 - pnorm(y_grid, mean = 2.5, sd = 1) + 0.05)
      ),
      type = "survival",
      grid = y_grid
    )
  } else if (type == "location") {
    out <- list(
      fit = data.frame(
        mean = 5.0, mean_lower = 4.5, mean_upper = 5.5,
        median = 4.9, median_lower = 4.4, median_upper = 5.4
      ),
      type = "location",
      grid = NULL
    )
  } else {
    stop("Unknown mock predict type")
  }
  class(out) <- "mixgpd_predict"
  out
}

make_mock_causal_predict_plots <- function() {
  out <- list(
    trt_control = ggplot2::ggplot() + ggplot2::ggtitle("trt_control"),
    treatment_effect = ggplot2::ggplot() + ggplot2::ggtitle("treatment_effect")
  )
  class(out) <- c("causalmixgpd_causal_predict_plots", "list")
  out
}

make_mock_mixgpd_predict_plots <- function() {
  p <- ggplot2::ggplot() + ggplot2::ggtitle("predict_plot")
  class(p) <- c("mixgpd_predict_plots", class(p))
  p
}

make_mock_mixgpd_fit_plots <- function() {
  out <- list(
    traceplot = ggplot2::ggplot() + ggplot2::ggtitle("traceplot"),
    density = ggplot2::ggplot() + ggplot2::ggtitle("density")
  )
  class(out) <- c("mixgpd_fit_plots", "list")
  out
}

make_mock_mixgpd_fitted <- function() {
  n <- 20
  fit_vals <- rnorm(n, mean = 5)
  y_obs <- fit_vals + rnorm(n, sd = 0.5)
  out <- data.frame(
    fit = fit_vals,
    lower = fit_vals - 0.5,
    upper = fit_vals + 0.5,
    residuals = y_obs - fit_vals
  )
  class(out) <- c("mixgpd_fitted", "data.frame")

  # Attach a minimal mock object

mock_obj <- list(
    data = list(y = y_obs),
    spec = list(meta = list(backend = "sb", kernel = "normal", GPD = FALSE))
  )
  attr(out, "object") <- mock_obj
  attr(out, "level") <- 0.95
  attr(out, "interval") <- "credible"
  out
}

make_mock_causal_fit_plots <- function() {
  treated_plots <- list(
    traceplot = ggplot2::ggplot() + ggplot2::ggtitle("treated_trace"),
    density = ggplot2::ggplot() + ggplot2::ggtitle("treated_density")
  )
  class(treated_plots) <- c("mixgpd_fit_plots", "list")

  control_plots <- list(
    traceplot = ggplot2::ggplot() + ggplot2::ggtitle("control_trace"),
    density = ggplot2::ggplot() + ggplot2::ggtitle("control_density")
  )
  class(control_plots) <- c("mixgpd_fit_plots", "list")

  out <- list(treated = treated_plots, control = control_plots)
  class(out) <- c("causalmixgpd_causal_fit_plots", "list")
  out
}

make_mock_summary_qte <- function() {
  qte_obj <- make_mock_qte()
  out <- list(
    overall = list(
      n_pred = 5,
      n_quantiles = 3,
      quantiles = c(0.25, 0.5, 0.75),
      level = 0.95,
      interval = "credible",
      has_covariates = TRUE,
      ps_used = TRUE
    ),
    quantile_summary = data.frame(
      quantile = c(0.25, 0.5, 0.75),
      mean_qte = c(0.1, 0.2, 0.3),
      median_qte = c(0.1, 0.2, 0.3),
      min_qte = c(-0.1, 0, 0.1),
      max_qte = c(0.3, 0.4, 0.5),
      sd_qte = c(0.1, 0.1, 0.1)
    ),
    ci_summary = list(
      mean_width = 0.5,
      median_width = 0.5,
      min_width = 0.3,
      max_width = 0.7
    ),
    meta = list(
      backend = list(trt = "sb", con = "sb"),
      kernel = list(trt = "normal", con = "normal"),
      GPD = list(trt = FALSE, con = FALSE)
    ),
    object = qte_obj
  )
  class(out) <- "summary.causalmixgpd_qte"
  out
}

make_mock_summary_ate <- function() {
  ate_obj <- make_mock_ate()
  out <- list(
    overall = list(
      n_pred = 5,
      level = 0.95,
      interval = "credible",
      nsim_mean = 200,
      has_covariates = TRUE,
      ps_used = TRUE
    ),
    ate_stats = list(
      mean_ate = 0.5,
      median_ate = 0.45,
      min_ate = 0.1,
      max_ate = 0.9,
      sd_ate = 0.2
    ),
    ci_summary = list(
      mean_width = 0.5,
      median_width = 0.5,
      min_width = 0.3,
      max_width = 0.7
    ),
    meta = list(
      backend = list(trt = "sb", con = "sb"),
      kernel = list(trt = "normal", con = "normal"),
      GPD = list(trt = FALSE, con = FALSE)
    ),
    object = ate_obj
  )
  class(out) <- "summary.causalmixgpd_ate"
  out
}

make_mock_mixgpd_fitted_plots <- function() {
  out <- list(
    observed_fitted_plot = ggplot2::ggplot() + ggplot2::ggtitle("obs_fit"),
    residual_plot = ggplot2::ggplot() + ggplot2::ggtitle("residuals")
  )
  class(out) <- c("mixgpd_fitted_plots", "list")
  out
}

# ======================================================================
# print.mixgpd_params tests
# ======================================================================

test_that("print.mixgpd_params works with non-empty params", {
  params <- make_mock_mixgpd_params()
  expect_output(print_mixgpd_params(params), "Posterior mean parameters")
  expect_output(print_mixgpd_params(params), "alpha")
  expect_output(print_mixgpd_params(params), "beta_mu")
})

test_that("print.mixgpd_params works with empty params", {
  params <- make_mock_mixgpd_params(empty = TRUE)
  expect_output(print_mixgpd_params(params), "empty")
})

test_that("print.mixgpd_params respects digits argument", {
  params <- make_mock_mixgpd_params()
  expect_output(print_mixgpd_params(params, digits = 2), regexp = "Posterior mean parameters")
})

# ======================================================================
# print.mixgpd_params_pair tests
# ======================================================================

test_that("print.mixgpd_params_pair works", {
  params_pair <- make_mock_mixgpd_params_pair()
  expect_output(print_mixgpd_params_pair(params_pair), "Posterior mean parameters \\(causal\\)")
  expect_output(print_mixgpd_params_pair(params_pair), "treated")
  expect_output(print_mixgpd_params_pair(params_pair), "control")
})

# ======================================================================
# print.causalmixgpd_qte tests
# ======================================================================

test_that("print.causalmixgpd_qte works", {
  qte <- make_mock_qte()
  expect_output(print_causalmixgpd_qte(qte), "QTE \\(Quantile Treatment Effect\\)")
  expect_output(print_causalmixgpd_qte(qte), "Prediction points")
  expect_output(print_causalmixgpd_qte(qte), "Quantile grid")
})

test_that("print.causalmixgpd_qte respects max_rows", {
  qte <- make_mock_qte()
  expect_output(print_causalmixgpd_qte(qte, max_rows = 2), "more rows")
})

test_that("print.causalmixgpd_qte uses estimand-specific labels", {
  qte <- make_mock_qte()
  qte$type <- "cqte"
  expect_output(print_causalmixgpd_qte(qte), "CQTE")
  qte$type <- "qtt"
  expect_output(print_causalmixgpd_qte(qte), "QTT")
})

# ======================================================================
# print.causalmixgpd_ate tests
# ======================================================================

test_that("print.causalmixgpd_ate works", {
  ate <- make_mock_ate()
  expect_output(print_causalmixgpd_ate(ate), "ATE \\(Average Treatment Effect\\)")
  expect_output(print_causalmixgpd_ate(ate), "Prediction points")
})

test_that("print.causalmixgpd_ate respects max_rows", {
  ate <- make_mock_ate()
  expect_output(print_causalmixgpd_ate(ate, max_rows = 2), "more rows")
})

test_that("print.causalmixgpd_ate uses estimand-specific labels", {
  ate <- make_mock_ate()
  ate$type <- "cate"
  expect_output(print_causalmixgpd_ate(ate), "CATE")
  ate$type <- "att"
  expect_output(print_causalmixgpd_ate(ate), "ATT")
})

# ======================================================================
# summary.causalmixgpd_qte tests
# ======================================================================

test_that("summary.causalmixgpd_qte returns proper structure", {
  qte <- make_mock_qte()
  summ <- summary_causalmixgpd_qte(qte)
  expect_s3_class(summ, "summary.causalmixgpd_qte")
  expect_true("overall" %in% names(summ))
  expect_true("quantile_summary" %in% names(summ))
})

test_that("summary.causalmixgpd_qte prefers profile over id for conditional tables", {
  qte <- make_mock_qte()
  qte$type <- "cqte"
  qte$profile <- paste("Profile", 1:5)
  qte$qte$fit$profile <- rep(qte$profile, each = 3)

  summ <- summary_causalmixgpd_qte(qte)
  expect_true("profile" %in% names(summ$effect_table))
  expect_false("id" %in% names(summ$effect_table))
})

# ======================================================================
# summary.causalmixgpd_ate tests
# ======================================================================

test_that("summary.causalmixgpd_ate returns proper structure", {
  ate <- make_mock_ate()
  summ <- summary_causalmixgpd_ate(ate)
  expect_s3_class(summ, "summary.causalmixgpd_ate")
  expect_true("overall" %in% names(summ))
  expect_true("ate_stats" %in% names(summ))
})

test_that("summary.causalmixgpd_ate prefers profile over id for conditional tables", {
  ate <- make_mock_ate()
  ate$type <- "cate"
  ate$profile <- paste("Profile", 1:5)
  ate$ate$fit$profile <- ate$profile

  summ <- summary_causalmixgpd_ate(ate)
  expect_true("profile" %in% names(summ$effect_table))
  expect_false("id" %in% names(summ$effect_table))
})

# ======================================================================
# print.summary.causalmixgpd_qte tests
# ======================================================================

test_that("print.summary.causalmixgpd_qte works", {
  summ <- make_mock_summary_qte()
  expect_output(print_summary_causalmixgpd_qte(summ), "QTE Summary")
  expect_output(print_summary_causalmixgpd_qte(summ), "Prediction points")
  expect_output(print_summary_causalmixgpd_qte(summ), "Quantile grid")
  expect_output(print_summary_causalmixgpd_qte(summ), "Model specification")
})

test_that("print.summary.causalmixgpd_qte respects digits", {
  summ <- make_mock_summary_qte()
  expect_output(print_summary_causalmixgpd_qte(summ, digits = 2), "QTE Summary")
})

test_that("print.summary.causalmixgpd_qte uses estimand-specific labels", {
  summ <- make_mock_summary_qte()
  summ$object$type <- "cqte"
  expect_output(print_summary_causalmixgpd_qte(summ), "CQTE Summary")
})

# ======================================================================
# print.summary.causalmixgpd_ate tests
# ======================================================================

test_that("print.summary.causalmixgpd_ate works", {
  summ <- make_mock_summary_ate()
  expect_output(print_summary_causalmixgpd_ate(summ), "ATE Summary")
  expect_output(print_summary_causalmixgpd_ate(summ), "Prediction points")
  expect_output(print_summary_causalmixgpd_ate(summ), "Model specification")
})

test_that("print.summary.causalmixgpd_ate respects digits", {
  summ <- make_mock_summary_ate()
  expect_output(print_summary_causalmixgpd_ate(summ, digits = 2), "ATE Summary")
})

test_that("print.summary.causalmixgpd_ate uses estimand-specific labels", {
  summ <- make_mock_summary_ate()
  summ$object$type <- "att"
  expect_output(print_summary_causalmixgpd_ate(summ), "ATT Summary")
})

# ======================================================================
# plot.mixgpd_predict tests
# ======================================================================

test_that("plot.mixgpd_predict works for quantile", {
  skip_if_not_installed("ggplot2")
  pred <- make_mock_mixgpd_predict("quantile")
  p <- plot_mixgpd_predict(pred)
  expect_true(inherits(p, "gg") || inherits(p, "ggplot"))
})

test_that("plot.mixgpd_predict works for sample", {
  skip_if_not_installed("ggplot2")
  pred <- make_mock_mixgpd_predict("sample")
  p <- plot_mixgpd_predict(pred)
  expect_true(inherits(p, "gg") || inherits(p, "ggplot"))
})

test_that("plot.mixgpd_predict works for mean", {
  skip_if_not_installed("ggplot2")
  pred <- make_mock_mixgpd_predict("mean")
  p <- plot_mixgpd_predict(pred)
  expect_true(inherits(p, "gg") || inherits(p, "ggplot"))
})

test_that("plot.mixgpd_predict works for density", {
  skip_if_not_installed("ggplot2")
  pred <- make_mock_mixgpd_predict("density")
  p <- plot_mixgpd_predict(pred)
  expect_true(inherits(p, "gg") || inherits(p, "ggplot"))
})

test_that("plot.mixgpd_predict works for survival", {
  skip_if_not_installed("ggplot2")
  pred <- make_mock_mixgpd_predict("survival")
  p <- plot_mixgpd_predict(pred)
  expect_true(inherits(p, "gg") || inherits(p, "ggplot"))
})

test_that("plot.mixgpd_predict works for location", {
  skip_if_not_installed("ggplot2")
  pred <- make_mock_mixgpd_predict("location")
  p <- plot_mixgpd_predict(pred)
  expect_true(inherits(p, "gg") || inherits(p, "ggplot") || is.list(p))
})

test_that("plot.mixgpd_predict errors on non-list input", {
  expect_error(plot_mixgpd_predict("not_a_list"), "must be a prediction object")
})

test_that("plot.mixgpd_predict errors on missing type", {
  pred <- list(fit = data.frame(x = 1))
  class(pred) <- "mixgpd_predict"
  expect_error(plot_mixgpd_predict(pred), "missing 'type'")
})

# ======================================================================
# print.causalmixgpd_causal_predict_plots tests
# ======================================================================

test_that("print.causalmixgpd_causal_predict_plots works", {
  skip_if_not_installed("ggplot2")
  plots <- make_mock_causal_predict_plots()
  expect_output(print_causalmixgpd_causal_predict_plots(plots), regexp = ".")
})

# ======================================================================
# print.mixgpd_predict_plots tests
# ======================================================================

test_that("print.mixgpd_predict_plots works", {
  skip_if_not_installed("ggplot2")
  p <- make_mock_mixgpd_predict_plots()
  # Should not error
  expect_silent(invisible(capture.output(print_mixgpd_predict_plots(p))))
})

# ======================================================================
# print.mixgpd_fit_plots tests
# ======================================================================

test_that("print.mixgpd_fit_plots works", {
  skip_if_not_installed("ggplot2")
  plots <- make_mock_mixgpd_fit_plots()
  expect_output(print_mixgpd_fit_plots(plots), "traceplot")
  expect_output(print_mixgpd_fit_plots(plots), "density")
})

# ======================================================================
# print.causalmixgpd_causal_fit_plots tests
# ======================================================================

test_that("print.causalmixgpd_causal_fit_plots works", {
  skip_if_not_installed("ggplot2")
  plots <- make_mock_causal_fit_plots()
  expect_output(print_causalmixgpd_causal_fit_plots(plots), "treated")
  expect_output(print_causalmixgpd_causal_fit_plots(plots), "control")
})

# ======================================================================
# plot.mixgpd_fitted tests
# ======================================================================

test_that("plot.mixgpd_fitted returns plot list", {
  skip_if_not_installed("ggplot2")
  fitted_obj <- make_mock_mixgpd_fitted()
  result <- plot_mixgpd_fitted(fitted_obj)
  expect_s3_class(result, "mixgpd_fitted_plots")
  expect_true("observed_fitted_plot" %in% names(result))
  expect_true("residual_plot" %in% names(result))
})

# ======================================================================
# print.mixgpd_fitted_plots tests
# ======================================================================

test_that("print.mixgpd_fitted_plots works", {
  skip_if_not_installed("ggplot2")
  plots <- make_mock_mixgpd_fitted_plots()
  # Should not error (print produces ggplot output but not text)
  expect_silent(invisible(capture.output(print_mixgpd_fitted_plots(plots))))
})

# ======================================================================
# plot.causalmixgpd_qte tests
# ======================================================================

test_that("plot.causalmixgpd_qte returns both plots by default", {
  skip_if_not_installed("ggplot2")
  qte <- make_mock_qte()
  result <- plot_causalmixgpd_qte(qte)
  expect_s3_class(result, "causalmixgpd_causal_predict_plots")
  expect_true("trt_control" %in% names(result))
  expect_true("treatment_effect" %in% names(result))
})

test_that("plot.causalmixgpd_qte defaults CQTE to effect plot faceted by id", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$type <- "cqte"
  qte$n_pred <- 3L
  qte$ps <- c(0.2, 0.5, 0.8)
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(
    0.1, 0.2, 0.3,
    0.2, 0.3, 0.4,
    0.3, 0.4, 0.5
  ), nrow = 3L, byrow = TRUE)
  qte$lower <- qte$fit - 0.1
  qte$upper <- qte$fit + 0.1

  result <- plot_causalmixgpd_qte(qte)
  expect_true(inherits(result, "gg") || inherits(result, "ggplot"))
  expect_true(is.factor(result$data$x_plot))
  expect_equal(levels(result$data$x_plot), c("0.25", "0.5", "0.75"))

  build <- ggplot2::ggplot_build(result)
  line_data <- build$data[[1]]
  panel_unique_x <- vapply(split(line_data$x, line_data$PANEL), function(z) length(unique(z)), integer(1))

  expect_true(all(panel_unique_x == length(probs)))
  expect_equal(sort(unique(as.numeric(line_data$x))), seq_along(probs), tolerance = 1e-12)
})

test_that("plot.causalmixgpd_qte effect type works", {
  skip_if_not_installed("ggplot2")
  qte <- make_mock_qte()
  result <- plot_causalmixgpd_qte(qte, type = "effect")
  expect_true(inherits(result, "gg") || inherits(result, "ggplot"))
})

test_that("plot.causalmixgpd_qte marginal effect uses quantile levels on x-axis", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$n_pred <- 1L
  qte$ps <- NULL
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(0.1, 0.2, 0.3), nrow = 1L)
  qte$lower <- matrix(c(-0.1, 0.0, 0.1), nrow = 1L)
  qte$upper <- matrix(c(0.3, 0.4, 0.5), nrow = 1L)
  qte$qte$fit <- data.frame(
    id = 1L,
    index = probs,
    estimate = c(0.1, 0.2, 0.3),
    lower = c(-0.1, 0.0, 0.1),
    upper = c(0.3, 0.4, 0.5)
  )
  qte$trt$fit <- data.frame(id = 1L, index = probs, estimate = c(1.0, 1.3, 1.6))
  qte$con$fit <- data.frame(id = 1L, index = probs, estimate = c(0.9, 1.1, 1.3))

  result <- plot_causalmixgpd_qte(qte, type = "effect")
  expect_true(is.factor(result$data$x_plot))
  expect_equal(levels(result$data$x_plot), c("0.25", "0.5", "0.75"))
  build <- ggplot2::ggplot_build(result)
  line_data <- build$data[[1]]

  expect_equal(as.numeric(line_data$x), seq_along(probs), tolerance = 1e-12)
  expect_equal(length(unique(line_data$PANEL)), 1L)
})

test_that("plot.causalmixgpd_qte effect uses pointwise error bars instead of ribbons", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$n_pred <- 1L
  qte$ps <- NULL
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(0.1, 0.2, 0.3), nrow = 1L)
  qte$lower <- matrix(c(-0.1, 0.0, 0.1), nrow = 1L)
  qte$upper <- matrix(c(0.3, 0.4, 0.5), nrow = 1L)
  qte$trt$fit <- data.frame(id = 1L, index = probs, estimate = c(1.0, 1.3, 1.6))
  qte$con$fit <- data.frame(id = 1L, index = probs, estimate = c(0.9, 1.1, 1.3))

  result <- plot_causalmixgpd_qte(qte, type = "effect")
  geom_classes <- vapply(result[["layers"]], function(layer) class(layer[["geom"]])[1], character(1))

  expect_true("GeomErrorbar" %in% geom_classes)
  expect_false("GeomRibbon" %in% geom_classes)
})

test_that("plot.causalmixgpd_qte effect facets by id using quantile levels on x-axis", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$n_pred <- 3L
  qte$ps <- c(0.2, 0.5, 0.8)
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(
    0.1, 0.2, 0.3,
    0.2, 0.3, 0.4,
    0.3, 0.4, 0.5
  ), nrow = 3L, byrow = TRUE)
  qte$lower <- qte$fit - 0.1
  qte$upper <- qte$fit + 0.1

  result <- plot_causalmixgpd_qte(qte, type = "effect", facet_by = "id")
  expect_true(is.factor(result$data$x_plot))
  expect_equal(levels(result$data$x_plot), c("0.25", "0.5", "0.75"))
  build <- ggplot2::ggplot_build(result)
  line_data <- build$data[[1]]
  panel_unique_x <- vapply(split(line_data$x, line_data$PANEL), function(z) length(unique(z)), integer(1))

  expect_true(all(panel_unique_x == length(probs)))
  expect_equal(sort(unique(as.numeric(line_data$x))), seq_along(probs), tolerance = 1e-12)
})

test_that("plot.causalmixgpd_qte arms type works", {
  skip_if_not_installed("ggplot2")
  qte <- make_mock_qte()
  result <- plot_causalmixgpd_qte(qte, type = "arms")
  expect_true(inherits(result, "gg") || inherits(result, "ggplot"))
})

test_that("plot.causalmixgpd_qte arms use pointwise error bars instead of ribbons", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$n_pred <- 1L
  qte$ps <- NULL
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(0.1, 0.2, 0.3), nrow = 1L)
  qte$lower <- matrix(c(-0.1, 0.0, 0.1), nrow = 1L)
  qte$upper <- matrix(c(0.3, 0.4, 0.5), nrow = 1L)
  qte$trt$fit <- data.frame(
    id = 1L,
    index = probs,
    estimate = c(1.0, 1.3, 1.6),
    lower = c(0.8, 1.1, 1.4),
    upper = c(1.2, 1.5, 1.8)
  )
  qte$con$fit <- data.frame(
    id = 1L,
    index = probs,
    estimate = c(0.9, 1.1, 1.3),
    lower = c(0.7, 0.9, 1.1),
    upper = c(1.1, 1.3, 1.5)
  )

  result <- plot_causalmixgpd_qte(qte, type = "arms")
  geom_classes <- vapply(result[["layers"]], function(layer) class(layer[["geom"]])[1], character(1))

  expect_true("GeomErrorbar" %in% geom_classes)
  expect_false("GeomRibbon" %in% geom_classes)
})

test_that("plot.causalmixgpd_qte arms facets by id using quantile levels on x-axis", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$n_pred <- 3L
  qte$ps <- c(0.2, 0.5, 0.8)
  qte$probs <- probs
  qte$grid <- probs
  qte$trt$fit <- data.frame(
    id = rep(1:3, each = 3),
    index = rep(probs, times = 3),
    estimate = c(1.0, 1.2, 1.4, 1.1, 1.3, 1.5, 1.2, 1.4, 1.6),
    lower = c(0.9, 1.1, 1.3, 1.0, 1.2, 1.4, 1.1, 1.3, 1.5),
    upper = c(1.1, 1.3, 1.5, 1.2, 1.4, 1.6, 1.3, 1.5, 1.7)
  )
  qte$con$fit <- data.frame(
    id = rep(1:3, each = 3),
    index = rep(probs, times = 3),
    estimate = c(0.8, 1.0, 1.2, 0.9, 1.1, 1.3, 1.0, 1.2, 1.4),
    lower = c(0.7, 0.9, 1.1, 0.8, 1.0, 1.2, 0.9, 1.1, 1.3),
    upper = c(0.9, 1.1, 1.3, 1.0, 1.2, 1.4, 1.1, 1.3, 1.5)
  )

  result <- plot_causalmixgpd_qte(qte, type = "arms", facet_by = "id")
  expect_true(is.factor(result$data$x_plot))
  expect_equal(levels(result$data$x_plot), c("0.25", "0.5", "0.75"))
  build <- ggplot2::ggplot_build(result)
  line_data <- build$data[[1]]
  panel_unique_x <- vapply(split(line_data$x, line_data$PANEL), function(z) length(unique(z)), integer(1))

  expect_true(all(panel_unique_x == length(probs)))
  expect_equal(sort(unique(as.numeric(line_data$x))), seq_along(probs), tolerance = 1e-12)
})

test_that("plot.causalmixgpd_qte single quantile mentions tau in x-axis label", {
  skip_if_not_installed("ggplot2")
  qte <- make_mock_qte()
  qte$type <- "cqte"
  qte$probs <- 0.5
  qte$grid <- 0.5
  qte$n_pred <- 3L
  qte$ps <- c(0.2, 0.5, 0.8)
  qte$fit <- matrix(c(0.1, 0.2, 0.3), nrow = 3L, ncol = 1L)
  qte$lower <- qte$fit - 0.1
  qte$upper <- qte$fit + 0.1
  qte$trt$fit <- data.frame(id = 1:3, index = 0.5, estimate = c(1.0, 1.2, 1.4))
  qte$con$fit <- data.frame(id = 1:3, index = 0.5, estimate = c(0.9, 1.0, 1.1))

  result <- plot_causalmixgpd_qte(qte)
  expect_match(result$labels$x, "\u03C4 = 0.5", fixed = TRUE)
})

test_that("plot.causalmixgpd_qte applies categorical quantile axis and error bars to QTT", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$type <- "qtt"
  qte$n_pred <- 2L
  qte$ps <- c(0.3, 0.7)
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(0.1, 0.15, 0.2, 0.12, 0.18, 0.24), nrow = 2L, byrow = TRUE)
  qte$lower <- qte$fit - 0.05
  qte$upper <- qte$fit + 0.05
  qte$qte$fit <- data.frame(
    id = rep(1:2, each = 3),
    index = rep(probs, times = 2),
    estimate = c(0.1, 0.15, 0.2, 0.12, 0.18, 0.24),
    lower = c(0.05, 0.10, 0.15, 0.07, 0.13, 0.19),
    upper = c(0.15, 0.20, 0.25, 0.17, 0.23, 0.29)
  )

  result <- plot_causalmixgpd_qte(qte, type = "effect", facet_by = "id")
  geom_classes <- vapply(result$layers, function(layer) class(layer$geom)[1], character(1))

  expect_true(is.factor(result$data$x_plot))
  expect_equal(levels(result$data$x_plot), c("0.25", "0.5", "0.75"))
  expect_equal(result$labels$title, "QTT")
  expect_equal(result$labels$y, "Quantile Treatment Effect on the Treated")
  expect_true("GeomErrorbar" %in% geom_classes)
  expect_false("GeomRibbon" %in% geom_classes)
})

test_that("plot.causalmixgpd_qte keeps effect intervals aligned with estimates", {
  skip_if_not_installed("ggplot2")
  probs <- c(0.25, 0.5, 0.75)
  qte <- make_mock_qte()
  qte$type <- "cqte"
  qte$n_pred <- 2L
  qte$ps <- c(0.2, 0.8)
  qte$probs <- probs
  qte$grid <- probs
  qte$fit <- matrix(c(11, 12, 13, 21, 22, 23), nrow = 2L, byrow = TRUE)
  qte$lower <- matrix(c(101, 102, 103, 201, 202, 203), nrow = 2L, byrow = TRUE)
  qte$upper <- qte$lower + 0.5
  qte$qte$fit <- data.frame(
    id = rep(1:2, each = 3),
    index = rep(probs, times = 2),
    estimate = c(11, 12, 13, 21, 22, 23),
    lower = c(101, 102, 103, 201, 202, 203),
    upper = c(101.5, 102.5, 103.5, 201.5, 202.5, 203.5)
  )

  result <- plot_causalmixgpd_qte(qte, type = "effect", facet_by = "id")
  plotted <- result$data[order(result$data$id, result$data$index), c("id", "index", "estimate", "lower", "upper")]
  expected <- qte$qte$fit[order(qte$qte$fit$id, qte$qte$fit$index), c("id", "index", "estimate", "lower", "upper")]

  expect_equal(unname(as.matrix(plotted)), unname(as.matrix(expected)))
})

# ======================================================================
# plot.causalmixgpd_ate tests
# ======================================================================

test_that("plot.causalmixgpd_ate returns both plots by default", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  result <- plot_causalmixgpd_ate(ate)
  expect_s3_class(result, "causalmixgpd_causal_predict_plots")
  expect_true("trt_control" %in% names(result))
  expect_true("treatment_effect" %in% names(result))
})

test_that("plot.causalmixgpd_ate defaults CATE to effect plot", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  ate$type <- "cate"
  ate$ps <- c(0.2, 0.4, 0.6, 0.8, 0.9)
  ate$fit <- c(0.1, 0.15, 0.2, 0.25, 0.3)
  ate$lower <- ate$fit - 0.05
  ate$upper <- ate$fit + 0.05

  result <- plot_causalmixgpd_ate(ate)
  expect_true(inherits(result, "gg") || inherits(result, "ggplot"))
  expect_true(is.factor(result$data$x_plot))
  expect_equal(result$labels$x, "Index")
})

test_that("plot.causalmixgpd_ate uses profile labels for CATE when available", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  ate$type <- "cate"
  ate$profile <- paste("Profile", 1:5)
  ate$ate$fit$profile <- ate$profile

  result <- plot_causalmixgpd_ate(ate)
  expect_true(is.factor(result$data$x_plot))
  expect_equal(levels(result$data$x_plot), ate$profile)
  expect_equal(result$labels$x, "Profile")
})

test_that("plot.causalmixgpd_ate effect type works", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  result <- plot_causalmixgpd_ate(ate, type = "effect")
  expect_true(inherits(result, "gg") || inherits(result, "ggplot"))
})

test_that("plot.causalmixgpd_ate marginal effect returns a single point", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  ate$n_pred <- 1L
  ate$ps <- NULL
  ate$fit <- 0.25
  ate$lower <- -0.10
  ate$upper <- 0.45
  ate$ate$fit <- data.frame(id = 1L, estimate = 0.25, lower = -0.10, upper = 0.45)
  ate$trt$fit <- data.frame(id = 1L, estimate = 1.20, lower = 1.00, upper = 1.35)
  ate$con$fit <- data.frame(id = 1L, estimate = 0.95, lower = 0.80, upper = 1.10)

  result <- plot_causalmixgpd_ate(ate, type = "effect")
  build <- ggplot2::ggplot_build(result)

  expect_equal(nrow(build$data[[1]]), 1L)
  expect_equal(as.numeric(build$data[[1]]$y), 0.25, tolerance = 1e-12)
})

test_that("plot.causalmixgpd_ate uses pointwise error bars instead of ribbons", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  ate$type <- "cate"
  ate$ps <- c(0.2, 0.4, 0.6, 0.8, 0.9)
  ate$fit <- c(0.1, 0.15, 0.2, 0.25, 0.3)
  ate$lower <- ate$fit - 0.05
  ate$upper <- ate$fit + 0.05
  ate$trt$fit <- data.frame(
    id = 1:5,
    estimate = c(1.0, 1.1, 1.2, 1.3, 1.4),
    lower = c(0.9, 1.0, 1.1, 1.2, 1.3),
    upper = c(1.1, 1.2, 1.3, 1.4, 1.5)
  )
  ate$con$fit <- data.frame(
    id = 1:5,
    estimate = c(0.9, 0.95, 1.0, 1.05, 1.1),
    lower = c(0.8, 0.85, 0.9, 0.95, 1.0),
    upper = c(1.0, 1.05, 1.1, 1.15, 1.2)
  )

  effect_plot <- plot_causalmixgpd_ate(ate, type = "effect")
  effect_geoms <- vapply(effect_plot[["layers"]], function(layer) class(layer[["geom"]])[1], character(1))
  expect_true("GeomErrorbar" %in% effect_geoms)
  expect_false("GeomRibbon" %in% effect_geoms)

  arms_plot <- plot_causalmixgpd_ate(ate, type = "arms")
  arms_geoms <- vapply(arms_plot[["layers"]], function(layer) class(layer[["geom"]])[1], character(1))
  expect_true("GeomErrorbar" %in% arms_geoms)
  expect_false("GeomRibbon" %in% arms_geoms)
})

test_that("plot.causalmixgpd_ate applies shared effect styling to ATT", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  ate$type <- "att"
  ate$ps <- c(0.25, 0.5, 0.75)
  ate$fit <- c(0.1, 0.18, 0.26)
  ate$lower <- ate$fit - 0.04
  ate$upper <- ate$fit + 0.04

  result <- plot_causalmixgpd_ate(ate, type = "effect")
  geom_classes <- vapply(result$layers, function(layer) class(layer$geom)[1], character(1))

  expect_equal(result$labels$title, "ATT")
  expect_equal(result$labels$y, "Average Treatment Effect on the Treated")
  expect_true("GeomErrorbar" %in% geom_classes)
  expect_false("GeomRibbon" %in% geom_classes)
})

test_that("plot.causalmixgpd_ate arms type works", {
  skip_if_not_installed("ggplot2")
  ate <- make_mock_ate()
  result <- plot_causalmixgpd_ate(ate, type = "arms")
  expect_true(inherits(result, "gg") || inherits(result, "ggplot"))
})

# ======================================================================
# Error handling
# ======================================================================

test_that("plot.causalmixgpd_qte errors on invalid object", {
  skip_if_not_installed("ggplot2")
  invalid <- list(x = 1)
  class(invalid) <- "causalmixgpd_qte"
  expect_error(plot_causalmixgpd_qte(invalid), "Invalid QTE object")
})

test_that("plot.causalmixgpd_ate errors on invalid object", {
  skip_if_not_installed("ggplot2")
  invalid <- list(x = 1)
  class(invalid) <- "causalmixgpd_ate"
  expect_error(plot_causalmixgpd_ate(invalid), "Invalid ATE object")
})

# ======================================================================
# Additional S3 method tests for bundles and fits
# ======================================================================

pkg <- "CausalMixGPD"
print_causalmixgpd_bundle <- getFromNamespace("print.causalmixgpd_bundle", pkg)
print_causalmixgpd_causal_bundle <- getFromNamespace("print.causalmixgpd_causal_bundle", pkg)
summary_causalmixgpd_bundle <- getFromNamespace("summary.causalmixgpd_bundle", pkg)
summary_causalmixgpd_causal_bundle <- getFromNamespace("summary.causalmixgpd_causal_bundle", pkg)
print_causalmixgpd_ps_bundle <- getFromNamespace("print.causalmixgpd_ps_bundle", pkg)
summary_causalmixgpd_ps_bundle <- getFromNamespace("summary.causalmixgpd_ps_bundle", pkg)
print_causalmixgpd_causal_fit <- getFromNamespace("print.causalmixgpd_causal_fit", pkg)
summary_causalmixgpd_causal_fit <- getFromNamespace("summary.causalmixgpd_causal_fit", pkg)
print_causalmixgpd_ps_fit <- getFromNamespace("print.causalmixgpd_ps_fit", pkg)
summary_causalmixgpd_ps_fit <- getFromNamespace("summary.causalmixgpd_ps_fit", pkg)
print_mixgpd_fit <- getFromNamespace("print.mixgpd_fit", pkg)
summary_mixgpd_fit <- getFromNamespace("summary.mixgpd_fit", pkg)
print_mixgpd_summary <- getFromNamespace("print.mixgpd_summary", pkg)

# ======================================================================
# print/summary.causalmixgpd_bundle tests
# ======================================================================

test_that("print.causalmixgpd_bundle works with basic bundle", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )
  expect_output(print_causalmixgpd_bundle(bundle), "CausalMixGPD bundle")
  expect_output(print_causalmixgpd_bundle(bundle), "Stick-Breaking")
  expect_output(print_causalmixgpd_bundle(bundle), "Normal")
})

test_that("print.causalmixgpd_bundle shows code when requested", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )
  expect_output(print_causalmixgpd_bundle(bundle, code = TRUE), "Model code")
})

test_that("summary.causalmixgpd_bundle works", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )
  expect_output(summary_causalmixgpd_bundle(bundle), "CausalMixGPD bundle summary")
  expect_output(summary_causalmixgpd_bundle(bundle), "Parameter specification")
  expect_output(summary_causalmixgpd_bundle(bundle), "Monitors")
})

# ======================================================================
# print/summary.causalmixgpd_causal_bundle tests
# ======================================================================

test_that("print.causalmixgpd_causal_bundle works", {
  set.seed(1)
  y <- rnorm(20)
  X <- matrix(rnorm(40), ncol = 2)
  A <- c(rep(0, 10), rep(1, 10))

  bundle <- build_causal_bundle(
    y = y, X = X, A = A,
    backend = "sb", kernel = "normal",
    GPD = FALSE, components = 4
  )
  expect_output(print_causalmixgpd_causal_bundle(bundle), "CausalMixGPD causal bundle")
  expect_output(print_causalmixgpd_causal_bundle(bundle), "Outcome.*treated")
  expect_output(print_causalmixgpd_causal_bundle(bundle), "Outcome.*control")
})

test_that("print.causalmixgpd_causal_bundle shows PS info for observational", {
  set.seed(1)
  y <- rnorm(20)
  X <- matrix(rnorm(40), ncol = 2)
  A <- c(rep(0, 10), rep(1, 10))

  bundle <- build_causal_bundle(
    y = y, X = X, A = A,
    backend = "sb", kernel = "normal",
    GPD = FALSE, components = 4, PS = "logit"
  )
  expect_output(print_causalmixgpd_causal_bundle(bundle), "PS model.*logit")
})

test_that("summary.causalmixgpd_causal_bundle works", {
  set.seed(1)
  y <- rnorm(20)
  X <- matrix(rnorm(40), ncol = 2)
  A <- c(rep(0, 10), rep(1, 10))

  bundle <- build_causal_bundle(
    y = y, X = X, A = A,
    backend = "sb", kernel = "normal",
    GPD = FALSE, components = 4
  )
  expect_output(summary_causalmixgpd_causal_bundle(bundle), "CausalMixGPD causal bundle summary")
})

# ======================================================================
# Mock PS bundle tests
# ======================================================================

make_mock_ps_bundle <- function() {
  out <- list(
    spec = list(
      meta = list(
        type = "ps_logit",
        include_intercept = TRUE
      )
    ),
    code = quote({
      beta[1] ~ dnorm(0, 1)
    }),
    data = list(X = cbind(`(Intercept)` = 1, x1 = c(-1, 0, 1, 2))),
    monitors = "beta"
  )
  class(out) <- "causalmixgpd_ps_bundle"
  out
}

test_that("print.causalmixgpd_ps_bundle works", {
  bundle <- make_mock_ps_bundle()
  expect_output(print_causalmixgpd_ps_bundle(bundle), "PS bundle")
  expect_output(print_causalmixgpd_ps_bundle(bundle), "model: logit")
})

test_that("print.causalmixgpd_ps_bundle shows code when requested", {
  bundle <- make_mock_ps_bundle()
  expect_output(print_causalmixgpd_ps_bundle(bundle, code = TRUE), "beta")
})

test_that("summary.causalmixgpd_ps_bundle works", {
  bundle <- make_mock_ps_bundle()
  expect_output(summary_causalmixgpd_ps_bundle(bundle), "PS bundle")
})

# ======================================================================
# Mock PS fit tests
# ======================================================================

make_mock_ps_fit <- function() {
  out <- list(
    bundle = make_mock_ps_bundle(),
    mcmc = list(
      samples = matrix(
        c(-0.5, 0.2,
          -0.1, 0.4,
          0.2, 0.6,
          0.5, 0.8),
        ncol = 2,
        byrow = TRUE,
        dimnames = list(NULL, c("beta[1]", "beta[2]"))
      )
    )
  )
  class(out) <- "causalmixgpd_ps_fit"
  out
}

test_that("print.causalmixgpd_ps_fit works", {
  fit <- make_mock_ps_fit()
  expect_output(print_causalmixgpd_ps_fit(fit), "CausalMixGPD PS fit")
  expect_output(print_causalmixgpd_ps_fit(fit), "model: logit")
})

test_that("summary.causalmixgpd_ps_fit works", {
  fit <- make_mock_ps_fit()
  summ <- summary_causalmixgpd_ps_fit(fit)
  expect_s3_class(summ, "summary.causalmixgpd_ps_fit")
  expect_true(all(c("parameter", "mean", "sd", "ess") %in% names(summ$table)))
  expect_output(print_summary_causalmixgpd_ps_fit(summ), "CausalMixGPD PS fit summary")
  expect_output(print_summary_causalmixgpd_ps_fit(summ), "Summary table")
  expect_false(any(grepl("\\bess\\b", capture.output(print_summary_causalmixgpd_ps_fit(summ)))))
  expect_true(any(grepl("\\bess\\b", capture.output(print_summary_causalmixgpd_ps_fit(summ, show_ess = TRUE)))))
})

test_that("params.causalmixgpd_ps_fit returns named beta values", {
  fit <- make_mock_ps_fit()
  out <- params(fit)
  expect_s3_class(out, "mixgpd_params")
  expect_true("beta" %in% names(out))
  expect_equal(names(out$beta), c("(Intercept)", "x1"))
  expect_equal(unname(out$beta), c(0.025, 0.5), tolerance = 1e-8)
})

test_that("print.mixgpd_summary hides ess by default", {
  summ <- structure(
    list(
      model = list(
        backend = "sb",
        kernel = "normal",
        gpd = FALSE,
        epsilon = 0.025,
        truncation = list(Kt = 2L),
        n = 20L,
        components = 3L
      ),
      waic = NULL,
      table = data.frame(
        parameter = c("alpha", "weights[1]"),
        mean = c(1.2, 0.7),
        sd = c(0.3, 0.1),
        ess = c(42, 55)
      )
    ),
    class = "mixgpd_summary"
  )

  expect_false(any(grepl("\\bess\\b", capture.output(print.mixgpd_summary(summ)))))
  expect_true(any(grepl("\\bess\\b", capture.output(print.mixgpd_summary(summ, show_ess = TRUE)))))
})

# ======================================================================
# Mock causal fit tests
# ======================================================================

make_mock_mixgpd_fit <- function() {
  out <- list(
    spec = list(
      meta = list(
        backend = "sb",
        kernel = "normal",
        GPD = FALSE,
        N = 50,
        components = 4
      ),
      dispatch = list(backend = "sb")
    ),
    data = list(y = rnorm(50)),
    mcmc = list(niter = 100, nburnin = 20, thin = 1, nchains = 1),
    epsilon = 0.025
  )
  class(out) <- "mixgpd_fit"
  out
}

make_mock_causal_fit <- function() {
  out <- list(
    bundle = list(
      meta = list(
        backend = list(trt = "sb", con = "sb"),
        kernel = list(trt = "normal", con = "normal"),
        GPD = list(trt = FALSE, con = FALSE),
        ps = list(enabled = FALSE, model_type = FALSE)
      ),
      data = list(X = NULL)
    ),
    outcome_fit = list(
      trt = make_mock_mixgpd_fit(),
      con = make_mock_mixgpd_fit()
    ),
    ps_fit = NULL
  )
  class(out) <- "causalmixgpd_causal_fit"
  out
}

test_that("print.causalmixgpd_causal_fit works", {
  fit <- make_mock_causal_fit()
  expect_output(print_causalmixgpd_causal_fit(fit), "CausalMixGPD causal fit")
  expect_output(print_causalmixgpd_causal_fit(fit), "Outcome.*treated")
  expect_output(print_causalmixgpd_causal_fit(fit), "Outcome.*control")
})

test_that("print.causalmixgpd_causal_fit shows PS info", {
  fit <- make_mock_causal_fit()
  expect_output(print_causalmixgpd_causal_fit(fit), "PS model")
})

test_that("summary.causalmixgpd_causal_fit returns arm summaries", {
  fit <- make_mock_causal_fit()
  fit$bundle$data$X <- cbind(x1 = c(-1, 1))
  fit$bundle$meta$ps$enabled <- TRUE
  fit$ps_fit <- make_mock_ps_fit()

  testthat::local_mocked_bindings(
    summary.mixgpd_fit = function(object, pars = NULL, probs = c(0.025, 0.5, 0.975), ...) {
      out <- list(
        model = list(
          backend = "sb",
          kernel = object$spec$meta$kernel,
          gpd = FALSE,
          n = 50,
          components = 4,
          epsilon = 0.025,
          truncation = list(Kt = 2)
        ),
        waic = list(WAIC = 123.4, lppd = -60.1, pWAIC = 1.6),
        table = data.frame(
          parameter = paste0(c("w[1]", "w[2]", "mean[1]", "mean[2]"), "_", object$spec$meta$kernel),
          mean = c(0.6, 0.4, 1.5, 2.5),
          sd = c(0.1, 0.1, 0.2, 0.3)
        )
      )
      class(out) <- "mixgpd_summary"
      out
    },
    .package = "CausalMixGPD"
  )

  summ <- summary_causalmixgpd_causal_fit(fit)
  expect_s3_class(summ, "summary.causalmixgpd_causal_fit")
  expect_s3_class(summ$ps, "summary.causalmixgpd_ps_fit")
  expect_s3_class(summ$outcome$control, "mixgpd_summary")
  expect_s3_class(summ$outcome$treated, "mixgpd_summary")
  expect_output(print_summary_causalmixgpd_causal_fit(summ), "Outcome fits")
  expect_output(print_summary_causalmixgpd_causal_fit(summ), "\\[control\\]")
  expect_output(print_summary_causalmixgpd_causal_fit(summ), "MixGPD summary")
})

test_that("params.causalmixgpd_causal_fit includes PS block when available", {
  fit <- make_mock_causal_fit()
  fit$bundle$data$X <- cbind(x1 = c(-1, 1))
  fit$bundle$meta$ps$enabled <- TRUE
  fit$ps_fit <- make_mock_ps_fit()

  testthat::local_mocked_bindings(
    params.mixgpd_fit = function(object, ...) {
      structure(list(alpha = 1, w = c(0.7, 0.3)), class = "mixgpd_params")
    },
    .package = "CausalMixGPD"
  )

  out <- params(fit)
  expect_s3_class(out, "mixgpd_params_pair")
  expect_true(all(c("ps", "treated", "control") %in% names(out)))
  expect_equal(names(out$ps$beta), c("(Intercept)", "x1"))
  expect_output(print_mixgpd_params_pair(out), "\\[ps\\]")
})

# ======================================================================
# Mock mixgpd_fit print tests
# ======================================================================

test_that("print.mixgpd_fit works", {
  fit <- make_mock_mixgpd_fit()
  expect_output(print_mixgpd_fit(fit), "MixGPD fit")
  expect_output(print_mixgpd_fit(fit), "Stick-Breaking")
  expect_output(print_mixgpd_fit(fit), "Normal")
})

# ======================================================================
# Mock summary output tests
# ======================================================================

make_mock_mixgpd_summary <- function() {
  out <- list(
    model = list(
      backend = "sb",
      kernel = "normal",
      gpd = FALSE,
      n = 50,
      components = 4,
      epsilon = 0.025,
      truncation = list(Kt = 2)
    ),
    waic = list(WAIC = 123.4, lppd = -60.1, pWAIC = 1.6),
    table = data.frame(
      parameter = c("w[1]", "w[2]", "mean[1]", "mean[2]"),
      mean = c(0.6, 0.4, 1.5, 2.5),
      sd = c(0.1, 0.1, 0.2, 0.3)
    )
  )
  class(out) <- "mixgpd_summary"
  out
}

test_that("print.mixgpd_summary works", {
  summ <- make_mock_mixgpd_summary()
  expect_output(print_mixgpd_summary(summ), "MixGPD summary")
  expect_output(print_mixgpd_summary(summ), "Stick-Breaking")
  expect_output(print_mixgpd_summary(summ), "Normal")
})

# ===== END unit/test-s3-methods.R =====

# ===== BEGIN unit/test-theory.R =====
test_that("GPD splice is continuous at the threshold", {
  w <- c(0.3, 0.7)
  mean <- c(-1, 1)
  sd <- c(0.5, 1)
  threshold <- 0.2
  tail_scale <- 1.0
  tail_shape <- 0.2

  p_bulk <- pNormMix(threshold, w = w, mean = mean, sd = sd,
                     lower.tail = TRUE, log.p = FALSE)
  p_splice <- pNormMixGpd(threshold, w = w, mean = mean, sd = sd,
                          threshold = threshold, tail_scale = tail_scale,
                          tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE)
  expect_equal(p_splice, p_bulk, tolerance = 1e-6)
})

test_that("Tail dominates at high quantiles with GPD", {
  w <- c(0.3, 0.7)
  mean <- c(-1, 1)
  sd <- c(0.5, 1)
  threshold <- 0.2
  # Make the tail appreciably heavier so the GPD splice should inflate extremes
  tail_scale <- 2.5
  tail_shape <- 0.5

  q_bulk <- qNormMix(0.99, w = w, mean = mean, sd = sd)
  q_tail <- qNormMixGpd(0.99, w = w, mean = mean, sd = sd,
                        threshold = threshold, tail_scale = tail_scale,
                        tail_shape = tail_shape)
  expect_gt(q_tail, q_bulk)
})

test_that("SB and CRP predictions agree in a small synthetic case", {
  skip_if_not_test_level("ci")  # MCMC tests are slow, skip at cran level

  set.seed(42)
  y <- abs(stats::rnorm(20)) + 0.1
  mcmc_cfg <- list(niter = 40, nburnin = 10, thin = 1, nchains = 1, seed = 1)

  sb_bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 6,
    mcmc = mcmc_cfg
  )
  crp_bundle <- build_nimble_bundle(
    y = y,
    backend = "crp",
    kernel = "normal",
    GPD = FALSE,
    components = 6,
    mcmc = mcmc_cfg
  )

  sb_fit <- run_mcmc_bundle_manual(sb_bundle, show_progress = FALSE)
  crp_fit <- run_mcmc_bundle_manual(crp_bundle, show_progress = FALSE)

  q_sb <- predict(sb_fit, type = "quantile", index = 0.5)$fit[1, 1]
  q_crp <- predict(crp_fit, type = "quantile", index = 0.5)$fit[1, 1]
  expect_lt(abs(q_sb - q_crp), 1.0)
})

test_that("GPD on/off changes high-quantile behavior", {
  skip_if_not_test_level("ci")  # MCMC tests are slow, skip at cran level

  set.seed(123)
  y <- abs(stats::rnorm(20)) + 0.1
  mcmc_cfg <- list(niter = 40, nburnin = 10, thin = 1, nchains = 1, seed = 2)

  fit_off <- run_mcmc_bundle_manual(build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 6,
    mcmc = mcmc_cfg
  ), show_progress = FALSE)

  fit_on <- run_mcmc_bundle_manual(build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = TRUE,
    components = 6,
    mcmc = mcmc_cfg
  ), show_progress = FALSE)

  q_off <- predict(fit_off, type = "quantile", index = 0.99)$fit[1, 1]
  q_on <- predict(fit_on, type = "quantile", index = 0.99)$fit[1, 1]
  pr_on <- params(fit_on)
  ts <- as.numeric(pr_on$tail_scale %||% NA_real_)[1]
  tsh <- as.numeric(pr_on$tail_shape %||% NA_real_)[1]
  if (!is.finite(ts) || !is.finite(tsh) || abs(ts) < 1e-6) {
    skip("Tail parameters not identified in short MCMC; skip GPD effect check.")
  }
  # Turning on the GPD tail should alter high-quantile predictions when tail parameters move.
  expect_gt(abs(q_on - q_off), 0.1)
})

# ===== END unit/test-theory.R =====

# ===== BEGIN unit/test-vectorization.R =====
test_that("lowercase d/p/q functions vectorize over first argument", {

  w <- c(0.60, 0.40)

  shape <- c(2.0, 5.0)
  scale <- c(1.0, 2.0)
  x <- seq(0.5, 2.0, length.out = 5)

  # Use lowercase vectorized wrappers
  d_vec <- dgammamix(x, w = w, shape = shape, scale = scale, log = FALSE)
  d_scalar <- vapply(x, function(xx) {
    as.numeric(dGammaMix(xx, w = w, shape = shape, scale = scale, log = 0L))
  }, numeric(1))
  expect_equal(d_vec, d_scalar, tolerance = 1e-10)

  p_vec <- pgammamix(q = x, w = w, shape = shape, scale = scale,
                     lower.tail = TRUE, log.p = FALSE)
  p_scalar <- vapply(x, function(xx) {
    as.numeric(pGammaMix(xx, w = w, shape = shape, scale = scale,
              lower.tail = 1L, log.p = 0L))
  }, numeric(1))
  expect_equal(p_vec, p_scalar, tolerance = 1e-10)

  probs <- c(0.1, 0.5, 0.9)
  q_vec <- qgammamix(probs, w = w, shape = shape, scale = scale)
  q_scalar <- vapply(probs, function(pp) {
    qGammaMix(pp, w = w, shape = shape, scale = scale)
  }, numeric(1))
  expect_equal(q_vec, q_scalar, tolerance = 1e-10)
})

test_that("lowercase r functions support n > 1", {

  w <- c(0.60, 0.40)
  shape <- c(2.0, 5.0)
  scale <- c(1.0, 2.0)

  set.seed(123)
  draws <- rgammamix(10, w = w, shape = shape, scale = scale)
  expect_length(draws, 10)
  expect_true(all(is.finite(draws)))
})

test_that("lowercase base wrappers vectorize over first argument", {

  x <- c(1.2, 1.5, 2.0)
  probs <- c(0.2, 0.5, 0.9)

  threshold <- 1.0
  scale <- 0.8
  shape <- 0.2

  expect_equal(
    dgpd(x, threshold = threshold, scale = scale, shape = shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dGpd(xx, threshold, scale, shape, 0L)), numeric(1))
  )
  expect_equal(
    pgpd(x, threshold = threshold, scale = scale, shape = shape,
         lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pGpd(xx, threshold, scale, shape, 1L, 0L)), numeric(1))
  )
  expect_equal(
    qgpd(probs, threshold = threshold, scale = scale, shape = shape,
         lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qGpd(pp, threshold, scale, shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  mean <- 2.0
  shape_ig <- 3.0

  expect_equal(
    dinvgauss(x, mean = mean, shape = shape_ig, log = FALSE),
    vapply(x, function(xx) as.numeric(dInvGauss(xx, mean, shape_ig, 0L)), numeric(1))
  )
  expect_equal(
    pinvgauss(x, mean = mean, shape = shape_ig, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pInvGauss(xx, mean, shape_ig, 1L, 0L)), numeric(1))
  )
  expect_equal(
    qinvgauss(probs, mean = mean, shape = shape_ig, lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qInvGauss(pp, mean, shape_ig, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  loc <- 0.0
  scale_am <- 1.2
  shape1 <- 2.0
  shape2 <- 1.5

  expect_equal(
    damoroso(x, loc = loc, scale = scale_am, shape1 = shape1, shape2 = shape2, log = FALSE),
    vapply(x, function(xx) as.numeric(dAmoroso(xx, loc, scale_am, shape1, shape2, 0L)), numeric(1))
  )
  expect_equal(
    pamoroso(x, loc = loc, scale = scale_am, shape1 = shape1, shape2 = shape2,
             lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pAmoroso(xx, loc, scale_am, shape1, shape2, 1L, 0L)),
           numeric(1))
  )
  expect_equal(
    qamoroso(probs, loc = loc, scale = scale_am, shape1 = shape1, shape2 = shape2,
             lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qAmoroso(pp, loc, scale_am, shape1, shape2,
                                        lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  location <- 0.0
  scale_c <- 1.1

  expect_equal(
    dcauchy_vec(x, location = location, scale = scale_c, log = FALSE),
    vapply(x, function(xx) as.numeric(dCauchy(xx, location, scale_c, 0L)), numeric(1))
  )
  expect_equal(
    pcauchy_vec(x, location = location, scale = scale_c, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pCauchy(xx, location, scale_c, 1L, 0L)), numeric(1))
  )
  expect_equal(
    qcauchy_vec(probs, location = location, scale = scale_c, lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qCauchy(pp, location, scale_c, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase base r wrappers support n > 1", {

  set.seed(456)
  draws_gpd <- rgpd(5, threshold = 1.0, scale = 0.8, shape = 0.2)
  expect_length(draws_gpd, 5)
  expect_true(all(is.finite(draws_gpd)))

  draws_ig <- rinvgauss(5, mean = 2.0, shape = 3.0)
  expect_length(draws_ig, 5)
  expect_true(all(is.finite(draws_ig)))

  draws_am <- ramoroso(5, loc = 0.0, scale = 1.2, shape1 = 2.0, shape2 = 1.5)
  expect_length(draws_am, 5)
  expect_true(all(is.finite(draws_am)))

  draws_c <- rcauchy_vec(5, location = 0.0, scale = 1.1)
  expect_length(draws_c, 5)
  expect_true(all(is.finite(draws_c)))
})

test_that("lowercase mix wrappers vectorize over first argument", {

  x <- c(-1.0, 0.5, 1.2)
  probs <- c(0.2, 0.5, 0.9)

  w <- c(0.7, 0.3)
  location <- c(-0.5, 0.8)
  scale <- c(1.0, 1.5)

  expect_equal(
    dcauchymix(x, w = w, location = location, scale = scale, log = FALSE),
    vapply(x, function(xx) as.numeric(dCauchyMix(xx, w = w, location = location, scale = scale, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pcauchymix(x, w = w, location = location, scale = scale, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pCauchyMix(xx, w = w, location = location, scale = scale,
                                                 lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qcauchymix(probs, w = w, location = location, scale = scale, lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qCauchyMix(pp, w = w, location = location, scale = scale,
                                          lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase mixgpd and gpd wrappers vectorize over first argument", {

  x <- c(-0.5, 0.5, 2.5)
  probs <- c(0.1, 0.5, 0.9)

  w <- c(0.6, 0.4)
  mean <- c(0.0, 1.5)
  sd <- c(1.0, 1.8)
  threshold <- 1.2
  tail_scale <- 1.0
  tail_shape <- 0.15

  expect_equal(
    dnormmixgpd(x, w = w, mean = mean, sd = sd, threshold = threshold,
                tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dNormMixGpd(xx, w = w, mean = mean, sd = sd,
                                                 threshold = threshold, tail_scale = tail_scale,
                                                 tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pnormmixgpd(x, w = w, mean = mean, sd = sd, threshold = threshold,
                tail_scale = tail_scale, tail_shape = tail_shape,
                lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pNormMixGpd(xx, w = w, mean = mean, sd = sd,
                                                 threshold = threshold, tail_scale = tail_scale,
                                                 tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qnormmixgpd(probs, w = w, mean = mean, sd = sd, threshold = threshold,
                tail_scale = tail_scale, tail_shape = tail_shape,
                lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qNormMixGpd(pp, w = w, mean = mean, sd = sd,
                                          threshold = threshold, tail_scale = tail_scale,
                                          tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dnormgpd(x, mean = mean[1], sd = sd[1], threshold = threshold,
             tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dNormGpd(xx, mean = mean[1], sd = sd[1],
                                              threshold = threshold, tail_scale = tail_scale,
                                              tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pnormgpd(x, mean = mean[1], sd = sd[1], threshold = threshold,
             tail_scale = tail_scale, tail_shape = tail_shape,
             lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pNormGpd(xx, mean = mean[1], sd = sd[1],
                                              threshold = threshold, tail_scale = tail_scale,
                                              tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qnormgpd(probs, mean = mean[1], sd = sd[1], threshold = threshold,
             tail_scale = tail_scale, tail_shape = tail_shape,
             lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qNormGpd(pp, mean = mean[1], sd = sd[1],
                                        threshold = threshold, tail_scale = tail_scale,
                                        tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase gamma mixgpd and gpd wrappers vectorize over first argument", {
  x <- c(1.0, 2.2, 4.5)
  probs <- c(0.2, 0.6, 0.9)

  w <- c(0.6, 0.4)
  shape <- c(2.0, 4.0)
  scale <- c(1.0, 2.0)
  threshold <- 2.5
  tail_scale <- 0.9
  tail_shape <- 0.2

  expect_equal(
    dgammamixgpd(x, w = w, shape = shape, scale = scale, threshold = threshold,
                 tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dGammaMixGpd(xx, w = w, shape = shape, scale = scale,
                                                  threshold = threshold, tail_scale = tail_scale,
                                                  tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pgammamixgpd(x, w = w, shape = shape, scale = scale, threshold = threshold,
                 tail_scale = tail_scale, tail_shape = tail_shape,
                 lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pGammaMixGpd(xx, w = w, shape = shape, scale = scale,
                                                  threshold = threshold, tail_scale = tail_scale,
                                                  tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qgammamixgpd(probs, w = w, shape = shape, scale = scale, threshold = threshold,
                 tail_scale = tail_scale, tail_shape = tail_shape,
                 lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qGammaMixGpd(pp, w = w, shape = shape, scale = scale,
                                           threshold = threshold, tail_scale = tail_scale,
                                           tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dgammagpd(x, shape = shape[1], scale = scale[1], threshold = threshold,
              tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dGammaGpd(xx, shape = shape[1], scale = scale[1],
                                               threshold = threshold, tail_scale = tail_scale,
                                               tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pgammagpd(x, shape = shape[1], scale = scale[1], threshold = threshold,
              tail_scale = tail_scale, tail_shape = tail_shape,
              lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pGammaGpd(xx, shape = shape[1], scale = scale[1],
                                               threshold = threshold, tail_scale = tail_scale,
                                               tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qgammagpd(probs, shape = shape[1], scale = scale[1], threshold = threshold,
              tail_scale = tail_scale, tail_shape = tail_shape,
              lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qGammaGpd(pp, shape = shape[1], scale = scale[1],
                                        threshold = threshold, tail_scale = tail_scale,
                                        tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase laplace wrappers vectorize over first argument", {
  x <- c(-1.0, 0.5, 2.0)
  probs <- c(0.2, 0.6, 0.9)

  w <- c(0.7, 0.3)
  location <- c(-0.5, 1.0)
  scale <- c(1.0, 1.4)
  threshold <- 0.8
  tail_scale <- 0.9
  tail_shape <- 0.15

  expect_equal(
    dlaplacemix(x, w = w, location = location, scale = scale, log = FALSE),
    vapply(x, function(xx) as.numeric(dLaplaceMix(xx, w = w, location = location, scale = scale, log = 0L)),
           numeric(1))
  )
  expect_equal(
    plaplacemix(x, w = w, location = location, scale = scale, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pLaplaceMix(xx, w = w, location = location, scale = scale,
                                                 lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qlaplacemix(probs, w = w, location = location, scale = scale,
                lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qLaplaceMix(pp, w = w, location = location, scale = scale,
                                          lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dlaplacemixgpd(x, w = w, location = location, scale = scale, threshold = threshold,
                   tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dLaplaceMixGpd(xx, w = w, location = location, scale = scale,
                                                    threshold = threshold, tail_scale = tail_scale,
                                                    tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    plaplacemixgpd(x, w = w, location = location, scale = scale, threshold = threshold,
                   tail_scale = tail_scale, tail_shape = tail_shape,
                   lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pLaplaceMixGpd(xx, w = w, location = location, scale = scale,
                                                    threshold = threshold, tail_scale = tail_scale,
                                                    tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qlaplacemixgpd(probs, w = w, location = location, scale = scale, threshold = threshold,
                   tail_scale = tail_scale, tail_shape = tail_shape,
                   lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qLaplaceMixGpd(pp, w = w, location = location, scale = scale,
                                             threshold = threshold, tail_scale = tail_scale,
                                             tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dlaplacegpd(x, location = location[1], scale = scale[1], threshold = threshold,
               tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dLaplaceGpd(xx, location = location[1], scale = scale[1],
                                                 threshold = threshold, tail_scale = tail_scale,
                                                 tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    plaplacegpd(x, location = location[1], scale = scale[1], threshold = threshold,
               tail_scale = tail_scale, tail_shape = tail_shape,
               lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pLaplaceGpd(xx, location = location[1], scale = scale[1],
                                                 threshold = threshold, tail_scale = tail_scale,
                                                 tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qlaplacegpd(probs, location = location[1], scale = scale[1], threshold = threshold,
               tail_scale = tail_scale, tail_shape = tail_shape,
               lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qLaplaceGpd(pp, location = location[1], scale = scale[1],
                                          threshold = threshold, tail_scale = tail_scale,
                                          tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase lognormal wrappers vectorize over first argument", {
  x <- c(0.5, 1.2, 3.0)
  probs <- c(0.2, 0.6, 0.9)

  w <- c(0.7, 0.3)
  meanlog <- c(0.0, 0.4)
  sdlog <- c(0.3, 0.6)
  threshold <- 1.8
  tail_scale <- 0.8
  tail_shape <- 0.1

  expect_equal(
    dlognormalmix(x, w = w, meanlog = meanlog, sdlog = sdlog, log = FALSE),
    vapply(x, function(xx) as.numeric(dLognormalMix(xx, w = w, meanlog = meanlog, sdlog = sdlog, log = 0L)),
           numeric(1))
  )
  expect_equal(
    plognormalmix(x, w = w, meanlog = meanlog, sdlog = sdlog, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pLognormalMix(xx, w = w, meanlog = meanlog, sdlog = sdlog,
                                                   lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qlognormalmix(probs, w = w, meanlog = meanlog, sdlog = sdlog,
                  lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qLognormalMix(pp, w = w, meanlog = meanlog, sdlog = sdlog,
                                            lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dlognormalmixgpd(x, w = w, meanlog = meanlog, sdlog = sdlog, threshold = threshold,
                     tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dLognormalMixGpd(xx, w = w, meanlog = meanlog, sdlog = sdlog,
                                                      threshold = threshold, tail_scale = tail_scale,
                                                      tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    plognormalmixgpd(x, w = w, meanlog = meanlog, sdlog = sdlog, threshold = threshold,
                     tail_scale = tail_scale, tail_shape = tail_shape,
                     lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pLognormalMixGpd(xx, w = w, meanlog = meanlog, sdlog = sdlog,
                                                      threshold = threshold, tail_scale = tail_scale,
                                                      tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qlognormalmixgpd(probs, w = w, meanlog = meanlog, sdlog = sdlog, threshold = threshold,
                     tail_scale = tail_scale, tail_shape = tail_shape,
                     lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qLognormalMixGpd(pp, w = w, meanlog = meanlog, sdlog = sdlog,
                                               threshold = threshold, tail_scale = tail_scale,
                                               tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dlognormalgpd(x, meanlog = meanlog[1], sdlog = sdlog[1], threshold = threshold,
                  tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dLognormalGpd(xx, meanlog = meanlog[1], sdlog = sdlog[1],
                                                   threshold = threshold, tail_scale = tail_scale,
                                                   tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    plognormalgpd(x, meanlog = meanlog[1], sdlog = sdlog[1], threshold = threshold,
                  tail_scale = tail_scale, tail_shape = tail_shape,
                  lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pLognormalGpd(xx, meanlog = meanlog[1], sdlog = sdlog[1],
                                                   threshold = threshold, tail_scale = tail_scale,
                                                   tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qlognormalgpd(probs, meanlog = meanlog[1], sdlog = sdlog[1], threshold = threshold,
                  tail_scale = tail_scale, tail_shape = tail_shape,
                  lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qLognormalGpd(pp, meanlog = meanlog[1], sdlog = sdlog[1],
                                            threshold = threshold, tail_scale = tail_scale,
                                            tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase invgauss wrappers vectorize over first argument", {
  x <- c(0.8, 1.5, 3.0)
  probs <- c(0.2, 0.6, 0.9)

  w <- c(0.6, 0.4)
  mean <- c(1.2, 2.0)
  shape <- c(2.5, 4.0)
  threshold <- 1.6
  tail_scale <- 0.7
  tail_shape <- 0.15

  expect_equal(
    dinvgaussmix(x, w = w, mean = mean, shape = shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dInvGaussMix(xx, w = w, mean = mean, shape = shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pinvgaussmix(x, w = w, mean = mean, shape = shape, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pInvGaussMix(xx, w = w, mean = mean, shape = shape,
                                                  lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qinvgaussmix(probs, w = w, mean = mean, shape = shape,
                 lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qInvGaussMix(pp, w = w, mean = mean, shape = shape,
                                           lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dinvgaussmixgpd(x, w = w, mean = mean, shape = shape, threshold = threshold,
                    tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dInvGaussMixGpd(xx, w = w, mean = mean, shape = shape,
                                                     threshold = threshold, tail_scale = tail_scale,
                                                     tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pinvgaussmixgpd(x, w = w, mean = mean, shape = shape, threshold = threshold,
                    tail_scale = tail_scale, tail_shape = tail_shape,
                    lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pInvGaussMixGpd(xx, w = w, mean = mean, shape = shape,
                                                     threshold = threshold, tail_scale = tail_scale,
                                                     tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qinvgaussmixgpd(probs, w = w, mean = mean, shape = shape, threshold = threshold,
                    tail_scale = tail_scale, tail_shape = tail_shape,
                    lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qInvGaussMixGpd(pp, w = w, mean = mean, shape = shape,
                                              threshold = threshold, tail_scale = tail_scale,
                                              tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    dinvgaussgpd(x, mean = mean[1], shape = shape[1], threshold = threshold,
                 tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dInvGaussGpd(xx, mean = mean[1], shape = shape[1],
                                                  threshold = threshold, tail_scale = tail_scale,
                                                  tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pinvgaussgpd(x, mean = mean[1], shape = shape[1], threshold = threshold,
                 tail_scale = tail_scale, tail_shape = tail_shape,
                 lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pInvGaussGpd(xx, mean = mean[1], shape = shape[1],
                                                  threshold = threshold, tail_scale = tail_scale,
                                                  tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qinvgaussgpd(probs, mean = mean[1], shape = shape[1], threshold = threshold,
                 tail_scale = tail_scale, tail_shape = tail_shape,
                 lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qInvGaussGpd(pp, mean = mean[1], shape = shape[1],
                                           threshold = threshold, tail_scale = tail_scale,
                                           tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase amoroso wrappers vectorize over first argument", {
  x <- c(0.8, 1.4, 3.0)
  probs <- c(0.2, 0.6, 0.9)

  w <- c(0.5, 0.3, 0.2)
  loc <- c(0.0, 0.5, 1.0)
  scale <- c(1.0, 1.2, 1.5)
  shape1 <- c(2.0, 3.0, 4.0)
  shape2 <- c(1.1, 1.2, 1.3)
  threshold <- 1.6
  tail_scale <- 0.8
  tail_shape <- 0.15

  expect_equal(
    damorosomix(x, w = w, loc = loc, scale = scale, shape1 = shape1, shape2 = shape2, log = FALSE),
    vapply(x, function(xx) as.numeric(dAmorosoMix(xx, w = w, loc = loc, scale = scale,
                                                 shape1 = shape1, shape2 = shape2, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pamorosomix(x, w = w, loc = loc, scale = scale, shape1 = shape1, shape2 = shape2,
                lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pAmorosoMix(xx, w = w, loc = loc, scale = scale,
                                                 shape1 = shape1, shape2 = shape2,
                                                 lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qamorosomix(probs, w = w, loc = loc, scale = scale, shape1 = shape1, shape2 = shape2,
                lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qAmorosoMix(pp, w = w, loc = loc, scale = scale,
                                          shape1 = shape1, shape2 = shape2,
                                          lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    damorosomixgpd(x, w = w, loc = loc, scale = scale, shape1 = shape1, shape2 = shape2,
                   threshold = threshold, tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dAmorosoMixGpd(xx, w = w, loc = loc, scale = scale,
                                                    shape1 = shape1, shape2 = shape2,
                                                    threshold = threshold, tail_scale = tail_scale,
                                                    tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pamorosomixgpd(x, w = w, loc = loc, scale = scale, shape1 = shape1, shape2 = shape2,
                   threshold = threshold, tail_scale = tail_scale, tail_shape = tail_shape,
                   lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pAmorosoMixGpd(xx, w = w, loc = loc, scale = scale,
                                                    shape1 = shape1, shape2 = shape2,
                                                    threshold = threshold, tail_scale = tail_scale,
                                                    tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qamorosomixgpd(probs, w = w, loc = loc, scale = scale, shape1 = shape1, shape2 = shape2,
                   threshold = threshold, tail_scale = tail_scale, tail_shape = tail_shape,
                   lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qAmorosoMixGpd(pp, w = w, loc = loc, scale = scale,
                                             shape1 = shape1, shape2 = shape2,
                                             threshold = threshold, tail_scale = tail_scale,
                                             tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )

  expect_equal(
    damorosogpd(x, loc = loc[1], scale = scale[1], shape1 = shape1[1], shape2 = shape2[1],
                threshold = threshold, tail_scale = tail_scale, tail_shape = tail_shape, log = FALSE),
    vapply(x, function(xx) as.numeric(dAmorosoGpd(xx, loc = loc[1], scale = scale[1],
                                                 shape1 = shape1[1], shape2 = shape2[1],
                                                 threshold = threshold, tail_scale = tail_scale,
                                                 tail_shape = tail_shape, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pamorosogpd(x, loc = loc[1], scale = scale[1], shape1 = shape1[1], shape2 = shape2[1],
                threshold = threshold, tail_scale = tail_scale, tail_shape = tail_shape,
                lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pAmorosoGpd(xx, loc = loc[1], scale = scale[1],
                                                 shape1 = shape1[1], shape2 = shape2[1],
                                                 threshold = threshold, tail_scale = tail_scale,
                                                 tail_shape = tail_shape, lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qamorosogpd(probs, loc = loc[1], scale = scale[1], shape1 = shape1[1], shape2 = shape2[1],
                threshold = threshold, tail_scale = tail_scale, tail_shape = tail_shape,
                lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qAmorosoGpd(pp, loc = loc[1], scale = scale[1],
                                          shape1 = shape1[1], shape2 = shape2[1],
                                          threshold = threshold, tail_scale = tail_scale,
                                          tail_shape = tail_shape, lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase normal mix wrappers vectorize over first argument", {
  x <- c(-1.0, 0.2, 1.5)
  probs <- c(0.2, 0.6, 0.9)

  w <- c(0.7, 0.3)
  mean <- c(0.0, 1.2)
  sd <- c(1.0, 1.6)

  expect_equal(
    dnormmix(x, w = w, mean = mean, sd = sd, log = FALSE),
    vapply(x, function(xx) as.numeric(dNormMix(xx, w = w, mean = mean, sd = sd, log = 0L)),
           numeric(1))
  )
  expect_equal(
    pnormmix(x, w = w, mean = mean, sd = sd, lower.tail = TRUE, log.p = FALSE),
    vapply(x, function(xx) as.numeric(pNormMix(xx, w = w, mean = mean, sd = sd,
                                              lower.tail = 1L, log.p = 0L)),
           numeric(1))
  )
  expect_equal(
    qnormmix(probs, w = w, mean = mean, sd = sd,
             lower.tail = TRUE, log.p = FALSE),
    vapply(probs, function(pp) qNormMix(pp, w = w, mean = mean, sd = sd,
                                       lower.tail = TRUE, log.p = FALSE),
           numeric(1))
  )
})

test_that("lowercase mixgpd and gpd r wrappers support n > 1", {

  w <- c(0.6, 0.4)
  mean <- c(0.0, 1.5)
  sd <- c(1.0, 1.8)
  threshold <- 1.2
  tail_scale <- 1.0
  tail_shape <- 0.15

  set.seed(789)
  draws_mixgpd <- rnormmixgpd(7, w = w, mean = mean, sd = sd, threshold = threshold,
                             tail_scale = tail_scale, tail_shape = tail_shape)
  expect_length(draws_mixgpd, 7)
  expect_true(all(is.finite(draws_mixgpd)))

  draws_gpd <- rnormgpd(7, mean = mean[1], sd = sd[1], threshold = threshold,
                        tail_scale = tail_scale, tail_shape = tail_shape)
  expect_length(draws_gpd, 7)
  expect_true(all(is.finite(draws_gpd)))
})

test_that("vectorized d/p/q options match scalar versions", {
  d_specs <- list(
    list(vec = dgpd, scalar = dGpd, x = c(1.2, 1.6, 2.0),
         args = list(threshold = 1.0, scale = 0.8, shape = 0.2)),
    list(vec = dnormmix, scalar = dNormMix, x = c(-1.0, 0.5, 1.2),
         args = list(w = c(0.7, 0.3), mean = c(0.0, 1.2), sd = c(1.0, 1.6))),
    list(vec = dnormmixgpd, scalar = dNormMixGpd, x = c(-0.5, 0.5, 2.5),
         args = list(w = c(0.6, 0.4), mean = c(0.0, 1.5), sd = c(1.0, 1.8),
                     threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = dnormgpd, scalar = dNormGpd, x = c(-0.5, 0.5, 2.5),
         args = list(mean = 0.0, sd = 1.0, threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = dgammamix, scalar = dGammaMix, x = c(1.0, 2.2, 4.5),
         args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0))),
    list(vec = dgammamixgpd, scalar = dGammaMixGpd, x = c(1.0, 2.2, 4.5),
         args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0),
                     threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = dgammagpd, scalar = dGammaGpd, x = c(1.0, 2.2, 4.5),
         args = list(shape = 2.0, scale = 1.0, threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = dlognormalmix, scalar = dLognormalMix, x = c(0.5, 1.2, 3.0),
         args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6))),
    list(vec = dlognormalmixgpd, scalar = dLognormalMixGpd, x = c(0.5, 1.2, 3.0),
         args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6),
                     threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = dlognormalgpd, scalar = dLognormalGpd, x = c(0.5, 1.2, 3.0),
         args = list(meanlog = 0.0, sdlog = 0.3, threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = dlaplacemix, scalar = dLaplaceMix, x = c(-1.0, 0.5, 2.0),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4))),
    list(vec = dlaplacemixgpd, scalar = dLaplaceMixGpd, x = c(-1.0, 0.5, 2.0),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4),
                     threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = dlaplacegpd, scalar = dLaplaceGpd, x = c(-1.0, 0.5, 2.0),
         args = list(location = -0.5, scale = 1.0, threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = dinvgaussmix, scalar = dInvGaussMix, x = c(0.8, 1.5, 3.0),
         args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0))),
    list(vec = dinvgaussmixgpd, scalar = dInvGaussMixGpd, x = c(0.8, 1.5, 3.0),
         args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0),
                     threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = dinvgaussgpd, scalar = dInvGaussGpd, x = c(0.8, 1.5, 3.0),
         args = list(mean = 1.2, shape = 2.5, threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = dinvgauss, scalar = dInvGauss, x = c(1.2, 1.5, 2.0),
         args = list(mean = 2.0, shape = 3.0)),
    list(vec = damorosomix, scalar = dAmorosoMix, x = c(0.8, 1.4, 3.0),
         args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                     shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3))),
    list(vec = damorosomixgpd, scalar = dAmorosoMixGpd, x = c(0.8, 1.4, 3.0),
         args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                     shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3),
                     threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = damorosogpd, scalar = dAmorosoGpd, x = c(0.8, 1.4, 3.0),
         args = list(loc = 0.0, scale = 1.0, shape1 = 2.0, shape2 = 1.1,
                     threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = damoroso, scalar = dAmoroso, x = c(1.2, 1.5, 2.0),
         args = list(loc = 0.0, scale = 1.2, shape1 = 2.0, shape2 = 1.5)),
    list(vec = dcauchy_vec, scalar = dCauchy, x = c(-1.0, 0.0, 1.0),
         args = list(location = 0.0, scale = 1.1)),
    list(vec = dcauchymix, scalar = dCauchyMix, x = c(-1.0, 0.5, 1.2),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 0.8), scale = c(1.0, 1.5)))
  )

  for (spec in d_specs) {
    vec <- do.call(spec$vec, c(list(spec$x), spec$args, list(log = TRUE)))
    scalar <- vapply(
      spec$x,
      function(xx) as.numeric(do.call(spec$scalar, c(list(xx), spec$args, list(log = 1L)))),
      numeric(1)
    )
    expect_equal(vec, scalar, tolerance = 1e-10)
  }

  p_specs <- list(
    list(vec = pgpd, scalar = pGpd, q = c(1.2, 1.6, 2.0),
         args = list(threshold = 1.0, scale = 0.8, shape = 0.2)),
    list(vec = pnormmix, scalar = pNormMix, q = c(-1.0, 0.5, 1.2),
         args = list(w = c(0.7, 0.3), mean = c(0.0, 1.2), sd = c(1.0, 1.6))),
    list(vec = pnormmixgpd, scalar = pNormMixGpd, q = c(-0.5, 0.5, 2.5),
         args = list(w = c(0.6, 0.4), mean = c(0.0, 1.5), sd = c(1.0, 1.8),
                     threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = pnormgpd, scalar = pNormGpd, q = c(-0.5, 0.5, 2.5),
         args = list(mean = 0.0, sd = 1.0, threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = pgammamix, scalar = pGammaMix, q = c(1.0, 2.2, 4.5),
         args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0))),
    list(vec = pgammamixgpd, scalar = pGammaMixGpd, q = c(1.0, 2.2, 4.5),
         args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0),
                     threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = pgammagpd, scalar = pGammaGpd, q = c(1.0, 2.2, 4.5),
         args = list(shape = 2.0, scale = 1.0, threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = plognormalmix, scalar = pLognormalMix, q = c(0.5, 1.2, 3.0),
         args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6))),
    list(vec = plognormalmixgpd, scalar = pLognormalMixGpd, q = c(0.5, 1.2, 3.0),
         args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6),
                     threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = plognormalgpd, scalar = pLognormalGpd, q = c(0.5, 1.2, 3.0),
         args = list(meanlog = 0.0, sdlog = 0.3, threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = plaplacemix, scalar = pLaplaceMix, q = c(-1.0, 0.5, 2.0),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4))),
    list(vec = plaplacemixgpd, scalar = pLaplaceMixGpd, q = c(-1.0, 0.5, 2.0),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4),
                     threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = plaplacegpd, scalar = pLaplaceGpd, q = c(-1.0, 0.5, 2.0),
         args = list(location = -0.5, scale = 1.0, threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = pinvgaussmix, scalar = pInvGaussMix, q = c(0.8, 1.5, 3.0),
         args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0))),
    list(vec = pinvgaussmixgpd, scalar = pInvGaussMixGpd, q = c(0.8, 1.5, 3.0),
         args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0),
                     threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = pinvgaussgpd, scalar = pInvGaussGpd, q = c(0.8, 1.5, 3.0),
         args = list(mean = 1.2, shape = 2.5, threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = pinvgauss, scalar = pInvGauss, q = c(1.2, 1.5, 2.0),
         args = list(mean = 2.0, shape = 3.0)),
    list(vec = pamorosomix, scalar = pAmorosoMix, q = c(0.8, 1.4, 3.0),
         args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                     shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3))),
    list(vec = pamorosomixgpd, scalar = pAmorosoMixGpd, q = c(0.8, 1.4, 3.0),
         args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                     shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3),
                     threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = pamorosogpd, scalar = pAmorosoGpd, q = c(0.8, 1.4, 3.0),
         args = list(loc = 0.0, scale = 1.0, shape1 = 2.0, shape2 = 1.1,
                     threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = pamoroso, scalar = pAmoroso, q = c(1.2, 1.5, 2.0),
         args = list(loc = 0.0, scale = 1.2, shape1 = 2.0, shape2 = 1.5)),
    list(vec = pcauchy_vec, scalar = pCauchy, q = c(-1.0, 0.0, 1.0),
         args = list(location = 0.0, scale = 1.1)),
    list(vec = pcauchymix, scalar = pCauchyMix, q = c(-1.0, 0.5, 1.2),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 0.8), scale = c(1.0, 1.5)))
  )

  for (spec in p_specs) {
    vec <- do.call(spec$vec, c(list(spec$q), spec$args, list(lower.tail = FALSE, log.p = TRUE)))
    scalar <- vapply(
      spec$q,
      function(xx) as.numeric(do.call(spec$scalar, c(list(xx), spec$args, list(lower.tail = 0L, log.p = 1L)))),
      numeric(1)
    )
    expect_equal(vec, scalar, tolerance = 1e-10)
  }

  q_specs <- list(
    list(vec = qgpd, scalar = qGpd, p = c(0.1, 0.5, 0.9),
         args = list(threshold = 1.0, scale = 0.8, shape = 0.2)),
    list(vec = qnormmix, scalar = qNormMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.7, 0.3), mean = c(0.0, 1.2), sd = c(1.0, 1.6))),
    list(vec = qnormmixgpd, scalar = qNormMixGpd, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.6, 0.4), mean = c(0.0, 1.5), sd = c(1.0, 1.8),
                     threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = qnormgpd, scalar = qNormGpd, p = c(0.2, 0.5, 0.8),
         args = list(mean = 0.0, sd = 1.0, threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = qgammamix, scalar = qGammaMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0))),
    list(vec = qgammamixgpd, scalar = qGammaMixGpd, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0),
                     threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = qgammagpd, scalar = qGammaGpd, p = c(0.2, 0.5, 0.8),
         args = list(shape = 2.0, scale = 1.0, threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = qlognormalmix, scalar = qLognormalMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6))),
    list(vec = qlognormalmixgpd, scalar = qLognormalMixGpd, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6),
                     threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = qlognormalgpd, scalar = qLognormalGpd, p = c(0.2, 0.5, 0.8),
         args = list(meanlog = 0.0, sdlog = 0.3, threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = qlaplacemix, scalar = qLaplaceMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4))),
    list(vec = qlaplacemixgpd, scalar = qLaplaceMixGpd, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4),
                     threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = qlaplacegpd, scalar = qLaplaceGpd, p = c(0.2, 0.5, 0.8),
         args = list(location = -0.5, scale = 1.0, threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = qinvgaussmix, scalar = qInvGaussMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0))),
    list(vec = qinvgaussmixgpd, scalar = qInvGaussMixGpd, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0),
                     threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = qinvgaussgpd, scalar = qInvGaussGpd, p = c(0.2, 0.5, 0.8),
         args = list(mean = 1.2, shape = 2.5, threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = qinvgauss, scalar = qInvGauss, p = c(0.2, 0.5, 0.8),
         args = list(mean = 2.0, shape = 3.0)),
    list(vec = qamorosomix, scalar = qAmorosoMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                     shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3))),
    list(vec = qamorosomixgpd, scalar = qAmorosoMixGpd, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                     shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3),
                     threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = qamorosogpd, scalar = qAmorosoGpd, p = c(0.2, 0.5, 0.8),
         args = list(loc = 0.0, scale = 1.0, shape1 = 2.0, shape2 = 1.1,
                     threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = qamoroso, scalar = qAmoroso, p = c(0.2, 0.5, 0.8),
         args = list(loc = 0.0, scale = 1.2, shape1 = 2.0, shape2 = 1.5)),
    list(vec = qcauchy_vec, scalar = qCauchy, p = c(0.2, 0.5, 0.8),
         args = list(location = 0.0, scale = 1.1)),
    list(vec = qcauchymix, scalar = qCauchyMix, p = c(0.2, 0.5, 0.8),
         args = list(w = c(0.7, 0.3), location = c(-0.5, 0.8), scale = c(1.0, 1.5)))
  )

  for (spec in q_specs) {
    vec <- do.call(spec$vec, c(list(log(spec$p)), spec$args, list(lower.tail = TRUE, log.p = TRUE)))
    scalar <- vapply(
      log(spec$p),
      function(pp) do.call(spec$scalar, c(list(pp), spec$args, list(lower.tail = TRUE, log.p = TRUE))),
      numeric(1)
    )
    expect_equal(vec, scalar, tolerance = 1e-10)
  }
})

test_that("vectorized r functions match scalar draws with same seed", {
  r_specs <- list(
    list(vec = rgpd, scalar = rGpd, args = list(threshold = 1.0, scale = 0.8, shape = 0.2)),
    list(vec = rnormmix, scalar = rNormMix, args = list(w = c(0.7, 0.3), mean = c(0.0, 1.2), sd = c(1.0, 1.6))),
    list(vec = rnormmixgpd, scalar = rNormMixGpd, args = list(w = c(0.6, 0.4), mean = c(0.0, 1.5), sd = c(1.0, 1.8),
                                                               threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = rnormgpd, scalar = rNormGpd, args = list(mean = 0.0, sd = 1.0, threshold = 1.2, tail_scale = 1.0, tail_shape = 0.15)),
    list(vec = rgammamix, scalar = rGammaMix, args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0))),
    list(vec = rgammamixgpd, scalar = rGammaMixGpd, args = list(w = c(0.6, 0.4), shape = c(2.0, 4.0), scale = c(1.0, 2.0),
                                                                 threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = rgammagpd, scalar = rGammaGpd, args = list(shape = 2.0, scale = 1.0, threshold = 2.5, tail_scale = 0.9, tail_shape = 0.2)),
    list(vec = rlognormalmix, scalar = rLognormalMix, args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6))),
    list(vec = rlognormalmixgpd, scalar = rLognormalMixGpd, args = list(w = c(0.7, 0.3), meanlog = c(0.0, 0.4), sdlog = c(0.3, 0.6),
                                                                       threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = rlognormalgpd, scalar = rLognormalGpd, args = list(meanlog = 0.0, sdlog = 0.3, threshold = 1.8, tail_scale = 0.8, tail_shape = 0.1)),
    list(vec = rlaplacemix, scalar = rLaplaceMix, args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4))),
    list(vec = rlaplacemixgpd, scalar = rLaplaceMixGpd, args = list(w = c(0.7, 0.3), location = c(-0.5, 1.0), scale = c(1.0, 1.4),
                                                                     threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = rlaplacegpd, scalar = rLaplaceGpd, args = list(location = -0.5, scale = 1.0, threshold = 0.8, tail_scale = 0.9, tail_shape = 0.15)),
    list(vec = rinvgaussmix, scalar = rInvGaussMix, args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0))),
    list(vec = rinvgaussmixgpd, scalar = rInvGaussMixGpd, args = list(w = c(0.6, 0.4), mean = c(1.2, 2.0), shape = c(2.5, 4.0),
                                                                     threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = rinvgaussgpd, scalar = rInvGaussGpd, args = list(mean = 1.2, shape = 2.5, threshold = 1.6, tail_scale = 0.7, tail_shape = 0.15)),
    list(vec = rinvgauss, scalar = rInvGauss, args = list(mean = 2.0, shape = 3.0)),
    list(vec = ramorosomix, scalar = rAmorosoMix, args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                                                             shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3))),
    list(vec = ramorosomixgpd, scalar = rAmorosoMixGpd, args = list(w = c(0.5, 0.3, 0.2), loc = c(0.0, 0.5, 1.0), scale = c(1.0, 1.2, 1.5),
                                                                   shape1 = c(2.0, 3.0, 4.0), shape2 = c(1.1, 1.2, 1.3),
                                                                   threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = ramorosogpd, scalar = rAmorosoGpd, args = list(loc = 0.0, scale = 1.0, shape1 = 2.0, shape2 = 1.1,
                                                             threshold = 1.6, tail_scale = 0.8, tail_shape = 0.15)),
    list(vec = ramoroso, scalar = rAmoroso, args = list(loc = 0.0, scale = 1.2, shape1 = 2.0, shape2 = 1.5)),
    list(vec = rcauchy_vec, scalar = rCauchy, args = list(location = 0.0, scale = 1.1)),
    list(vec = rcauchymix, scalar = rCauchyMix, args = list(w = c(0.7, 0.3), location = c(-0.5, 0.8), scale = c(1.0, 1.5)))
  )

  for (spec in r_specs) {
    set.seed(123)
    vec <- do.call(spec$vec, c(list(n = 7), spec$args))
    set.seed(123)
    scalar <- vapply(
      seq_len(7),
      function(i) do.call(spec$scalar, c(list(n = 1), spec$args)),
      numeric(1)
    )
    expect_equal(vec, scalar, tolerance = 1e-10)
  }
})

# ===== END unit/test-vectorization.R =====

# ===== BEGIN unit/test-vectorization-contract.R =====
test_that("lowercase r* supports n = 0", {

  w <- c(0.60, 0.40)
  shape <- c(2.0, 5.0)
  scale <- c(1.0, 2.0)

  # Use lowercase vectorized wrapper
  draws <- rgammamix(0, w = w, shape = shape, scale = scale)
  expect_true(is.numeric(draws))
  expect_length(draws, 0)
})

# ===== END unit/test-vectorization-contract.R =====

# ===== BEGIN unit/test-visualization-helpers.R =====
# tests/testthat/test-visualization-helpers.R
# Unit tests for internal visualization helper functions (06-visualization-helpers.R)

# ======================================================================
# .plot_palette tests
# ======================================================================

test_that(".plot_palette returns correct number of colors", {
  expect_length(.plot_palette(3), 3)
  expect_length(.plot_palette(8), 8)
  expect_length(.plot_palette(1), 1)
})

test_that(".plot_palette recycles when n > 8", {
  pal <- .plot_palette(12)
  expect_length(pal, 12)
  # First 8 should repeat in positions 9-12
  expect_equal(pal[9], pal[1])
  expect_equal(pal[10], pal[2])
})

test_that(".plot_palette returns valid hex colors", {
  pal <- .plot_palette(8)
  expect_true(all(grepl("^#[0-9A-Fa-f]{6}$", pal)))
})

test_that(".plot_palette handles NULL input", {
  pal <- .plot_palette(NULL)
  expect_length(pal, 8)  # default length
})

# ======================================================================
# .plot_theme tests
# ======================================================================

test_that(".plot_theme returns a ggplot theme", {
  skip_if_not_installed("ggplot2")
  theme <- .plot_theme()
  expect_s3_class(theme, "theme")
})

# ======================================================================
# .plot_quantile_pred tests
# ======================================================================

test_that(".plot_quantile_pred works with data.frame fit", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(
      index = c(0.25, 0.5, 0.75),
      estimate = c(1, 2, 3),
      lower = c(0.5, 1.5, 2.5),
      upper = c(1.5, 2.5, 3.5)
    )
  )
  p <- .plot_quantile_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_quantile_pred works with id column", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(
      id = rep(1:2, each = 3),
      index = rep(c(0.25, 0.5, 0.75), 2),
      estimate = rnorm(6),
      lower = rnorm(6) - 0.5,
      upper = rnorm(6) + 0.5
    )
  )
  p <- .plot_quantile_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_quantile_pred errors on non-data.frame fit", {
  pred <- list(fit = c(1, 2, 3))
  expect_error(.plot_quantile_pred(pred), "data frame")
})

# ======================================================================
# .plot_sample_pred tests
# ======================================================================

test_that(".plot_sample_pred works with numeric samples", {
  skip_if_not_installed("ggplot2")
  pred <- list(fit = rnorm(100))
  p <- .plot_sample_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_sample_pred errors on non-numeric", {
  pred <- list(fit = data.frame(x = 1:5))
  expect_error(.plot_sample_pred(pred), "numeric vector")
})

# ======================================================================
# .plot_mean_pred tests
# ======================================================================

test_that(".plot_mean_pred works with data.frame fit and draws", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(estimate = 5.0, lower = 4.5, upper = 5.5),
    draws = rnorm(100, mean = 5)
  )
  p <- .plot_mean_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_mean_pred works with data.frame fit without draws", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(estimate = 5.0, lower = 4.5, upper = 5.5),
    draws = NULL
  )
  p <- .plot_mean_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_mean_pred works with numeric fit (old format)", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = c(4.5, 5.0, 5.5),
    draws = rnorm(100, mean = 5)
  )
  p <- .plot_mean_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_mean_pred works without CI bounds", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(estimate = 5.0),
    draws = rnorm(100, mean = 5)
  )
  p <- .plot_mean_pred(pred)
  expect_s3_class(p, "gg")
})

# ======================================================================
# .plot_density_pred tests
# ======================================================================

test_that(".plot_density_pred works with data.frame fit", {
  skip_if_not_installed("ggplot2")
  y_grid <- seq(0, 5, length.out = 20)
  pred <- list(
    fit = data.frame(
      y = y_grid,
      density = dnorm(y_grid, mean = 2.5)
    )
  )
  p <- .plot_density_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_density_pred works with id column", {
  skip_if_not_installed("ggplot2")
  y_grid <- seq(0, 5, length.out = 10)
  pred <- list(
    fit = data.frame(
      id = rep(1:2, each = 10),
      y = rep(y_grid, 2),
      density = c(dnorm(y_grid, 2), dnorm(y_grid, 3))
    )
  )
  p <- .plot_density_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_density_pred works with numeric vector fit", {
  skip_if_not_installed("ggplot2")
  y_grid <- seq(0, 5, length.out = 20)
  pred <- list(
    fit = dnorm(y_grid, mean = 2.5),
    grid = y_grid
  )
  p <- .plot_density_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_density_pred handles grid column name", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(
      grid = seq(0, 5, length.out = 10),
      density = dnorm(seq(0, 5, length.out = 10), mean = 2.5)
    )
  )
  p <- .plot_density_pred(pred)
  expect_s3_class(p, "gg")
})

# ======================================================================
# .plot_survival_pred tests
# ======================================================================

test_that(".plot_survival_pred works with data.frame fit", {
  skip_if_not_installed("ggplot2")
  y_grid <- seq(0, 5, length.out = 20)
  pred <- list(
    fit = data.frame(
      y = y_grid,
      survival = 1 - pnorm(y_grid, mean = 2.5)
    )
  )
  p <- .plot_survival_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_survival_pred works with id column", {
  skip_if_not_installed("ggplot2")
  y_grid <- seq(0, 5, length.out = 10)
  pred <- list(
    fit = data.frame(
      id = rep(1:2, each = 10),
      y = rep(y_grid, 2),
      survival = c(1 - pnorm(y_grid, 2), 1 - pnorm(y_grid, 3))
    )
  )
  p <- .plot_survival_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_survival_pred works with numeric vector fit", {
  skip_if_not_installed("ggplot2")
  y_grid <- seq(0, 5, length.out = 20)
  pred <- list(
    fit = 1 - pnorm(y_grid, mean = 2.5),
    grid = y_grid
  )
  p <- .plot_survival_pred(pred)
  expect_s3_class(p, "gg")
})

# ======================================================================
# .plot_location_pred tests
# ======================================================================

test_that(".plot_location_pred works with single row (no id)", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(
      mean = 5.0, mean_lower = 4.5, mean_upper = 5.5,
      median = 4.9, median_lower = 4.4, median_upper = 5.4
    )
  )
  p <- .plot_location_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_location_pred works with id column", {
  skip_if_not_installed("ggplot2")
  pred <- list(
    fit = data.frame(
      id = 1:3,
      mean = c(5.0, 5.5, 6.0),
      mean_lower = c(4.5, 5.0, 5.5),
      mean_upper = c(5.5, 6.0, 6.5),
      median = c(4.9, 5.4, 5.9),
      median_lower = c(4.4, 4.9, 5.4),
      median_upper = c(5.4, 5.9, 6.4)
    )
  )
  p <- .plot_location_pred(pred)
  expect_s3_class(p, "gg")
})

test_that(".plot_location_pred errors on non-data.frame fit", {
  pred <- list(fit = c(1, 2, 3))
  expect_error(.plot_location_pred(pred), "data frame")
})

# ===== END unit/test-visualization-helpers.R =====

# ===== BEGIN unit/test-formatting-and-wrapper-coverage.R =====

test_that("formatting helpers handle numeric, matrix, and data.frame inputs", {
  expect_equal(fmt3(c(1, 1.2345)), c("1", "1.234"))
  expect_true(grepl("e", fmt3_sci(c(1, 100000), big = 1000)[2], fixed = TRUE))
  expect_equal(fmt3_vec(numeric(0)), "")
  expect_match(fmt3_vec(c(1, 2.5)), "1, 2.5")

  df <- data.frame(a = c(1.2345, 2.3456), b = c("x", "y"), stringsAsFactors = FALSE)
  df_fmt <- format_df3(df)
  expect_equal(df_fmt$a, c("1.234", "2.346"))
  expect_equal(df_fmt$b, c("x", "y"))

  df_sci <- format_df3_sci(data.frame(a = c(1, 1e5)))
  expect_true(grepl("e", df_sci$a[2], fixed = TRUE))

  mat <- matrix(c(1.2345, 200000, 3.4567, 4.5678), nrow = 2, byrow = TRUE)
  colnames(mat) <- c("x", "y")
  mat_fmt <- format_mat3(mat)
  mat_sci <- format_mat3_sci(mat, big = 1000)
  expect_equal(dim(mat_fmt), dim(mat))
  expect_equal(dimnames(mat_fmt), dimnames(mat))
  expect_true(grepl("e", mat_sci[1, 2], fixed = TRUE))
})

test_that("knitr formatting helpers return stable defaults outside knitr", {
  old_opt <- getOption("knitr.in.progress")
  options(knitr.in.progress = FALSE)
  on.exit(options(knitr.in.progress = old_opt), add = TRUE)

  expect_false(.is_knitr_output())
  expect_equal(.kable_fmt(), "markdown")
})

test_that("wrapper helpers normalize treatment and formula inputs", {
  expect_equal(.coerce_treat(factor(c("control", "treated", "control"))), c(0L, 1L, 0L))
  expect_equal(.coerce_treat(c(FALSE, TRUE, TRUE)), c(0L, 1L, 1L))
  expect_equal(.coerce_treat(c("a", "b", "a")), c(0L, 1L, 0L))
  expect_error(.coerce_treat(c(0, 1, NA)), "cannot contain NA")
  expect_true(.treat_arg_supplied(list(treat = quote(trt)), quote(trt)))
  expect_false(.treat_arg_supplied(list(), quote(NULL)))
  expect_false(.treat_arg_supplied(list(treat = quote(NULL)), quote(NULL)))

  dat <- data.frame(
    y = c(1, 2, 3),
    x1 = c(0.1, 0.2, 0.3),
    id = c("a", "b", "c"),
    trt = c(0, 1, 0)
  )
  parsed <- .parse_formula_yX(y ~ x1 + id, data = dat)
  expect_equal(parsed$y, c(1, 2, 3))
  expect_true(is.matrix(parsed$X))
  expect_equal(colnames(parsed$X), "x1")
  expect_false(parsed$is_unconditional)

  mf <- stats::model.frame(y ~ x1, data = dat)
  expect_equal(.extract_treat_from_data(mf, dat, "trt"), c(0L, 1L, 0L))
  expect_error(.extract_treat_from_data(mf, dat, c(0, 1)), "length does not match")
})

test_that("mcmc override helpers split runner controls and update bundles", {
  parsed <- .normalize_mcmc_inputs(list(niter = 50, nburn = 10, show_progress = FALSE, timing = TRUE))
  expect_equal(parsed$overrides$niter, 50)
  expect_equal(parsed$overrides$nburnin, 10)
  expect_false("nburn" %in% names(parsed$overrides))
  expect_false(parsed$runner$show_progress)
  expect_true(parsed$runner$timing)
  expect_error(.normalize_mcmc_inputs(list(10)), "must be named")
  expect_error(.normalize_mcmc_inputs(list(foo = 1)), "Unknown mcmc argument")

  bundle_obj <- structure(list(mcmc = list(niter = 20, nburnin = 5)), class = "causalmixgpd_bundle")
  updated <- .apply_mcmc_overrides(bundle_obj, list(niter = 100, seed = 9))
  expect_equal(updated$mcmc$niter, 100)
  expect_equal(updated$mcmc$seed, 9)
})

test_that("build-run helpers validate inputs and expose structural metadata", {
  y <- c(1, 2, 3)
  X <- data.frame(x1 = c(0.1, 0.2, 0.3), x2 = c(0.4, 0.5, 0.6))
  built_data <- build_data_from_inputs(y = y, X = X, ps = c(0.2, 0.4, 0.6))
  expect_equal(names(built_data), c("y", "X", "ps"))
  expect_true(is.matrix(built_data$X))
  expect_error(build_data_from_inputs(y = y, X = X[1:2, ]), "nrow\\(X\\) == length\\(y\\)")

  spec <- compile_model_spec(
    y = abs(stats::rnorm(10)) + 0.1,
    X = as.matrix(X[rep(1:3, length.out = 10), ]),
    backend = "sb",
    kernel = "normal",
    GPD = TRUE,
    components = 4
  )
  dims <- build_dimensions_from_spec(spec)
  mons <- build_monitors_from_spec(spec, monitor_v = TRUE, monitor_latent = TRUE)
  expect_true(is.list(dims))
  expect_true(all(c("v", "w", "z") %in% names(dims)))
  expect_true(any(startsWith(names(dims), "beta_")))
  expect_true(any(grepl("^w\\[", mons)))
  expect_true(any(grepl("^v\\[", mons)))
})

test_that("compile_model_spec aligns backend default links and validates link-mode priors", {
  y <- abs(stats::rnorm(8)) + 0.1
  X <- cbind(
    x1 = seq(-1, 1, length.out = 8),
    x2 = rep(c(-0.5, 0.5), length.out = 8)
  )

  specs <- lapply(c("sb", "crp", "spliced"), function(be) {
    compile_model_spec(
      y = y,
      X = X,
      backend = be,
      kernel = "normal",
      GPD = TRUE,
      components = 3L
    )
  })
  names(specs) <- c("sb", "crp", "spliced")

  expect_true(all(vapply(specs, function(sp) identical(sp$plan$bulk$mean$mode, "link"), logical(1))))
  expect_true(all(vapply(specs, function(sp) identical(sp$plan$bulk$sd$mode, "dist"), logical(1))))
  expect_true(all(vapply(specs, function(sp) identical(sp$plan$gpd$threshold$mode, "link"), logical(1))))
  expect_true(all(vapply(specs, function(sp) identical(sp$plan$gpd$tail_scale$mode, "link"), logical(1))))
  expect_true(all(vapply(specs, function(sp) identical(sp$plan$gpd$tail_shape$mode, "dist"), logical(1))))
  expect_true(all(vapply(specs, function(sp) identical(sp$plan$bulk$mean$beta_prior$dist, "normal"), logical(1))))
  expect_true(all(vapply(specs, function(sp) identical(sp$plan$gpd$tail_scale$beta_prior$dist, "normal"), logical(1))))
  expect_equal(specs$sb$plan$gpd$threshold$link_dist$dist, "lognormal")
  expect_equal(specs$crp$plan$gpd$threshold$link_dist$dist, "lognormal")
  expect_equal(specs$spliced$plan$gpd$threshold$link_dist$dist, "lognormal")

  expect_error(
    compile_model_spec(
      y = y,
      X = X,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 3L,
      param_specs = list(
        bulk = list(
          sd = list(mode = "link", link = "identity")
        )
      )
    ),
    "positive-support parameter 'sd'"
  )

  expect_error(
    compile_model_spec(
      y = y,
      X = X,
      backend = "sb",
      kernel = "normal",
      GPD = TRUE,
      components = 3L,
      param_specs = list(
        gpd = list(
          tail_scale = list(mode = "link", link = "identity")
        )
      )
    ),
    "positive-support parameter 'tail_scale'"
  )

  expect_error(
    compile_model_spec(
      y = y,
      X = X,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 3L,
      param_specs = list(
        bulk = list(
          mean = list(
            mode = "link",
            beta_prior = list(dist = "lognormal", args = list(meanlog = 0, sdlog = 0.5))
          )
        )
      )
    ),
    "must use dist = 'normal'"
  )

  spliced_link_dist <- compile_model_spec(
    y = y,
    X = X,
    backend = "spliced",
    kernel = "normal",
    GPD = TRUE,
    components = 3L,
    param_specs = list(
      gpd = list(
        threshold = list(mode = "link", link_dist = list(dist = "lognormal"))
      )
    )
  )
  expect_equal(spliced_link_dist$plan$gpd$threshold$link_dist$dist, "lognormal")
  expect_equal(spliced_link_dist$plan$gpd$sdlog_u$dist, "invgamma")

  spliced_link_dims <- build_dimensions_from_spec(spliced_link_dist)
  spliced_link_mons <- build_monitors_from_spec(spliced_link_dist)
  spliced_link_inits <- build_inits_from_spec(spliced_link_dist, y = y)

  expect_true(all(c("beta_threshold", "threshold_i") %in% names(spliced_link_dims)))
  expect_true(any(grepl("^beta_threshold\\[", spliced_link_mons)))
  expect_true(any(grepl("^threshold_i\\[", spliced_link_mons)))
  expect_true("sdlog_u" %in% spliced_link_mons)
  expect_true(all(c("beta_threshold", "threshold_i", "sdlog_u") %in% names(spliced_link_inits)))

  bundle_spliced <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "spliced",
    kernel = "normal",
    GPD = TRUE,
    components = 3L,
    param_specs = list(
      gpd = list(
        threshold = list(mode = "link", link_dist = list(dist = "lognormal"))
      )
    )
  )
  code_spliced <- paste(deparse(.extract_nimble_code(bundle_spliced$code)), collapse = "\n")
  expect_true(any(grepl("^threshold_i\\[", bundle_spliced$monitors)))
  expect_match(code_spliced, "threshold_i\\[i, k\\] ~ dlnorm")
  expect_match(code_spliced, "threshold_i\\[i, z\\[i\\]\\]")
})

test_that(".configure_samplers preserves conjugate beta_threshold updates", {
  set.seed(1)
  y <- exp(stats::rnorm(12))
  X <- cbind(
    x1 = as.numeric(scale(stats::rnorm(12))),
    x2 = as.numeric(scale(stats::rnorm(12)))
  )

  for (be in c("sb", "crp")) {
    bundle <- build_nimble_bundle(
      y = y,
      X = X,
      backend = be,
      kernel = "lognormal",
      GPD = TRUE,
      components = 3L,
      mcmc = list(niter = 60, nburnin = 20, thin = 2, nchains = 1, seed = 1)
    )
    code <- .extract_nimble_code(bundle$code)
    model <- nimble::nimbleModel(
      code = code,
      data = bundle$data,
      constants = bundle$constants,
      inits = bundle$inits,
      dimensions = bundle$dimensions,
      check = FALSE,
      calculate = FALSE
    )
    conf <- nimble::configureMCMC(model, monitors = bundle$monitors, enableWAIC = TRUE)
    conf <- .configure_samplers(
      conf,
      spec = bundle$spec,
      data_info = list(constants = bundle$constants, dimensions = bundle$dimensions),
      z_update_every = 1L
    )

    has_threshold_sampler <- any(vapply(conf$getSamplers(), function(sc) {
      tgt <- sc$target %||% character(0)
      any(grepl("^beta_threshold\\[", tgt))
    }, logical(1)))

    expect_true(has_threshold_sampler, info = be)
  }
})

test_that("check_glue_validity covers validation and continuity branches", {
  skip_if_not_test_level("ci")

  fit <- .load_fixture("fit_gpd_small.rds")
  expect_error(
    check_glue_validity(fit, grid = c(0.1, NA_real_)),
    "'grid' must be finite numeric"
  )
  expect_error(
    check_glue_validity(fit, n_draws = 0L),
    "'n_draws' must be >= 1"
  )

  out <- check_glue_validity(
    fit,
    grid = seq(0.1, 4, length.out = 25),
    n_draws = 3L,
    check_continuity = TRUE
  )
  expect_true(is.list(out$violations))
  expect_true("continuity" %in% names(out$pass))

  expect_error(
    check_glue_validity(fit, x = matrix(1, ncol = 1)),
    "Unconditional model: 'x' not allowed"
  )
})

# ===== END unit/test-formatting-and-wrapper-coverage.R =====

# ===== BEGIN unit/test-causal-and-cluster-helper-coverage.R =====

predict.fake_mixfit <- function(object,
                                newdata = NULL,
                                x = NULL,
                                y = NULL,
                                ps = NULL,
                                id = NULL,
                                type = c("mean", "quantile", "density", "survival", "prob", "fit", "rmean", "location"),
                                index = NULL,
                                interval = NULL,
                                probs = c(0.025, 0.5, 0.975),
                                store_draws = FALSE,
                                nsim_mean = 200L,
                                ...) {
  type <- match.arg(type)
  arm_shift <- if (identical(object$arm, "trt")) 1 else 0
  x_mat <- if (!is.null(newdata)) as.matrix(newdata) else if (!is.null(x)) as.matrix(x) else NULL
  n_pred <- if (!is.null(x_mat)) nrow(x_mat) else if (!is.null(y)) length(y) else 1L
  id_vals <- if (!is.null(id)) id else seq_len(n_pred)

  if (type %in% c("density", "survival")) {
    est <- arm_shift + ifelse(type == "density", 0.2, 0.8) + seq_len(n_pred) / 100
    return(list(
      fit = data.frame(
        id = id_vals,
        estimate = est,
        lower = est - 0.05,
        upper = est + 0.05
      )
    ))
  }

  if (type %in% c("mean", "rmean")) {
    est <- arm_shift + seq_len(n_pred)
    out <- list(
      fit = data.frame(
        id = id_vals,
        estimate = est,
        lower = est - 0.25,
        upper = est + 0.25
      ),
      diagnostics = if (type == "mean") list(mean_method = "analytic") else list()
    )
    if (isTRUE(store_draws)) {
      out$draws <- cbind(est - 0.2, est + 0.2)
      out$draws <- t(out$draws)
    }
    return(out)
  }

  if (type == "fit") {
    vals <- matrix(
      rep(arm_shift + seq_len(n_pred), each = 4L) + rep(c(-0.2, -0.05, 0.05, 0.2), times = n_pred),
      nrow = 4L
    )
    return(list(
      fit = data.frame(id = id_vals, estimate = colMeans(vals)),
      draws = vals
    ))
  }

  if (type == "quantile") {
    p <- as.numeric(index %||% 0.5)
    m <- length(p)
    if (isTRUE(object$no_quantile_draws)) {
      fit <- data.frame(
        id = rep(id_vals, each = m),
        index = rep(p, times = n_pred),
        estimate = rep(arm_shift + seq_len(n_pred), each = m) + rep(p, times = n_pred),
        lower = rep(arm_shift + seq_len(n_pred), each = m) + rep(p, times = n_pred) - 0.1,
        upper = rep(arm_shift + seq_len(n_pred), each = m) + rep(p, times = n_pred) + 0.1
      )
      return(list(fit = fit))
    }

    draws <- array(NA_real_, dim = c(3L, n_pred, m))
    for (j in seq_len(m)) {
      base <- arm_shift + seq_len(n_pred) + p[j]
      draws[, , j] <- rbind(base - 0.1, base, base + 0.1)
    }
    est <- apply(draws, c(2, 3), mean)
    lower <- apply(draws, c(2, 3), min)
    upper <- apply(draws, c(2, 3), max)
    fit <- data.frame(
      id = rep(id_vals, each = m),
      index = rep(p, times = n_pred),
      estimate = as.vector(t(est)),
      lower = as.vector(t(lower)),
      upper = as.vector(t(upper))
    )
    out <- list(fit = fit)
    if (isTRUE(store_draws)) out$draws <- draws
    out
  } else {
    stop("Unsupported fake predict type.", call. = FALSE)
  }
}

.make_fake_causal_fit <- function(has_X = TRUE, ps_enabled = FALSE, ps_fit = FALSE, no_quantile_draws = FALSE) {
  X <- if (isTRUE(has_X)) cbind(x1 = c(-1, 0, 1), x2 = c(0.5, 1.5, -0.5)) else NULL
  bundle <- list(
    data = list(X = X),
    index = list(trt = c(2L, 3L)),
    meta = list(
      has_x = isTRUE(has_X),
      ps = list(enabled = ps_enabled, include_in_outcome = ps_enabled, model_type = if (ps_enabled) "logit" else FALSE),
      ps_scale = "prob",
      ps_summary = "mean",
      ps_clamp = 1e-6,
      backend = list(trt = "sb", con = "sb"),
      kernel = list(trt = "normal", con = "normal"),
      GPD = list(trt = FALSE, con = FALSE)
    )
  )

  design <- NULL
  ps_fit_obj <- NULL
  if (isTRUE(ps_enabled) && isTRUE(has_X)) {
    design <- .build_ps_bundle(
      A = c(0L, 1L, 1L),
      X = X,
      spec = list(model = "logit", prior = list(mean = 0, sd = 1), include_intercept = TRUE),
      mcmc = mcmc_fast(seed = 11L)
    )
    bundle$design <- design
    if (isTRUE(ps_fit)) {
      ps_fit_obj <- list(mcmc = list(samples = matrix(
        c(0.1, 0.2, -0.1, -0.2, 0.3, 0.1),
        nrow = 2,
        byrow = TRUE,
        dimnames = list(NULL, c("beta[1]", "beta[2]", "beta[3]"))
      )))
    }
  }

  fit <- list(
    bundle = bundle,
    ps_fit = ps_fit_obj,
    ps_hat = if (isTRUE(ps_enabled) && !isTRUE(ps_fit) && isTRUE(has_X)) c(0.3, 0.6, 0.8) else NULL,
    outcome_fit = list(
      trt = structure(list(arm = "trt", no_quantile_draws = no_quantile_draws), class = "fake_mixfit"),
      con = structure(list(arm = "con", no_quantile_draws = no_quantile_draws), class = "fake_mixfit")
    )
  )
  class(fit) <- "causalmixgpd_causal_fit"
  fit
}

test_that("cluster formula parsing and override helpers cover validation branches", {
  dat <- data.frame(
    y = c(1, 2, 3),
    x1 = c(0.1, 0.2, 0.3),
    g = factor(c("a", "b", "a")),
    id = c("u1", "u2", "u3")
  )

  parsed <- .cluster_parse_formula(y ~ x1 + g + id, data = dat)
  expect_equal(parsed$y, c(1, 2, 3))
  expect_true(parsed$has_X)
  expect_true(all(c("x1", "gb") %in% colnames(parsed$X)))
  expect_false(any(grepl("^id", colnames(parsed$X))))
  expect_error(.cluster_parse_formula(NULL, dat), "'formula' is required")
  expect_error(.cluster_parse_formula(y ~ x1, NULL), "'data' is required")

  expect_equal(.cluster_validate_type(NULL, "param"), "param")
  expect_error(
    .cluster_validate_type_requirements("weights", has_X = FALSE, components_missing = FALSE),
    "requires covariates"
  )
  expect_error(
    .cluster_validate_type_requirements("both", has_X = TRUE, components_missing = TRUE),
    "explicit 'components'"
  )

  parsed_overrides <- .cluster_split_overrides(
    list(mean = "identity", tail_scale = "exp", concentration = 2),
    bulk_names = c("mean", "sd"),
    gpd_names = c("threshold", "tail_scale", "tail_shape")
  )
  expect_equal(parsed_overrides$bulk$mean, "identity")
  expect_equal(parsed_overrides$gpd$tail_scale, "exp")
  expect_equal(parsed_overrides$concentration, 2)
  parsed_overrides_mixed <- .cluster_split_overrides(
    list(
      bulk = list(mean = "identity"),
      tail_scale = "exp",
      alpha = 3
    ),
    bulk_names = c("mean", "sd"),
    gpd_names = c("threshold", "tail_scale", "tail_shape")
  )
  expect_equal(parsed_overrides_mixed$bulk$mean, "identity")
  expect_equal(parsed_overrides_mixed$gpd$tail_scale, "exp")
  expect_equal(parsed_overrides_mixed$concentration, 3)
  expect_error(
    .cluster_split_overrides(list(bad = 1), bulk_names = "mean", gpd_names = "tail_scale"),
    "Unknown override names"
  )

  expect_equal(.cluster_normalize_link_entry("softplus")$link, "softplus")
  expect_error(.cluster_normalize_link_entry(1), "strings or lists")

  spec <- list(
    plan = list(
      bulk = list(mean = list(mode = "link", link = "identity")),
      gpd = list(tail_scale = list(mode = "dist")),
      concentration = list(mode = "dist", dist = "gamma")
    )
  )
  spec2 <- .cluster_apply_link_overrides(
    spec,
    link = list(
      mean = "softplus",
      tail_scale = list(mode = "link", link = "exp", beta_prior = list(dist = "normal", args = list(mean = 0, sd = 1)))
    ),
    has_X = TRUE
  )
  expect_equal(spec2$plan$bulk$mean$link, "softplus")
  expect_equal(spec2$plan$gpd$tail_scale$link, "exp")
  expect_error(.cluster_apply_link_overrides(spec, link = list(mean = "identity"), has_X = FALSE), "require covariates")

  expect_equal(.cluster_normalize_prior_entry(2, current_mode = "dist")$mode, "fixed")
  expect_equal(.cluster_normalize_prior_entry("gamma", current_mode = "dist")$dist, "gamma")
  expect_equal(
    .cluster_normalize_prior_entry(list(dist = "normal", args = list(mean = 0, sd = 1)), current_mode = "link", param_name = "mean")$mode,
    "link"
  )
  expect_error(.cluster_normalize_prior_entry(TRUE, current_mode = "dist"), "numeric, character, or list")

  spec3 <- .cluster_apply_prior_overrides(
    spec2,
    priors = list(
      bulk = list(mean = list(value = 1.5)),
      gpd = list(tail_scale = list(dist = "normal", args = list(mean = 0, sd = 0.5))),
      concentration = list(value = 3)
    )
  )
  expect_equal(spec3$plan$bulk$mean$value, 1.5)
  expect_equal(spec3$plan$gpd$tail_scale$mode, "link")
  expect_equal(spec3$plan$gpd$tail_scale$beta_prior$dist, "normal")
  expect_equal(spec3$plan$gpd$tail_scale$beta_prior$args$sd, 0.5)
  expect_equal(spec3$plan$concentration$value, 3)
})

test_that("cluster posterior helper utilities cover draw and density logic", {
  skip_if_not_installed("coda")

  expect_equal(.cluster_draw_indices(5, burnin = 1, thin = 2), c(2L, 4L))
  expect_error(.cluster_draw_indices(5, burnin = 5), "too large")
  expect_error(.cluster_draw_indices(5, thin = 0), ">= 1")

  m <- coda::mcmc(matrix(
    c(1, 1, 2, 2, 2, 1, 1, 2),
    nrow = 2,
    byrow = TRUE,
    dimnames = list(NULL, c("z[1]", "z[2]", "z[3]", "z[4]"))
  ))
  ml <- coda::mcmc.list(m)
  draw_mat <- .cluster_samples_to_matrix(ml)
  z <- .cluster_extract_z_from_matrix(draw_mat)
  expect_equal(dim(z), c(2, 4))
  expect_error(.cluster_samples_to_matrix(matrix(1:4, nrow = 2)), "Expected 'mcmc' or 'mcmc.list'")
  expect_error(.cluster_extract_z_from_matrix(matrix(1:4, nrow = 2)), "No z\\[i\\] columns found")

  psm <- compute_psm(z)
  dahl <- dahl_labels(z, psm)
  scores <- .cluster_compute_scores(z, dahl$labels, psm)
  expect_equal(dim(psm), c(4, 4))
  expect_true(all(diag(psm) == 1))
  expect_equal(length(dahl$labels), 4)
  expect_equal(nrow(scores), 4)
  expect_equal(rowSums(scores), rep(1, 4), tolerance = 1e-8)

  expect_equal(.cluster_link_apply(2, "identity"), 2)
  expect_equal(.cluster_link_apply(0, "exp"), 1)
  expect_equal(.cluster_link_apply(exp(2), "log"), 2)
  expect_equal(.cluster_link_apply(2, "power", link_power = 3), 8)
  expect_error(.cluster_link_apply(2, "power", link_power = NA_real_), "Invalid power link exponent")
  expect_error(.cluster_link_apply(2, "weird"), "Unsupported link")

  expect_equal(.cluster_softmax(c(0, 0)), c(0.5, 0.5))
  gd <- .cluster_extract_gating_draw(
    c("eta[1]" = 0.5, "B[1,1]" = 1.0, "B[1,2]" = -1.0),
    K = 2,
    P = 2
  )
  expect_true(is.list(gd))
  gw <- .cluster_gating_weights(gd, x_row = c(1, 0))
  expect_equal(length(gw), 2)
  expect_equal(sum(gw), 1, tolerance = 1e-8)
  expect_null(.cluster_extract_gating_draw(c("foo" = 1), K = 2, P = 1))
  expect_null(.cluster_gating_weights(NULL, x_row = 1))

  spec <- list(
    meta = list(kernel = "normal", GPD = FALSE, P = 1, components = 2),
    plan = list(
      bulk = list(
        mean = list(mode = "link", link = "identity"),
        sd = list(mode = "fixed", value = 1)
      ),
      gpd = list()
    ),
    signatures = list(bulk = list(args = c("mean", "sd")))
  )
  dens_fun <- .cluster_resolve_density_fun(spec)
  expect_true(is.function(dens_fun))
  dens <- .cluster_component_density(
    spec = spec,
    draw_row = c("beta_mean[1,1]" = 0),
    k = 1L,
    x_row = 2,
    y_val = 0,
    density_fun = dens_fun
  )
  expect_true(is.finite(dens))
  expect_gt(dens, 0)
})

test_that("causal aggregation helpers cover conditional and marginal branches", {
  expect_null(.causal_ate_draws_matrix(NULL))
  expect_equal(dim(.causal_ate_draws_matrix(c(1, 2, 3))), c(3, 1))
  mat <- matrix(1:6, nrow = 2)
  expect_identical(.causal_ate_draws_matrix(mat), mat)
  expect_error(.causal_ate_draws_matrix(array(1:8, dim = c(2, 2, 2))), "Unexpected mean draw dimensions")

  fit_stub <- list(bundle = list(index = list(trt = c(2L, 4L, 4L))))
  expect_equal(.causal_effect_subset_index(fit_stub, n_pred = 1L, subset = "all"), 1L)
  expect_equal(.causal_effect_subset_index(fit_stub, n_pred = 5L, subset = "treated"), c(2L, 4L))
  expect_error(
    .causal_effect_subset_index(list(bundle = list(index = list(trt = integer(0)))), n_pred = 3L, subset = "treated"),
    "No treated rows"
  )

  summ <- .causal_summarize_draw_cols(matrix(c(1, 2, 3, 2, 3, 4), nrow = 3), compute_interval = TRUE, interval = "credible", level = 0.8)
  expect_equal(length(summ$estimate), 2)
  expect_true(all(is.finite(summ$lower)))
  expect_true(all(is.finite(summ$upper)))

  qobj <- list(
    trt = list(draws = matrix(c(1, 2, 3, 4, 2, 3, 4, 5), nrow = 4, byrow = TRUE)),
    con = list(draws = matrix(c(0, 1, 2, 3, 1, 2, 3, 4), nrow = 4, byrow = TRUE)),
    probs = c(0.25, 0.75),
    level = 0.9,
    interval = "credible",
    n_pred = 2L,
    meta = list(ps_enabled = FALSE)
  )
  qte_out <- .causal_aggregate_qte(qobj, idx = c(1L, 2L), effect_type = "qte")
  expect_s3_class(qte_out, "causalmixgpd_qte")
  expect_equal(nrow(qte_out$fit), 2)
  expect_true(all(c("estimate", "lower", "upper") %in% names(qte_out$fit)))
  expect_s3_class(qte_out$fit_df, "data.frame")

  qobj_uncond <- list(
    trt = list(draws = matrix(c(1, 2, 3, 4), nrow = 2)),
    con = list(draws = matrix(c(0, 1, 1, 2), nrow = 2)),
    probs = 0.5,
    level = 0.95,
    interval = "none",
    n_pred = 1L,
    meta = list(ps_enabled = FALSE)
  )
  qtt_out <- .causal_aggregate_qte(qobj_uncond, idx = 1L, effect_type = "qtt")
  expect_s3_class(qtt_out, "causalmixgpd_qte")
  expect_equal(qtt_out$type, "qtt")

  aobj <- list(
    trt = list(draws = matrix(c(2, 4, 6, 8), nrow = 2)),
    con = list(draws = matrix(c(1, 2, 3, 4), nrow = 2)),
    level = 0.9,
    interval = "credible",
    nsim_mean = 20L,
    meta = list(ps_enabled = FALSE)
  )
  ate_out <- .causal_aggregate_ate(aobj, idx = c(1L, 2L), effect_type = "ate")
  expect_s3_class(ate_out, "causalmixgpd_ate")
  expect_length(ate_out$fit, 1)
  expect_true(is.finite(ate_out$fit))
  expect_s3_class(ate_out$fit_df, "data.frame")
  expect_error(.causal_aggregate_ate(list(trt = list(draws = NULL), con = list(draws = NULL)), idx = 1L), "requires treated/control mean draws")
})

test_that("causal validation and bundle helpers cover conditional and input branches", {
  fit_cond <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = FALSE)
  fit_uncond <- .make_fake_causal_fit(has_X = FALSE, ps_enabled = FALSE)

  expect_true(.causal_is_conditional_model(fit_cond))
  expect_false(.causal_is_conditional_model(fit_uncond))
  expect_silent(.causal_require_conditional(fit_cond, fn = "cate"))
  expect_error(.causal_require_conditional(fit_uncond, fn = "cate"), "cate\\(\\) is available only")
  expect_error(.causal_require_conditional(fit_uncond, fn = "cqte"), "cqte\\(\\) is available only")

  iv_none <- .causal_validate_interval(interval = "none", level = 0.9)
  expect_false(iv_none$compute_interval)
  expect_equal(iv_none$interval, "credible")
  expect_error(.causal_validate_interval(interval = "credible", level = 1), "'level' must be a numeric value between 0 and 1")

  expect_warning(.causal_warn_ignored_marginal_inputs("ate", newdata = matrix(1, ncol = 1), y = 1, conditional_fn = "cate"), "ignored")
  expect_false(.causal_warn_ignored_marginal_inputs("ate", conditional_fn = "cate"))

  expect_equal(.causal_resolve_conditional_x(fit_cond, fn = "cate"), fit_cond$bundle$data$X)
  expect_error(.causal_resolve_conditional_x(fit_uncond, fn = "cate"), "requires covariates")

  x <- cbind(x1 = c(-1, 0, 1, 2), x2 = c(0.5, 1.5, -0.5, 2.5))
  y <- abs(c(1.2, 0.8, 1.5, 2.2))
  A <- c(0L, 1L, 0L, 1L)
  cb <- build_causal_bundle(
    y = y,
    X = x,
    A = A,
    backend = c("sb", "crp"),
    kernel = c("normal", "gamma"),
    GPD = c(FALSE, TRUE),
    components = c(3, 4),
    PS = FALSE,
    monitor = "full",
    mcmc_outcome = mcmc_fast(seed = 5L)
  )
  expect_s3_class(cb, "causalmixgpd_causal_bundle")
  expect_true(cb$meta$monitor_policy$monitor_latent)
  expect_true(cb$meta$monitor_policy$monitor_v)
  expect_false(cb$meta$ps$enabled)

  cb_no_x <- build_causal_bundle(
    y = y,
    X = NULL,
    A = A,
    backend = "sb",
    kernel = "normal",
    components = 3,
    PS = TRUE,
    mcmc_outcome = mcmc_fast(seed = 6L)
  )
  expect_false(cb_no_x$meta$ps$enabled)

  expect_error(build_causal_bundle(y = numeric(0), X = x, A = A, backend = "sb", kernel = "normal", components = 3), "non-empty numeric vector")
  expect_error(build_causal_bundle(y = y, X = x[1:3, , drop = FALSE], A = A, backend = "sb", kernel = "normal", components = 3), "same number of rows")
  expect_error(build_causal_bundle(y = y, X = x, A = c(0L, 1L, 2L, 1L), backend = "sb", kernel = "normal", components = 3), "binary")
  expect_error(build_causal_bundle(y = y, X = x, A = rep(1L, 4), backend = "sb", kernel = "normal", components = 3), "Both treatment arms")
  expect_error(build_causal_bundle(y = y, X = x, A = A, backend = "sb", kernel = "normal", components = c(1, 3)), "treated")
  expect_error(build_causal_bundle(y = y, X = x, A = A, backend = "sb", kernel = "normal", components = c(3, 1)), "control")
  expect_error(build_causal_bundle(y = y, X = x, A = A, backend = "sb", kernel = "normal", components = 3, epsilon = c(1, 0.2)), "treated")
})

test_that("propensity score helper utilities cover bundle construction and prediction", {
  skip_if_not_installed("nimble")

  A <- c(0L, 1L, 0L, 1L)
  X <- cbind(x1 = c(-1, 0, 1, 2), x2 = c(0.5, 1.5, -0.5, 2.5))

  ps_logit <- .build_ps_bundle(
    A = A,
    X = X,
    spec = list(model = "logit", prior = list(mean = 0, sd = 1), include_intercept = TRUE),
    mcmc = mcmc_fast(seed = 1L)
  )
  expect_s3_class(ps_logit, "causalmixgpd_ps_bundle")
  expect_equal(ps_logit$spec$meta$type, "ps_logit")
  design <- .ps_design_matrix(ps_logit, X)
  expect_equal(dim(design), c(4, 3))
  expect_error(.ps_design_matrix(ps_logit, X[, "x1", drop = FALSE]), "do not match PS design")

  ps_probit <- .build_ps_bundle(
    A = A,
    X = X,
    spec = list(model = "probit", prior = list(mean = 0, sd = 1), include_intercept = FALSE),
    mcmc = mcmc_fast(seed = 2L)
  )
  expect_equal(ps_probit$spec$meta$type, "ps_probit")

  ps_naive <- .build_ps_bundle(
    A = A,
    X = X,
    spec = list(model = "naive", prior = list(mean = 0, sd = 1), include_intercept = FALSE),
    mcmc = mcmc_fast(seed = 3L)
  )
  expect_equal(ps_naive$spec$meta$type, "ps_naive")
  expect_error(
    .build_ps_bundle(A = A, X = X, spec = list(model = "bad", prior = list(mean = 0, sd = 1)), mcmc = mcmc_fast(seed = 4L)),
    "Unsupported PS model"
  )

  expect_equal(.apply_ps_scale(c(0.2, 0.8), scale = "prob"), c(0.2, 0.8))
  expect_equal(round(.apply_ps_scale(0.5, scale = "logit"), 6), 0)

  logit_fit <- list(mcmc = list(samples = coda::mcmc(matrix(
    c(0, 1, -1, 0.5, 0.5, -0.5),
    nrow = 2,
    byrow = TRUE,
    dimnames = list(NULL, c("beta[1]", "beta[2]", "beta[3]"))
  ))))
  ps_hat <- .compute_ps_from_fit(logit_fit, ps_logit, X, summary = "mean", clamp = 1e-6)
  expect_equal(length(ps_hat), nrow(X))
  expect_true(all(ps_hat > 0 & ps_hat < 1))

  probit_fit <- list(mcmc = list(samples = coda::mcmc(matrix(
    c(0, 0.5, -0.5, 0.25, 0.25, -0.25),
    nrow = 3,
    byrow = TRUE,
    dimnames = list(NULL, c("beta[1]", "beta[2]"))
  ))))
  probit_hat <- .compute_ps_from_fit(probit_fit, ps_probit, X, summary = "median", clamp = 1e-6)
  expect_equal(length(probit_hat), nrow(X))

  naive_fit <- list(mcmc = list(samples = coda::mcmc(matrix(
    c(
      0.5,
      -0.5, 0.5, 0.0, 1.0,
      1.0, 1.0, 1.0, 1.0
    ),
    nrow = 1,
    dimnames = list(NULL, c(
      "pi_prior",
      "mu[1,1]", "mu[2,1]", "mu[1,2]", "mu[2,2]",
      "sigma[1,1]", "sigma[2,1]", "sigma[1,2]", "sigma[2,2]"
    ))
  ))))
  naive_hat <- .compute_ps_from_fit(naive_fit, ps_naive, X, summary = "mean", clamp = 1e-6)
  expect_equal(length(naive_hat), nrow(X))
  expect_true(all(naive_hat > 0 & naive_hat < 1))
})

test_that("cluster label and psm S3 helpers cover summary and plotting branches", {
  lbl <- structure(
    list(
      labels = c(1L, 1L, 2L),
      components = 2L,
      source = "newdata",
      scores = matrix(c(0.9, 0.1, 0.8, 0.2, 0.25, 0.75), nrow = 3, byrow = TRUE)
    ),
    class = c("dpmixgpd_cluster_labels", "list")
  )
  lbl_sum <- summary(lbl)
  expect_s3_class(lbl_sum, "summary.dpmixgpd_cluster_labels")
  expect_equal(lbl_sum$n, 3)
  expect_silent(plot(lbl, type = "sizes"))
  expect_silent(plot(lbl, type = "certainty"))

  lbl_no_scores <- structure(list(labels = c(1L, 2L), components = 2L), class = c("dpmixgpd_cluster_labels", "list"))
  expect_warning(plot(lbl_no_scores, type = "certainty"), "No score matrix")

  psm_obj <- structure(
    list(psm = matrix(c(1, 0.2, 0.2, 1), 2), components = 2L, draw_index = 1L, psm_max_n = 2L),
    class = c("dpmixgpd_cluster_psm", "list")
  )
  psm_sum <- summary(psm_obj)
  expect_s3_class(psm_sum, "summary.dpmixgpd_cluster_psm")
  expect_equal(psm_sum$diagonal_mean, 1)
  expect_silent(plot(psm_obj, psm_max_n = 2L))
  expect_error(plot(psm_obj, psm_max_n = 1L), "PSM plot blocked")
  expect_error(plot(psm_obj, psm_max_n = 0L), "must be an integer >= 1")
})

test_that("cluster handcrafted fit covers newdata label prediction helpers", {
  skip_if_not_installed("coda")

  train <- data.frame(y = c(0.2, 2.8), x1 = c(0, 2))
  trm <- stats::terms(y ~ x1, data = train)
  mf <- stats::model.frame(trm, data = train)
  mm <- stats::model.matrix(stats::delete.response(trm), data = mf)
  formula_meta <- list(
    terms = trm,
    xlevels = stats::.getXlevels(trm, mf),
    contrasts = attr(mm, "contrasts"),
    X_cols = setdiff(colnames(mm), "(Intercept)"),
    response = "y"
  )

  samp <- coda::mcmc(matrix(
    c(
      1, 2, 0.5, 2.0, 0.0, 3.0,
      1, 2, 0.5, 2.0, 0.0, 3.0
    ),
    nrow = 2,
    byrow = TRUE,
    dimnames = list(NULL, c("z[1]", "z[2]", "eta[1]", "B[1,1]", "beta_mean[1,1]", "beta_mean[2,1]"))
  ))
  fit <- structure(
    list(
      spec = list(
        meta = list(kernel = "normal", GPD = FALSE, P = 1L, components = 2L),
        plan = list(
          bulk = list(
            mean = list(mode = "link", link = "identity"),
            sd = list(mode = "fixed", value = 1)
          ),
          gpd = list()
        ),
        signatures = list(bulk = list(args = c("mean", "sd"))),
        cluster = list(type = "weights", formula_meta = formula_meta)
      ),
      samples = coda::mcmc.list(samp),
      cache_env = new.env(parent = emptyenv())
    ),
    class = c("dpmixgpd_cluster_fit", "list")
  )

  pred <- predict_labels_newdata(fit, newdata = data.frame(y = c(0.1, 2.9), x1 = c(0, 2)))
  expect_equal(length(pred$labels), 2)
  expect_true(is.matrix(pred$scores))
  expect_equal(rowSums(pred$scores), rep(1, 2), tolerance = 1e-8)

  z_draws <- extract_z_draws(fit$samples)
  expect_equal(dim(z_draws), c(2, 2))
})

test_that("causal wrapper functions cover marginal branches with mocked outcome fits", {
  testthat::local_mocked_bindings(
    predict = function(object, ...) {
      if (inherits(object, "fake_mixfit")) {
        return(predict.fake_mixfit(object, ...))
      }
      stats::predict(object, ...)
    },
    .package = "CausalMixGPD"
  )

  fit_cond <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = TRUE, ps_fit = TRUE)
  fit_uncond <- .make_fake_causal_fit(has_X = FALSE, ps_enabled = FALSE)

  qte_out <- NULL
  expect_warning(
    qte_out <- qte(fit_cond, probs = c(0.25, 0.75), newdata = matrix(1, ncol = 1), y = 1, show_progress = FALSE),
    "are ignored"
  )
  expect_s3_class(qte_out, "causalmixgpd_qte")

  qtt_out <- qtt(fit_cond, probs = c(0.25, 0.75), interval = "hpd", show_progress = FALSE)
  expect_s3_class(qtt_out, "causalmixgpd_qte")
  expect_equal(qtt_out$type, "qtt")

  ate_out <- ate(fit_cond, type = "mean", nsim_mean = 10L, show_progress = FALSE)
  att_out <- att(fit_cond, type = "rmean", cutoff = 3, nsim_mean = 10L, show_progress = FALSE)
  expect_s3_class(ate_out, "causalmixgpd_ate")
  expect_s3_class(att_out, "causalmixgpd_ate")
  expect_equal(att_out$type, "att")
  expect_true(is.na(ate_out$nsim_mean))
  expect_equal(att_out$nsim_mean, 10L)
  expect_s3_class(ate_out$fit_df, "data.frame")
  expect_s3_class(att_out$fit_df, "data.frame")
  expect_false(any(grepl("Posterior mean draws", utils::capture.output(print(ate_out)))))
  expect_true(any(grepl("Posterior mean draws", utils::capture.output(print(att_out)))))

  rm_cond <- ate_rmean(fit_cond, cutoff = 2, nsim_mean = 10L, show_progress = FALSE)
  rm_uncond <- ate_rmean(fit_uncond, cutoff = 2, nsim_mean = 10L, show_progress = FALSE)
  expect_s3_class(rm_cond, "causalmixgpd_ate")
  expect_s3_class(rm_uncond, "causalmixgpd_ate")
  expect_s3_class(rm_cond$fit_df, "data.frame")
  expect_s3_class(rm_uncond$fit_df, "data.frame")
})

test_that("conditional causal wrappers cover cate and cqte branches with mocked fits", {
  testthat::local_mocked_bindings(
    predict = function(object, ...) {
      if (inherits(object, "fake_mixfit")) {
        return(predict.fake_mixfit(object, ...))
      }
      stats::predict(object, ...)
    },
    .package = "CausalMixGPD"
  )

  fit_cond <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = TRUE, ps_fit = TRUE)
  fit_no_ps <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = TRUE, ps_fit = FALSE)
  fit_uncond <- .make_fake_causal_fit(has_X = FALSE, ps_enabled = FALSE)

  cate_out <- cate(fit_cond, type = "mean", nsim_mean = 10L, show_progress = FALSE)
  cate_rm <- cate(fit_cond, type = "rmean", cutoff = 2, interval = "hpd", nsim_mean = 10L, show_progress = FALSE)
  cqte_out <- cqte(fit_cond, probs = c(0.25, 0.75), show_progress = FALSE)
  expect_s3_class(cate_out, "causalmixgpd_ate")
  expect_s3_class(cate_rm, "causalmixgpd_ate")
  expect_s3_class(cqte_out, "causalmixgpd_qte")
  expect_equal(length(cate_out$fit), 3)
  expect_equal(dim(cqte_out$fit), c(3, 2))
  expect_s3_class(cate_out$fit_df, "data.frame")
  expect_s3_class(cate_rm$fit_df, "data.frame")
  expect_s3_class(cqte_out$fit_df, "data.frame")
  expect_true(is.na(cate_out$nsim_mean))
  expect_equal(cate_rm$nsim_mean, 10L)

  expect_warning(
    cate(fit_no_ps, type = "mean", nsim_mean = 10L, show_progress = FALSE),
    "missing PS model"
  )
  expect_warning(
    cqte(fit_no_ps, probs = c(0.25, 0.75), show_progress = FALSE),
    "missing PS model"
  )
  expect_error(cate(fit_uncond, show_progress = FALSE), "available only for conditional causal models")
  expect_error(cqte(fit_uncond, show_progress = FALSE), "available only for conditional causal models")
})

test_that("run_mcmc_causal orchestration covers PS, fallback, and validation branches", {
  x <- cbind(x1 = c(-1, 0, 1, 2), x2 = c(0.5, 1.5, -0.5, 2.5))
  y <- abs(c(1.2, 0.8, 1.5, 2.2))
  A <- c(0L, 1L, 0L, 1L)

  bundle_ps <- build_causal_bundle(
    y = y,
    X = x,
    A = A,
    backend = "sb",
    kernel = "normal",
    components = 3,
    PS = "logit",
    mcmc_outcome = mcmc_fast(seed = 12L),
    mcmc_ps = mcmc_fast(seed = 13L)
  )
  bundle_no_ps <- build_causal_bundle(
    y = y,
    X = x,
    A = A,
    backend = "sb",
    kernel = "normal",
    components = 3,
    PS = FALSE,
    mcmc_outcome = mcmc_fast(seed = 14L)
  )

  testthat::local_mocked_bindings(
    .run_ps_mcmc_bundle = function(bundle, ...) {
      list(mcmc = list(samples = matrix(c(0.1, 0.2, -0.1), nrow = 1, dimnames = list(NULL, c("beta[1]", "beta[2]", "beta[3]")))), timing = list(total = 0.1), bundle = bundle)
    },
    .compute_ps_from_fit = function(ps_fit, ps_bundle, X_new, ...) rep(0.4, nrow(as.matrix(X_new))),
    .apply_ps_scale = function(ps_prob, scale = c("prob", "logit"), clamp = 1e-6) as.numeric(ps_prob),
    build_code_from_spec = function(spec) list(code = "fake"),
    build_constants_from_spec = function(spec) list(N = 1L),
    build_monitors_from_spec = function(spec, ...) c("alpha"),
    run_mcmc_bundle_manual = function(bundle, ...) {
      list(samples = NULL, mcmc = list(samples = NULL), timing = list(mcmc = 0.2), bundle = bundle)
    },
    .package = "CausalMixGPD"
  )

  fit_ps <- run_mcmc_causal(bundle_ps, show_progress = FALSE, quiet = TRUE, timing = TRUE)
  fit_no_ps <- run_mcmc_causal(bundle_no_ps, show_progress = FALSE, quiet = TRUE, timing = TRUE)
  expect_s3_class(fit_ps, "causalmixgpd_causal_fit")
  expect_s3_class(fit_no_ps, "causalmixgpd_causal_fit")
  expect_equal(fit_ps$ps_hat, rep(0.4, 4))
  expect_equal(fit_ps$bundle$outcome$con$data$ps, c(0.4, 0.4))
  expect_null(fit_no_ps$ps_fit)
  expect_error(run_mcmc_causal(bundle_ps, z_update_every = 0L, show_progress = FALSE, quiet = TRUE), ">= 1")
})

test_that("cluster wrapper functions cover bundle-fit orchestration with mocked MCMC", {
  fake_cluster_fit <- structure(
    list(samples = NULL, mcmc = list(samples = NULL), timing = list(mcmc = 0.1)),
    class = "fake_cluster_fit"
  )

  testthat::local_mocked_bindings(
    run_cluster_mcmc = function(bundle, ...) {
      structure(
        list(bundle = bundle, spec = bundle$spec, samples = NULL, mcmc = list(samples = NULL), cache_env = new.env(parent = emptyenv())),
        class = c("dpmixgpd_cluster_fit", "list")
      )
    },
    run_mcmc_bundle_manual = function(bundle, ...) fake_cluster_fit,
    .package = "CausalMixGPD"
  )

  dat <- data.frame(y = abs(stats::rnorm(8)) + 0.1, x1 = stats::rnorm(8))
  b_mix <- build_cluster_bundle(y ~ x1, data = dat, kernel = "normal", GPD = FALSE, type = "weights", components = 3, mcmc = mcmc_fast(seed = 15L))
  b_gpd <- build_cluster_bundle(y ~ x1, data = dat, kernel = "normal", GPD = TRUE, type = "weights", components = 3, mcmc = mcmc_fast(seed = 16L))
  expect_s3_class(b_mix, "dpmixgpd_cluster_bundle")
  expect_s3_class(b_gpd, "dpmixgpd_cluster_bundle")

  fit1 <- dpmix.cluster(y ~ x1, data = dat, kernel = "normal", type = "weights", components = 3, mcmc = mcmc_fast(seed = 17L))
  fit2 <- dpmgpd.cluster(y ~ x1, data = dat, kernel = "normal", type = "weights", components = 3, mcmc = mcmc_fast(seed = 18L))
  expect_s3_class(fit1, "dpmixgpd_cluster_fit")
  expect_s3_class(fit2, "dpmixgpd_cluster_fit")
})

test_that("cluster bundle and fit S3 methods cover print summary and plot branches", {
  bundle <- structure(
    list(
      spec = list(
        meta = list(kernel = "normal", GPD = FALSE, N = 4L, P = 2L, components = 3L),
        cluster = list(type = "weights", param_link = TRUE)
      ),
      monitors = c("z[1]", "z[2]")
    ),
    class = c("dpmixgpd_cluster_bundle", "list")
  )

  expect_output(print(bundle), "Cluster bundle")
  sm_bundle <- summary(bundle)
  expect_s3_class(sm_bundle, "summary.dpmixgpd_cluster_bundle")
  expect_equal(sm_bundle$type, "weights")
  expect_silent(plot(bundle))

  fit <- structure(
    list(
      spec = bundle$spec,
      samples = matrix(0, nrow = 3, ncol = 3),
      cache_env = new.env(parent = emptyenv())
    ),
    class = c("dpmixgpd_cluster_fit", "list")
  )

  testthat::local_mocked_bindings(
    predict = function(object, type = c("label", "psm"), ...) {
      type <- match.arg(type)
      if (identical(type, "psm")) {
        return(structure(
          list(
            psm = matrix(c(1, 0.2, 0.2, 1), nrow = 2),
            labels = c(1L, 2L),
            components = 2L,
            draw_index = 1L,
            burnin = 0L,
            thin = 1L,
            psm_max_n = 2000L
          ),
          class = c("dpmixgpd_cluster_psm", "list")
        ))
      }
      structure(
        list(
          labels = c(1L, 1L, 2L),
          components = 2L,
          source = "train",
          burnin = 0L,
          thin = 1L
        ),
        class = c("dpmixgpd_cluster_labels", "list")
      )
    },
    extract_z_draws = function(samples, burnin = NULL, thin = NULL) {
      matrix(c(1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 1L), nrow = 3, byrow = TRUE)
    },
    .package = "CausalMixGPD"
  )

  expect_output(print(fit), "Cluster fit")
  sm_fit <- summary(fit)
  expect_s3_class(sm_fit, "summary.dpmixgpd_cluster_fit")
  expect_equal(sm_fit$K_star, 2L)
  psm_out <- plot(fit, which = "psm", psm_max_n = 50L)
  k_out <- plot(fit, which = "k")
  lbl_out <- plot(fit, which = "sizes")
  expect_s3_class(psm_out, "ggplot")
  expect_s3_class(k_out, "ggplot")
  expect_s3_class(lbl_out, "ggplot")
})

test_that("predict.dpmixgpd_cluster_fit covers label psm newdata and validation", {
  fit <- structure(
    list(
      spec = list(meta = list()),
      samples = "fake_samples",
      cache_env = new.env(parent = emptyenv())
    ),
    class = c("dpmixgpd_cluster_fit", "list")
  )

  testthat::local_mocked_bindings(
    .cluster_samples_to_matrix = function(samples) {
      matrix(
        c(1, 1, 2, 1, 2, 2),
        nrow = 2,
        byrow = TRUE,
        dimnames = list(NULL, c("z[1]", "z[2]", "z[3]"))
      )
    },
    .cluster_draw_indices = function(n_draws, burnin = NULL, thin = NULL) seq_len(n_draws),
    .cluster_extract_z_from_matrix = function(draw_sub) {
      matrix(c(1L, 1L, 2L, 1L, 2L, 2L), nrow = 2, byrow = TRUE)
    },
    compute_psm = function(z_draws) {
      matrix(
        c(1, 0.2, 0.1, 0.2, 1, 0.3, 0.1, 0.3, 1),
        nrow = 3,
        byrow = TRUE
      )
    },
    dahl_labels = function(z_draws, psm) list(labels = c(1L, 2L, 2L), K = 2L, draw_index = 1L),
    .cluster_compute_scores = function(z_draws, labels, psm) {
      matrix(c(0.9, 0.1, 0.2, 0.8, 0.1, 0.9), nrow = 3, byrow = TRUE)
    },
    predict_labels_newdata = function(fit, newdata, burnin = NULL, thin = NULL) {
      list(
        labels = c(1L, 2L),
        scores = matrix(c(0.7, 0.3, 0.4, 0.6), nrow = 2, byrow = TRUE),
        K = 2L
      )
    },
    .package = "CausalMixGPD"
  )

  train_lbl <- predict(fit, type = "label", return_scores = TRUE)
  train_lbl_cached <- predict(fit, type = "label", return_scores = FALSE)
  psm <- predict(fit, type = "psm", psm_max_n = 10L)
  nd_lbl <- predict(fit, newdata = data.frame(y = c(0.1, 0.2)), type = "label", return_scores = TRUE)

  expect_s3_class(train_lbl, "dpmixgpd_cluster_labels")
  expect_s3_class(train_lbl_cached, "dpmixgpd_cluster_labels")
  expect_true(is.matrix(train_lbl$scores))
  expect_s3_class(psm, "dpmixgpd_cluster_psm")
  expect_s3_class(nd_lbl, "dpmixgpd_cluster_labels")
  expect_equal(nd_lbl$source, "newdata")

  expect_error(predict(fit, type = "label", psm_max_n = 0L), "integer >= 1")
  expect_error(
    predict(fit, newdata = data.frame(y = 1), type = "psm"),
    "only available for training data"
  )
  expect_error(
    predict(fit, type = "psm", psm_max_n = 2L),
    "PSM is O\\(n\\^2\\)"
  )
})

test_that("causal prediction wrapper covers quantile, density, survival, prob, and validation branches", {
  testthat::local_mocked_bindings(
    predict = function(object, ...) {
      if (inherits(object, "fake_mixfit")) {
        return(predict.fake_mixfit(object, ...))
      }
      stats::predict(object, ...)
    },
    .package = "CausalMixGPD"
  )

  fit_ps <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = TRUE, ps_fit = TRUE)
  fit_no_ps_model <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = TRUE, ps_fit = FALSE)
  fit_nodraw <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = FALSE, no_quantile_draws = TRUE)

  x_new <- cbind(x1 = c(0, 1), x2 = c(1, 2))
  dens <- predict(fit_ps, newdata = x_new, y = c(0.1, 0.2), type = "density", interval = "credible", show_progress = FALSE)
  surv <- predict(fit_ps, newdata = x_new, y = c(0.1, 0.2), type = "survival", interval = NULL, show_progress = FALSE)
  prob <- predict(fit_ps, newdata = x_new, y = c(0.1, 0.2), type = "prob", interval = NULL, show_progress = FALSE)
  expect_s3_class(dens, "causalmixgpd_causal_predict")
  expect_s3_class(surv, "causalmixgpd_causal_predict")
  expect_s3_class(prob, "causalmixgpd_causal_predict")
  expect_true(all(c("trt_estimate", "con_estimate") %in% names(dens)))

  qdraw <- predict(fit_ps, newdata = x_new, type = "quantile", p = c(0.25, 0.75), show_progress = FALSE)
  qnodraw <- predict(fit_nodraw, newdata = x_new, type = "quantile", p = c(0.25, 0.75), show_progress = FALSE)
  qsingle <- predict(fit_ps, newdata = x_new, type = "quantile", p = 0.5, interval = "none", show_progress = FALSE)
  mean_out <- predict(fit_ps, newdata = x_new, type = "mean", ps = c(0.2, 0.8), show_progress = FALSE)
  expect_s3_class(qdraw, "causalmixgpd_causal_predict")
  expect_s3_class(qnodraw, "causalmixgpd_causal_predict")
  expect_s3_class(qsingle, "causalmixgpd_causal_predict")
  expect_s3_class(mean_out, "causalmixgpd_causal_predict")
  expect_true(all(c("index", "estimate") %in% names(qdraw)))

  expect_error(predict(fit_ps, newdata = x_new, type = "quantile", p = c(0, 0.5), show_progress = FALSE), "must be in \\(0, 1\\)")
  expect_error(predict(fit_ps, newdata = x_new, type = "density", show_progress = FALSE), "requires 'y'")
  expect_error(predict(fit_ps, newdata = x_new, y = 1, type = "density", show_progress = FALSE), "Length of 'y' must match")
  expect_error(predict(fit_ps, newdata = x_new, ps = 0.5, type = "mean", show_progress = FALSE), "Length of 'ps' must equal")
  expect_error(
    predict(fit_no_ps_model, newdata = x_new, type = "mean", show_progress = FALSE),
    "missing PS model"
  )
})

test_that("causal prediction sample branch exposes long-form data frames", {
  testthat::local_mocked_bindings(
    .extract_draws_matrix = function(object, ...) matrix(1:6, nrow = 3L),
    .predict_mixgpd = function(object, x = NULL, id = NULL, type = "sample", nsim = NULL, ...) {
      n_pred <- if (!is.null(x)) nrow(as.matrix(x)) else 1L
      id_vals <- if (!is.null(id) && length(id) == n_pred) id else seq_len(n_pred)
      fit <- if (n_pred == 1L) {
        c(1.1, 1.2, 1.3)
      } else {
        matrix(seq_len(n_pred * 3L) / 10, nrow = n_pred, ncol = 3L)
      }
      list(
        fit = fit,
        fit_df = .values_to_long_df(fit, id = id_vals, value_name = "sample"),
        type = type
      )
    },
    .package = "CausalMixGPD"
  )

  fit_ps <- .make_fake_causal_fit(has_X = TRUE, ps_enabled = TRUE, ps_fit = TRUE)
  x_new <- cbind(x1 = c(0, 1), x2 = c(1, 2))
  out <- predict(fit_ps, newdata = x_new, type = "sample", nsim = 3L, show_progress = FALSE)

  expect_s3_class(out, "causalmixgpd_causal_predict")
  expect_s3_class(out$fit_df, "data.frame")
  expect_s3_class(out$trt_fit_df, "data.frame")
  expect_s3_class(out$con_fit_df, "data.frame")
  expect_true(all(c("id", "draw", "sample") %in% names(out$fit_df)))
})

# ===== END unit/test-causal-and-cluster-helper-coverage.R =====