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

skip_if_not_test_level("ci")

# ===== BEGIN integration/test-bundle.R =====
test_that("Bundle generation: minimal y only", {

  set.seed(1)
  y <- abs(stats::rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 6
  )

  expect_true(inherits(bundle, "causalmixgpd_bundle"), info = "bundle class (y only)")
  expect_true(is.list(bundle$spec$meta), info = "bundle spec meta exists (y only)")
  expect_false(isTRUE(bundle$spec$meta$has_X), info = "has_X is FALSE (y only)")
})

test_that("Bundle generation: with X", {

  set.seed(1)
  y <- abs(stats::rnorm(20)) + 0.1
  X <- data.frame(x1 = stats::rnorm(20), x2 = stats::runif(20))
  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 6
  )

  expect_true(inherits(bundle, "causalmixgpd_bundle"), info = "bundle class (with X)")
  expect_true(isTRUE(bundle$spec$meta$has_X), info = "has_X is TRUE (with X)")
  expect_equal(bundle$spec$meta$P, ncol(X), info = "P matches ncol(X)")
})

test_that("Bundle generation: custom param_specs", {

  set.seed(1)
  y <- abs(stats::rnorm(20)) + 0.1
  X <- data.frame(x1 = stats::rnorm(20), x2 = stats::runif(20))

  param_specs <- list(
    bulk = list(
      mean = list(mode = "link", link = "identity")
    ),
    gpd = list(
      threshold = list(mode = "link", link = "identity"),
      tail_scale = list(mode = "link", link = "exp")
    )
  )

  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "sb",
    kernel = "normal",
    GPD = TRUE,
    components = 6,
    param_specs = param_specs
  )

  plan <- bundle$spec$plan
  expect_equal(plan$bulk$mean$mode, "link", info = "bulk mean link mode")
  expect_equal(plan$gpd$threshold$mode, "link", info = "threshold link mode")
  expect_equal(plan$gpd$tail_scale$mode, "link", info = "tail_scale link mode")
})

test_that("Bundle generation: custom mcmc args", {

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

  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 6,
    mcmc = mcmc_args
  )

  expect_equal(bundle$mcmc$niter, 40, info = "mcmc niter stored")
  expect_equal(bundle$mcmc$nburnin, 10, info = "mcmc nburnin stored")
})

# ===== END integration/test-bundle.R =====

# ===== BEGIN integration/test-bundle-validation.R =====
# tests/testthat/test-bundle-validation.R
# Unit tests for build_nimble_bundle input validation and edge cases

# ======================================================================
# Input validation tests for build_nimble_bundle
# ======================================================================

test_that("build_nimble_bundle errors on empty y", {
  expect_error(
    build_nimble_bundle(
      y = numeric(0),
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 4
    ),
    regexp = "non-empty"
  )
})

test_that("build_nimble_bundle errors on NULL y", {
  expect_error(
    build_nimble_bundle(
      y = NULL,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 4
    ),
    regexp = "non-empty|length"
  )
})

test_that("build_nimble_bundle errors on components < 2", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  expect_error(
    build_nimble_bundle(
      y = y,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 1
    ),
    regexp = ">= 2"
  )
})

test_that("build_nimble_bundle errors on invalid epsilon", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1

  expect_error(
    build_nimble_bundle(
      y = y,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 4,
      epsilon = -0.1
    ),
    regexp = "epsilon|\\[0.*1\\)"
  )

  expect_error(
    build_nimble_bundle(
      y = y,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 4,
      epsilon = 1.5
    ),
    regexp = "epsilon|\\[0.*1\\)"
  )

  expect_error(
    build_nimble_bundle(
      y = y,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 4,
      epsilon = NA
    ),
    regexp = "epsilon|numeric"
  )
})

test_that("build_nimble_bundle errors on ps length mismatch", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  ps <- runif(10)  # wrong length

  expect_error(
    build_nimble_bundle(
      y = y,
      ps = ps,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 4
    ),
    regexp = "same length"
  )
})

# ======================================================================
# Bundle structure tests
# ======================================================================

test_that("build_nimble_bundle returns correct class", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )
  expect_s3_class(bundle, "causalmixgpd_bundle")
})

test_that("build_nimble_bundle has required components", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )

  expect_true("spec" %in% names(bundle))
  expect_true("code" %in% names(bundle))
  expect_true("constants" %in% names(bundle))
  expect_true("data" %in% names(bundle))
  expect_true("inits" %in% names(bundle))
  expect_true("monitors" %in% names(bundle))
  expect_true("mcmc" %in% names(bundle))
  expect_true("epsilon" %in% names(bundle))
})

test_that("build_nimble_bundle stores epsilon value", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4,
    epsilon = 0.05
  )
  expect_equal(bundle$epsilon, 0.05)
})

test_that("build_nimble_bundle stores custom MCMC settings", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4,
    mcmc = list(niter = 100, nburnin = 20, thin = 2, nchains = 2, seed = 42)
  )

  expect_equal(bundle$mcmc$niter, 100)
  expect_equal(bundle$mcmc$nburnin, 20)
  expect_equal(bundle$mcmc$thin, 2)
  expect_equal(bundle$mcmc$nchains, 2)
  expect_equal(bundle$mcmc$seed, 42)
})

# ======================================================================
# Backend-specific tests
# ======================================================================

test_that("build_nimble_bundle works with CRP backend", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "crp",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_equal(bundle$spec$meta$backend, "crp")
})

test_that("build_nimble_bundle works with GPD tail", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = TRUE,
    components = 4
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_true(bundle$spec$meta$GPD)
})

# ======================================================================
# Kernel-specific tests
# ======================================================================

test_that("build_nimble_bundle works with all kernels (sb, no GPD)", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  kernels <- c("normal", "lognormal", "gamma", "invgauss", "laplace", "amoroso", "cauchy")

  for (k in kernels) {
    bundle <- build_nimble_bundle(
      y = y,
      backend = "sb",
      kernel = k,
      GPD = FALSE,
      components = 4
    )
    expect_s3_class(bundle, "causalmixgpd_bundle")
    expect_equal(bundle$spec$meta$kernel, k)
  }
})

test_that("build_nimble_bundle works with X matrix", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  X <- matrix(rnorm(40), ncol = 2)
  colnames(X) <- c("x1", "x2")

  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_true(bundle$spec$meta$has_X)
  expect_equal(bundle$spec$meta$P, 2)
})

test_that("build_nimble_bundle converts data.frame X to matrix", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  X <- data.frame(x1 = rnorm(20), x2 = runif(20))

  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_true(bundle$spec$meta$has_X)
})

# ======================================================================
# alpha_random parameter tests
# ======================================================================

test_that("build_nimble_bundle respects alpha_random = FALSE", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1

  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4,
    alpha_random = FALSE
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
})

test_that("build_nimble_bundle respects alpha_random = TRUE", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1

  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4,
    alpha_random = TRUE
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
})

# ======================================================================
# param_specs override tests
# ======================================================================

test_that("build_nimble_bundle applies param_specs overrides", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  X <- matrix(rnorm(40), ncol = 2)
  colnames(X) <- c("x1", "x2")

  param_specs <- list(
    bulk = list(
      mean = list(mode = "link", link = "identity"),
      sd = list(mode = "dist")
    )
  )

  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4,
    param_specs = param_specs
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_equal(bundle$spec$plan$bulk$mean$mode, "link")
})

# ======================================================================
# J parameter alias tests
# ======================================================================

test_that("build_nimble_bundle uses components parameter correctly", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1

  bundle <- build_nimble_bundle(
    y = y,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 6
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_equal(bundle$spec$meta$components, 6)
})

# ======================================================================
# PS augmentation tests
# ======================================================================

test_that("build_nimble_bundle works with propensity score vector", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  ps <- runif(20)

  bundle <- build_nimble_bundle(
    y = y,
    ps = ps,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  # PS should be included in data
  expect_true("ps" %in% names(bundle$data))
})

test_that("build_nimble_bundle works with X and ps combined", {
  set.seed(1)
  y <- abs(rnorm(20)) + 0.1
  X <- matrix(rnorm(40), ncol = 2)
  ps <- runif(20)

  bundle <- build_nimble_bundle(
    y = y,
    X = X,
    ps = ps,
    backend = "sb",
    kernel = "normal",
    GPD = FALSE,
    components = 4
  )

  expect_s3_class(bundle, "causalmixgpd_bundle")
  expect_true(bundle$spec$meta$has_X)
})

# ===== END integration/test-bundle-validation.R =====

# ===== BEGIN integration/test-fitted.R =====
# test-fitted.R
# Tests for fitted.mixgpd_fit, params(), and plot.mixgpd_fitted
#
# Tier B (ci): Requires MCMC fit

if (!exists(".cache_enabled")) {
  helper_path <- file.path("tests", "testthat", "helper-02-cache.R")
  if (file.exists(helper_path)) source(helper_path)
}

set.seed(123)
N <- 25
y <- abs(rnorm(N)) + 0.1
X <- data.frame(x1 = rnorm(N), x2 = runif(N))

mcmc_cfg <- list(niter = 60, nburnin = 20, thin = 1, nchains = 1, seed = 1)
cache_key <- NULL
if (exists(".cache_enabled") && isTRUE(.cache_enabled())) {
  key_str <- paste("fitted", "conditional", N, mcmc_cfg$niter, mcmc_cfg$nburnin,
                   mcmc_cfg$thin, mcmc_cfg$nchains, mcmc_cfg$seed, sep = "|")
  cache_key <- .cache_hash(key_str)
}
cached <- if (!is.null(cache_key)) .cache_get(cache_key) else NULL

# Only build the fit if we're at ci level or higher
# This prevents expensive MCMC compilation at cran level
fit <- NULL
if (test_level_at_least("ci")) {
  if (!is.null(cached) && inherits(cached$fit, "mixgpd_fit")) {
    fit <- cached$fit
  } else {
    bundle <- build_nimble_bundle(
      y = y,
      X = X,
      backend = "sb",
      kernel = "normal",
      GPD = FALSE,
      components = 5,
      mcmc = mcmc_cfg
    )
    fit <- run_mcmc_bundle_manual(bundle, show_progress = FALSE)
    if (!is.null(cache_key)) .cache_set(cache_key, list(fit = fit))
  }
}

test_that("Conditional fitted returns one value per observation", {
  skip_if_not_test_level("ci")

  ftd <- fitted(fit, type = "median", level = 0.9, seed = 1)
  expect_s3_class(ftd, "mixgpd_fitted")
  expect_equal(nrow(ftd), N)
  expect_true(all(is.finite(ftd$fit)))
})

test_that("Conditional fitted supports quantile type", {
  skip_if_not_test_level("ci")

  ftd <- fitted(fit, type = "quantile", p = 0.8, level = 0.9, seed = 1)
  expect_equal(nrow(ftd), N)
  expect_true(all(is.finite(ftd$fit)))
})

# =============================================================================
# Tests for params() extractor
# =============================================================================
test_that("params() returns mixgpd_params object with expected structure", {
  skip_if_not_test_level("ci")

  p <- params(fit)

  expect_s3_class(p, "mixgpd_params")
  expect_true(is.list(p))

  # Should have alpha (concentration parameter)
  expect_true("alpha" %in% names(p))
  expect_true(is.numeric(p$alpha))
  expect_true(is.finite(p$alpha))

  # Should have weights
  expect_true("w" %in% names(p))
  expect_true(is.numeric(p$w))
  expect_true(all(p$w >= 0))
  expect_true(abs(sum(p$w) - 1) < 0.1)  # Weights should sum to ~1
})

test_that("params() print method works", {
  skip_if_not_test_level("ci")

  p <- params(fit)

  expect_output(print(p), "Posterior mean parameters")
  expect_output(print(p), "alpha")
})

# =============================================================================
# Tests for plot.mixgpd_fitted
# =============================================================================
test_that("plot.mixgpd_fitted returns diagnostic plots", {
  skip_if_not_test_level("ci")
  skip_if_not_installed("ggplot2")

  ftd <- fitted(fit, type = "mean", level = 0.9, seed = 1)

  expect_no_error({
    plots <- plot(ftd)
  })

  expect_s3_class(plots, "mixgpd_fitted_plots")
  expect_true(is.list(plots))
  expect_true("observed_fitted_plot" %in% names(plots))
  expect_true("residual_plot" %in% names(plots))

  # Both should be ggplot objects
  expect_true(inherits(plots$observed_fitted_plot, "ggplot"))
  expect_true(inherits(plots$residual_plot, "ggplot"))
})

# ===== END integration/test-fitted.R =====

# ===== BEGIN integration/test-pit-residuals.R =====
# test-pit-residuals.R
# Focused regression checks for PIT diagnostics and draw filtering.

test_that("PIT residuals include diagnostics and drop invalid draws", {
  skip_if_not_test_level("ci")

  fit <- .load_fixture("fit_cond_small.rds")

  smp <- fit$mcmc$samples %||% fit$samples
  if (is.null(smp) || !is.list(smp)) {
    skip("Fixture has no posterior samples to patch.")
  }

  kernel <- fit$spec$meta$kernel %||% fit$spec$kernel$key
  kdef <- CausalMixGPD::get_kernel_registry()[[kernel]] %||% list()
  bulk_support <- kdef$bulk_support %||% list()
  positive_params <- names(bulk_support)[bulk_support %in%
    c("positive_sd", "positive_scale", "positive_shape", "positive_location")]
  if (!length(positive_params)) {
    skip("No positive-support parameters found to invalidate.")
  }

  ch <- as.matrix(smp[[1]])
  target <- NULL
  for (nm in positive_params) {
    cand <- grep(paste0("^", nm, "\\["), colnames(ch), value = TRUE)
    if (length(cand)) {
      target <- cand[1]
      break
    }
  }
  if (is.null(target)) {
    skip("Could not find a bulk parameter column to invalidate.")
  }

  ch[1L, target] <- -abs(ch[1L, target]) - 1
  smp[[1]] <- coda::as.mcmc(ch)

  fit_patched <- fit
  if (!is.null(fit_patched$mcmc$samples)) fit_patched$mcmc$samples <- smp
  if (!is.null(fit_patched$samples)) fit_patched$samples <- smp

  res_mean <- residuals(fit_patched, type = "pit", pit = "bayes_mean", pit_seed = 1L)
  diag_mean <- attr(res_mean, "pit_diagnostics")
  expect_true(is.list(diag_mean))
  expect_true(diag_mean$n_draws_total >= diag_mean$n_draws_valid)
  expect_true(diag_mean$n_draws_dropped >= 1)
  expect_true(is.numeric(diag_mean$n_draws_used))
  expect_equal(length(diag_mean$n_draws_used), length(res_mean))
  expect_true(all(diag_mean$n_draws_used <= diag_mean$n_draws_total))
  expect_true(all(res_mean >= 0 & res_mean <= 1, na.rm = TRUE))

  res_draw <- residuals(fit_patched, type = "pit", pit = "bayes_draw", pit_seed = 1L)
  diag_draw <- attr(res_draw, "pit_diagnostics")
  expect_true(is.list(diag_draw))
  expect_true(diag_draw$n_draws_total >= diag_draw$n_draws_valid)
  expect_true(diag_draw$n_draws_dropped >= 1)
  expect_true(is.numeric(diag_draw$n_draws_used))
  expect_equal(length(diag_draw$n_draws_used), length(res_draw))
  expect_true(all(diag_draw$n_draws_used <= diag_draw$n_draws_total))
  expect_true(all(res_draw >= 0 & res_draw <= 1, na.rm = TRUE))
})

# ===== END integration/test-pit-residuals.R =====

# ===== BEGIN integration/test-predict-rmean.R =====
# test-predict-rmean.R
# Focused regression checks for mean vs restricted mean under heavy tails.

test_that("predict(type='mean') returns Inf when xi >= 1 and rmean stays finite", {
  skip_if_not_test_level("ci")

  fit <- .load_fixture("fit_gpd_small.rds")

  smp <- fit$mcmc$samples %||% fit$samples
  if (is.null(smp) || !is.list(smp)) {
    skip("Fixture has no posterior samples to patch.")
  }

  ch <- as.matrix(smp[[1]])
  if (!("tail_shape" %in% colnames(ch))) {
    skip("Fixture has no tail_shape draws (non-GPD or unexpected structure).")
  }

  ch[1L, "tail_shape"] <- 1.5
  smp[[1]] <- coda::as.mcmc(ch)

  fit_patched <- fit
  if (!is.null(fit_patched$mcmc$samples)) fit_patched$mcmc$samples <- smp
  if (!is.null(fit_patched$samples)) fit_patched$samples <- smp

  expect_warning(
    pred_mean <- predict(fit_patched, type = "mean", interval = "none", store_draws = FALSE),
    "infinite"
  )
  expect_equal(pred_mean$fit$estimate, Inf)

  pred_rmean <- predict(fit_patched, type = "rmean", cutoff = 50, interval = "none", store_draws = FALSE)
  expect_true(all(is.finite(pred_rmean$fit$estimate)))
})

# ===== END integration/test-predict-rmean.R =====

# ===== BEGIN integration/test-simulated-data.R =====
# Tests for simulated data functions (99-simulated-data.R)

# ============================================================================
# sim_bulk_tail() tests
# ============================================================================

test_that("sim_bulk_tail() returns numeric vector of correct length", {
  result <- sim_bulk_tail(n = 100)

  expect_true(is.numeric(result))
  expect_equal(length(result), 100L)
})
test_that("sim_bulk_tail() returns positive values", {
  result <- sim_bulk_tail(n = 200)

  expect_true(all(result > 0))
})

test_that("sim_bulk_tail() seed reproducibility works", {
  result1 <- sim_bulk_tail(n = 50, seed = 123)
  result2 <- sim_bulk_tail(n = 50, seed = 123)

  expect_equal(result1, result2)
})

test_that("sim_bulk_tail() different seeds produce different results", {
  result1 <- sim_bulk_tail(n = 50, seed = 123)
  result2 <- sim_bulk_tail(n = 50, seed = 456)

  expect_false(identical(result1, result2))
})

test_that("sim_bulk_tail() returns sorted values", {
  result <- sim_bulk_tail(n = 100, seed = 42)

  expect_equal(result, sort(result))
})

test_that("sim_bulk_tail() respects tail_prob parameter", {
  # With higher tail_prob, we expect more extreme values
  result_low <- sim_bulk_tail(n = 100, tail_prob = 0.05, seed = 1)
  result_high <- sim_bulk_tail(n = 100, tail_prob = 0.30, seed = 1)

  # Higher tail_prob should have larger maximum
  # (not always deterministic, but generally true)
  expect_true(is.numeric(result_low))
  expect_true(is.numeric(result_high))
})

# ============================================================================
# sim_causal_qte() tests
# ============================================================================

test_that("sim_causal_qte() returns list with required components", {
  result <- sim_causal_qte(n = 100)

  expect_true(is.list(result))
  expect_true("y" %in% names(result))
  expect_true("t" %in% names(result))
  expect_true("X" %in% names(result))
})

test_that("sim_causal_qte() y is numeric vector of correct length", {
  result <- sim_causal_qte(n = 150)

  expect_true(is.numeric(result$y))
  expect_equal(length(result$y), 150L)
})

test_that("sim_causal_qte() t is binary (0/1)", {
  result <- sim_causal_qte(n = 200, seed = 42)

  expect_true(is.numeric(result$t) || is.integer(result$t))
  expect_true(all(result$t %in% c(0L, 1L)))
  expect_equal(length(result$t), 200L)
})

test_that("sim_causal_qte() X is data frame with 3 columns", {
  result <- sim_causal_qte(n = 100)

  expect_true(is.data.frame(result$X))
  expect_equal(ncol(result$X), 3L)
  expect_equal(nrow(result$X), 100L)
  expect_equal(names(result$X), c("x1", "x2", "x3"))
})

test_that("sim_causal_qte() seed reproducibility works", {
  result1 <- sim_causal_qte(n = 50, seed = 789)
  result2 <- sim_causal_qte(n = 50, seed = 789)

  expect_equal(result1$y, result2$y)
  expect_equal(result1$t, result2$t)
  expect_equal(result1$X, result2$X)
})

test_that("sim_causal_qte() without seed produces different results", {
  set.seed(111)
  result1 <- sim_causal_qte(n = 50)
  set.seed(222)
  result2 <- sim_causal_qte(n = 50)

  expect_false(identical(result1$y, result2$y))
})

# ============================================================================
# sim_survival_tail() tests
# ============================================================================

test_that("sim_survival_tail() returns data frame", {
  result <- sim_survival_tail(n = 100)

  expect_true(is.data.frame(result))
})

test_that("sim_survival_tail() has correct columns", {
  result <- sim_survival_tail(n = 100)

  expect_true("time" %in% names(result))
  expect_true("status" %in% names(result))
  expect_true("x1" %in% names(result))
  expect_true("x2" %in% names(result))
})

test_that("sim_survival_tail() has correct number of rows", {
  result <- sim_survival_tail(n = 150)

  expect_equal(nrow(result), 150L)
})

test_that("sim_survival_tail() time is positive", {
  result <- sim_survival_tail(n = 200, seed = 42)

  expect_true(all(result$time > 0))
})

test_that("sim_survival_tail() status is binary", {
  result <- sim_survival_tail(n = 200, seed = 42)

  expect_true(all(result$status %in% c(0L, 1L)))
})

test_that("sim_survival_tail() x2 is binary", {
  result <- sim_survival_tail(n = 200, seed = 42)

  expect_true(all(result$x2 %in% c(0L, 1L)))
})

test_that("sim_survival_tail() seed reproducibility works", {
  result1 <- sim_survival_tail(n = 50, seed = 321)
  result2 <- sim_survival_tail(n = 50, seed = 321)

  expect_equal(result1, result2)
})

test_that("sim_survival_tail() different seeds produce different results", {
  result1 <- sim_survival_tail(n = 50, seed = 321)
  result2 <- sim_survival_tail(n = 50, seed = 654)

  expect_false(identical(result1$time, result2$time))
})

# ===== END integration/test-simulated-data.R =====

# ===== BEGIN integration/test-cluster-and-causal-coverage.R =====
test_that("cluster workflows cover fit, predict, and S3 methods from cluster.R", {
  skip_if_not_test_level("ci")
  skip_if_not_installed("nimble")
  skip_if(
    identical(Sys.getenv("COVERAGE"), "1") &&
      identical(Sys.getenv("DPMIXGPD_SKIP_COVR_CLUSTER_WORKFLOWS"), "1"),
    "Skipping cluster workflow coverage block under covr"
  )

  set.seed(123)
  dat <- data.frame(
    y = abs(stats::rnorm(24)) + 0.25,
    x1 = stats::rnorm(24),
    x2 = stats::runif(24)
  )

  for (tp in c("weights", "param", "both")) {
    fit <- dpmix.cluster(
      y ~ x1 + x2,
      data = dat,
      kernel = "normal",
      components = 4,
      type = tp,
      mcmc = mcmc_fast(seed = 10L)
    )
    expect_s3_class(fit, "dpmixgpd_cluster_fit")

    psm <- predict(fit, type = "psm")
    expect_s3_class(psm, "dpmixgpd_cluster_psm")
    expect_true(is.matrix(psm$psm))
    expect_equal(nrow(psm$psm), ncol(psm$psm))
    expect_equal(unname(diag(psm$psm)), rep(1, nrow(psm$psm)), tolerance = 1e-8)

    lbl_train <- predict(fit, type = "label")
    expect_s3_class(lbl_train, "dpmixgpd_cluster_labels")
    expect_equal(length(lbl_train$labels), nrow(dat))

    lbl_scores <- predict(fit, type = "label", return_scores = TRUE)
    expect_true(is.matrix(lbl_scores$scores))
    expect_equal(rowSums(lbl_scores$scores), rep(1, nrow(dat)), tolerance = 1e-8)

    nd <- dat[1:6, c("y", "x1", "x2")]
    lbl_new <- predict(fit, newdata = nd, type = "label", return_scores = TRUE)
    expect_equal(length(lbl_new$labels), nrow(nd))
    expect_equal(rowSums(lbl_new$scores), rep(1, nrow(nd)), tolerance = 1e-8)

    expect_output(print(fit), "Cluster fit")
    expect_silent(summary(fit))
    expect_silent(plot(fit, which = "psm"))
    expect_silent(plot(fit, which = "k"))
    expect_silent(plot(fit, which = "sizes"))
    expect_output(print(lbl_scores), "Cluster labels")
    expect_silent(summary(lbl_scores))
    expect_silent(plot(lbl_scores, type = "sizes"))
    expect_silent(plot(lbl_scores, type = "certainty"))
    expect_output(print(psm), "Cluster PSM")
    expect_silent(summary(psm))
    expect_silent(plot(psm, psm_max_n = nrow(psm$psm)))
  }
})

test_that("cluster helpers cover design parsing and override branches", {
  skip_if_not_test_level("ci")
  skip_if(
    identical(Sys.getenv("COVERAGE"), "1") ||
      identical(Sys.getenv("DPMIXGPD_SKIP_COVR_CLUSTER_HELPERS"), "1"),
    "Skipping cluster helper coverage block under covr"
  )

  set.seed(111)
  dat <- data.frame(
    y = abs(stats::rnorm(20)) + 0.2,
    x1 = stats::rnorm(20),
    x2 = stats::runif(20)
  )

  expect_error(
    dpmix.cluster(y ~ 1, data = dat, kernel = "normal", type = "weights", components = 4, mcmc = mcmc_fast(seed = 1L)),
    "requires covariates"
  )
  expect_error(
    dpmix.cluster(y ~ x1 + x2, data = dat, kernel = "normal", type = "weights", mcmc = mcmc_fast(seed = 1L)),
    "explicit 'components'"
  )

  expect_warning(
    fit_param <- dpmix.cluster(y ~ 1, data = dat, kernel = "normal", type = "param", mcmc = mcmc_fast(seed = 2L)),
    "using default components"
  )
  expect_error(
    predict(fit_param, type = "psm", psm_max_n = 10L),
    "PSM is O\\(n\\^2\\)"
  )

  b <- build_cluster_bundle(
    y ~ x1 + x2,
    data = dat,
    kernel = "normal",
    GPD = TRUE,
    type = "both",
    components = 4,
    link = list(
      bulk = list(mean = "identity"),
      gpd = list(tail_scale = list(link = "exp"))
    ),
    priors = list(
      bulk = list(sd = list(dist = "gamma", args = list(shape = 2, rate = 1))),
      gpd = list(
        tail_shape = list(dist = "normal", args = list(mean = 0, sd = 0.25)),
        tail_scale = list(dist = "normal", args = list(mean = 0, sd = 0.7))
      ),
      concentration = list(dist = "gamma", args = list(shape = 3, rate = 2))
    ),
    mcmc = mcmc_fast(seed = 7L)
  )
  expect_equal(b$spec$plan$bulk$mean$link, "identity")
  expect_equal(b$spec$plan$gpd$tail_scale$link, "exp")
  expect_equal(b$spec$plan$concentration$dist, "gamma")

  fn_beta <- getFromNamespace(".cluster_extract_beta_auto", "CausalMixGPD")
  draw_row <- c(
    "beta_tail_scale[1,1]" = 2.0,
    "beta_tail_scale[1,2]" = -1.0,
    "beta_tail_scale[2,1]" = 0.5,
    "beta_tail_scale[2,2]" = 0.25,
    "beta_tail_scale[1]" = 99.0,
    "beta_tail_scale[2]" = 98.0
  )
  expect_equal(as.numeric(fn_beta(draw_row = draw_row, base = "beta_tail_scale", comp = 1L, P = 2L)), c(2.0, -1.0))

  fn_design <- getFromNamespace(".cluster_build_design", "CausalMixGPD")
  train <- data.frame(y = c(1, 2, 3), g = factor(c("a", "b", "a")))
  trm <- stats::terms(y ~ g, data = train)
  mf <- stats::model.frame(trm, data = train)
  mm <- stats::model.matrix(stats::delete.response(trm), data = mf)
  meta <- list(
    terms = trm,
    xlevels = stats::.getXlevels(trm, mf),
    contrasts = attr(mm, "contrasts"),
    X_cols = setdiff(colnames(mm), "(Intercept)"),
    response = "y"
  )
  nd <- data.frame(y = 1, g = factor("c", levels = c("a", "b", "c")))
  expect_error(fn_design(meta = meta, newdata = nd), "unseen factor levels")
})

test_that("causal bundle and fit workflows cover PS and arm-specific branches", {
  skip_if_not_test_level("ci")
  skip_if_not_installed("nimble")
  skip_if(
    identical(Sys.getenv("COVERAGE"), "1") &&
      identical(Sys.getenv("DPMIXGPD_SKIP_COVR_CAUSAL_BRANCHES"), "1"),
    "Skipping causal coverage block under covr"
  )

  sim <- sim_causal_qte(n = 30, seed = 29)
  sim$y <- abs(sim$y) + 0.2

  bundle_ps <- build_causal_bundle(
    y = sim$y,
    X = as.matrix(sim$X),
    A = sim$t,
    backend = c("sb", "crp"),
    kernel = c("normal", "gamma"),
    GPD = c(FALSE, TRUE),
    components = c(3, 4),
    PS = "logit",
    mcmc_outcome = mcmc_fast(seed = 29L),
    mcmc_ps = mcmc_fast(seed = 31L)
  )
  expect_s3_class(bundle_ps, "causalmixgpd_causal_bundle")
  expect_equal(bundle_ps$meta$backend$trt, "sb")
  expect_equal(bundle_ps$meta$backend$con, "crp")

  bundle_no_ps <- build_causal_bundle(
    y = sim$y,
    X = as.matrix(sim$X),
    A = sim$t,
    backend = c("sb", "sb"),
    kernel = c("normal", "normal"),
    GPD = FALSE,
    components = c(3, 3),
    PS = FALSE,
    mcmc_outcome = mcmc_fast(seed = 33L)
  )
  expect_s3_class(bundle_no_ps, "causalmixgpd_causal_bundle")

  fit <- dpmix.causal(
    x = sim$y,
    X = as.matrix(sim$X),
    treat = sim$t,
    backend = c("sb", "crp"),
    kernel = c("normal", "gamma"),
    components = c(3, 3),
    PS = "logit",
    mcmc = c(mcmc_fast(seed = 35L), list(show_progress = FALSE))
  )
  expect_s3_class(fit, "causalmixgpd_causal_fit")
  expect_output(print(fit), "CausalMixGPD causal fit")
  expect_output(print(summary(fit)), "Outcome fits")
  expect_s3_class(params(fit), "mixgpd_params_pair")

  Xnew <- as.matrix(sim$X[1:4, , drop = FALSE])
  pred_mean <- predict(fit, newdata =Xnew, type = "mean", nsim_mean = 20L)
  pred_quant <- predict(fit, newdata =Xnew, type = "quantile", p = c(0.25, 0.75))
  expect_s3_class(pred_mean, "causalmixgpd_causal_predict")
  expect_s3_class(pred_quant, "causalmixgpd_causal_predict")

  expect_s3_class(cate(fit, newdata = Xnew, nsim_mean = 20L, interval = "credible"), "causalmixgpd_ate")
  expect_s3_class(cqte(fit, probs = c(0.25, 0.75), newdata = Xnew, interval = "credible"), "causalmixgpd_qte")
  expect_s3_class(ate(fit, nsim_mean = 20L, interval = "credible"), "causalmixgpd_ate")
  expect_s3_class(qte(fit, probs = c(0.25, 0.75), interval = "credible"), "causalmixgpd_qte")
})

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