# Tests for support generics, generic dist methods, and utility functions.
#
# This file tests:
#   - Support generics (has, infimum, supremum) polymorphic dispatch
#   - interval edge cases not covered in test-interval.R
#   - finite_set edge cases not covered in test-finite_set.R
#   - Generic dist methods: expectation.dist, conditional.dist, rmap.dist,
#     summary.dist, summary_dist, print.summary_dist
#   - Default methods: sampler.default, vcov.default
#   - is_dist predicate
#   - expectation_data utility

# ============================================================================
# Support generics: polymorphic dispatch across interval and finite_set
# ============================================================================

test_that("has generic dispatches correctly to interval and finite_set", {
  # Given: An interval and a finite_set with overlapping values
  iv <- interval$new(lower = 0, upper = 10,
                     lower_closed = TRUE, upper_closed = TRUE)
  fs <- finite_set$new(c(2, 5, 8))

  # When/Then: has() dispatches correctly to each class

  expect_true(has(iv, 5))
  expect_false(has(iv, -1))
  expect_true(has(fs, 5))
  expect_false(has(fs, 4))
})

test_that("infimum generic dispatches correctly to interval and finite_set", {
  iv <- interval$new(lower = -3, upper = 7)
  fs <- finite_set$new(c(10, 20, 30))

  expect_equal(infimum(iv), -3)
  expect_equal(infimum(fs), 10)
})

test_that("supremum generic dispatches correctly to interval and finite_set", {
  iv <- interval$new(lower = -3, upper = 7)
  fs <- finite_set$new(c(10, 20, 30))

  expect_equal(supremum(iv), 7)
  expect_equal(supremum(fs), 30)
})

# ============================================================================
# interval: additional edge cases
# ============================================================================

test_that("interval has() checks scalar membership for 1D interval", {
  # Note: For a 1D interval, has() accepts a scalar x and returns a scalar.
  # The ifelse() in has() uses the scalar lower_closed/upper_closed as the
  # test argument, so it returns a length-1 result.
  iv <- interval$new(lower = 0, upper = 10,
                     lower_closed = TRUE, upper_closed = TRUE)

  expect_true(iv$has(0))
  expect_true(iv$has(5))
  expect_true(iv$has(10))
  expect_false(iv$has(-1))
  expect_false(iv$has(11))
})

test_that("interval has() checks scalar membership for open interval", {
  iv <- interval$new(lower = 0, upper = 10,
                     lower_closed = FALSE, upper_closed = FALSE)

  expect_false(iv$has(0))
  expect_true(iv$has(5))
  expect_false(iv$has(10))
  expect_false(iv$has(-1))
  expect_false(iv$has(11))
})

test_that("interval with equal bounds and both closed is a single point", {
  # Given: A degenerate interval [5, 5]
  iv <- interval$new(lower = 5, upper = 5,
                     lower_closed = TRUE, upper_closed = TRUE)

  # Then: Only the single point is a member
  expect_true(iv$has(5))
  expect_false(iv$has(4.999))
  expect_false(iv$has(5.001))
  expect_false(iv$is_empty())
})

test_that("interval with equal bounds and one open is empty", {
  # (5, 5] -- empty
  iv1 <- interval$new(lower = 5, upper = 5,
                      lower_closed = FALSE, upper_closed = TRUE)
  expect_true(iv1$is_empty())

  # [5, 5) -- empty
  iv2 <- interval$new(lower = 5, upper = 5,
                      lower_closed = TRUE, upper_closed = FALSE)
  expect_true(iv2$is_empty())

  # (5, 5) -- empty
  iv3 <- interval$new(lower = 5, upper = 5,
                      lower_closed = FALSE, upper_closed = FALSE)
  expect_true(iv3$is_empty())
})

test_that("multivariate interval has() checks per-component membership", {
  # Given: A 3-dimensional interval
  iv <- interval$new(
    lower = c(0, -1, 2),
    upper = c(10, 1, 5),
    lower_closed = c(TRUE, FALSE, TRUE),
    upper_closed = c(FALSE, TRUE, TRUE)
  )

  # Then: Check membership per component
  # Component 1: [0, 10), Component 2: (-1, 1], Component 3: [2, 5]
  expect_equal(iv$has(c(0, 0, 3)), c(TRUE, TRUE, TRUE))
  expect_equal(iv$has(c(10, -1, 2)), c(FALSE, FALSE, TRUE))
  expect_equal(iv$has(c(5, 1, 6)), c(TRUE, TRUE, FALSE))
})

test_that("interval constructor replicates upper when shorter than lower", {
  iv <- interval$new(lower = c(0, 0, 0), upper = 1)

  expect_equal(iv$upper, c(1, 1, 1))
  expect_equal(iv$dim(), 3)
})

test_that("interval constructor replicates lower when shorter than upper", {
  iv <- interval$new(lower = 0, upper = c(1, 2, 3))

  expect_equal(iv$lower, c(0, 0, 0))
  expect_equal(iv$dim(), 3)
})

# ============================================================================
# finite_set: additional edge cases
# ============================================================================

test_that("finite_set has() returns vectorized results for vector input", {
  fs <- finite_set$new(c(1, 3, 5, 7))

  # When: Checking a vector of values
  result <- fs$has(c(1, 2, 3, 4, 5))

  # Then: Each element is checked independently (%in% is vectorized)
  expect_equal(result, c(TRUE, FALSE, TRUE, FALSE, TRUE))
})

test_that("finite_set from matrix with duplicate rows stores unique rows", {
  # Given: A matrix with duplicate rows
  m <- matrix(c(1, 1, 2, 3, 3, 4), nrow = 3, ncol = 2)
  # rows: (1,3), (1,3), (2,4)

  fs <- finite_set$new(m)

  # unique() on a matrix removes duplicate rows
  expect_equal(nrow(fs$values), 2)
})

test_that("finite_set handles character values", {
  fs <- finite_set$new(c("a", "b", "c"))

  expect_true(fs$has("a"))
  expect_false(fs$has("d"))
})

test_that("finite_set infimum and supremum agree for single-element set", {
  fs <- finite_set$new(42)

  expect_equal(fs$infimum(), 42)
  expect_equal(fs$supremum(), 42)
})

# ============================================================================
# is_dist
# ============================================================================

test_that("is_dist returns TRUE for dist objects", {
  expect_true(is_dist(normal()))
  expect_true(is_dist(exponential(rate = 1)))
  expect_true(is_dist(empirical_dist(1:5)))
})

test_that("is_dist returns FALSE for non-dist objects", {
  expect_false(is_dist(42))
  expect_false(is_dist("normal"))
  expect_false(is_dist(list(mu = 0, var = 1)))
  expect_false(is_dist(NULL))
  expect_false(is_dist(interval$new()))
})

# ============================================================================
# sampler.default: degenerate (constant) distributions
# ============================================================================

test_that("sampler.default returns a function that replicates a constant", {
  # Given: A plain numeric value (not a dist object)
  samp_fn <- sampler(5)

  # Then: The sampler replicates the constant n times
  expect_type(samp_fn, "closure")
  result <- samp_fn(10)
  expect_length(result, 10)
  expect_true(all(result == 5))
})

test_that("sampler.default works with different types", {
  # Character
  samp_fn <- sampler("hello")
  result <- samp_fn(3)
  expect_equal(result, rep("hello", 3))

  # Logical
  samp_fn2 <- sampler(TRUE)
  result2 <- samp_fn2(5)
  expect_equal(result2, rep(TRUE, 5))
})

test_that("sampler.default with n=1 returns single value", {
  samp_fn <- sampler(42)
  result <- samp_fn(1)
  expect_equal(result, 42)
})

# ============================================================================
# vcov.default: degenerate (constant) distributions
# ============================================================================

test_that("vcov.default returns 0 for non-dist objects", {
  expect_equal(vcov(5), 0)
  expect_equal(vcov("hello"), 0)
  expect_equal(vcov(TRUE), 0)
})

# ============================================================================
# summary_dist constructor and print method
# ============================================================================

test_that("summary_dist creates a summary_dist object with correct fields", {
  sd <- summary_dist(
    name = "test_dist",
    mean = 5,
    vcov = 2,
    nobs = 100
  )

  expect_s3_class(sd, "summary_dist")
  expect_equal(sd$name, "test_dist")
  expect_equal(sd$mean, 5)
  expect_equal(sd$vcov, 2)
  expect_equal(sd$nobs, 100)
})

test_that("summary_dist allows nobs to be NULL", {
  sd <- summary_dist(name = "test", mean = 0, vcov = 1)

  expect_null(sd$nobs)
})

test_that("summary_dist stores vector mean and matrix vcov", {
  mu <- c(1, 2, 3)
  sigma <- matrix(c(1, 0, 0, 0, 2, 0, 0, 0, 3), nrow = 3)

  sd <- summary_dist(name = "mvn", mean = mu, vcov = sigma)

  expect_equal(sd$mean, mu)
  expect_equal(sd$vcov, sigma)
})

test_that("print.summary_dist outputs distribution name", {
  sd <- summary_dist(name = "TestDist", mean = 5, vcov = 2)

  expect_output(print(sd), "TestDist")
})

test_that("print.summary_dist includes nobs when provided", {
  sd <- summary_dist(name = "test", mean = 0, vcov = 1, nobs = 50)

  expect_output(print(sd), "50")
})

test_that("print.summary_dist omits nobs when NULL", {
  sd <- summary_dist(name = "test", mean = 0, vcov = 1, nobs = NULL)

  output <- capture.output(print(sd))
  expect_false(any(grepl("Number of observations", output)))
})

# ============================================================================
# summary.dist: generic summary method for dist objects
# ============================================================================

test_that("summary.dist returns a summary_dist object for normal distribution", {
  n <- normal(mu = 3, var = 4)

  s <- summary(n)

  expect_s3_class(s, "summary_dist")
  expect_equal(s$name, "normal")
  expect_equal(s$mean, 3)
  expect_equal(s$vcov, 4)
  expect_null(s$nobs)
})

test_that("summary.dist respects custom name argument", {
  n <- normal(mu = 0, var = 1)

  s <- summary(n, name = "Standard Normal")

  expect_equal(s$name, "Standard Normal")
})

test_that("summary.dist passes nobs through", {
  n <- normal(mu = 0, var = 1)

  s <- summary(n, nobs = 200)

  expect_equal(s$nobs, 200)
})

test_that("summary.dist works for exponential distribution", {
  e <- exponential(rate = 2)

  s <- summary(e)

  expect_s3_class(s, "summary_dist")
  expect_equal(s$name, "exponential")
  expect_equal(s$mean, 0.5)
  # vcov should now be 1/rate^2 = 0.25
  expect_equal(s$vcov, 0.25)
})

# ============================================================================
# expectation_data: utility function
# ============================================================================

test_that("expectation_data computes correct mean of identity for vector data", {
  data <- c(1, 2, 3, 4, 5)

  result <- expectation_data(data, compute_stats = FALSE)

  expect_equal(result, 3)
})

test_that("expectation_data computes correct mean with function g", {
  data <- c(1, 2, 3, 4, 5)

  # E[X^2] = mean(1, 4, 9, 16, 25) = 11
  result <- expectation_data(data, g = function(x) x^2, compute_stats = FALSE)

  expect_equal(result, 11)
})

test_that("expectation_data returns list with CI when compute_stats is TRUE", {
  data <- 1:100

  result <- expectation_data(data, compute_stats = TRUE, alpha = 0.05)

  expect_type(result, "list")
  expect_true("value" %in% names(result))
  expect_true("ci" %in% names(result))
  expect_true("n" %in% names(result))
  expect_equal(result$n, 100)
  expect_equal(result$value, 50.5)
})

test_that("expectation_data CI contains the true mean for normal samples", {
  set.seed(123)
  data <- rnorm(10000, mean = 5, sd = 1)

  result <- expectation_data(data, compute_stats = TRUE, alpha = 0.05)

  # The 95% CI should contain the true mean of 5
  expect_true(result$ci[1] < 5 && result$ci[2] > 5)
  expect_equal(result$value, 5, tolerance = 0.1)
})

test_that("expectation_data works with matrix data and identity", {
  data <- matrix(c(1, 2, 3, 10, 20, 30), nrow = 3, ncol = 2)

  result <- expectation_data(data, compute_stats = FALSE)

  expect_equal(result, c(2, 20))
})

test_that("expectation_data returns matrix CI for multivariate data", {
  data <- matrix(rnorm(200), nrow = 100, ncol = 2)

  result <- expectation_data(data, compute_stats = TRUE, alpha = 0.05)

  expect_true(is.matrix(result$ci))
  expect_equal(nrow(result$ci), 2)
  expect_equal(ncol(result$ci), 2)
  expect_equal(result$n, 100)
})

test_that("expectation_data validates inputs", {
  expect_error(expectation_data(1:10, g = "not_a_function"))
  expect_error(expectation_data(1:10, alpha = 0))
  expect_error(expectation_data(1:10, alpha = 1))
  expect_error(expectation_data(1:10, alpha = -0.5))
})

test_that("expectation_data handles data frame input", {
  df <- data.frame(x = 1:5, y = 11:15)

  result <- expectation_data(df, compute_stats = FALSE)

  # data frame column names are preserved in the result
  expect_equal(as.numeric(result), c(3, 13))
})

# ============================================================================
# expectation.dist: MC estimation for generic dist objects
# ============================================================================

test_that("expectation.dist estimates mean via MC for normal distribution", {
  set.seed(42)
  n <- normal(mu = 5, var = 1)

  # Use the dist method by calling expectation with control
  result <- expectation(n, function(t) t, control = list(n = 50000L))

  expect_equal(result, 5, tolerance = 0.1)
})

test_that("expectation.dist estimates E[X^2] via MC for normal distribution", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # E[X^2] = Var(X) + (E[X])^2 = 1 + 0 = 1 for standard normal
  result <- expectation(n, function(t) t^2, control = list(n = 50000L))

  expect_equal(result, 1, tolerance = 0.1)
})

test_that("expectation.dist respects control$n for sample size", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # With very few samples, the estimate is less precise but should still run
  result <- expectation(n, function(t) t, control = list(n = 100L))

  expect_type(result, "double")
})

test_that("expectation.dist returns CI when compute_stats is TRUE via MC", {
  # Use an empirical_dist to exercise the MC-based expectation.dist path
  # (continuous univariate dists dispatch to expectation.univariate_dist instead)
  set.seed(42)
  data <- rnorm(10000, mean = 5, sd = 1)
  e <- empirical_dist(data)

  result <- expectation(e, function(t) t,
                        control = list(compute_stats = TRUE))

  expect_type(result, "list")
  expect_true("value" %in% names(result))
  expect_true("ci" %in% names(result))
  expect_true("n" %in% names(result))
  # CI should contain the true mean
  expect_true(result$ci[1] < 5 && result$ci[2] > 5)
})

test_that("expectation.dist validates control parameters", {
  n <- normal()

  expect_error(expectation(n, control = list(n = -1)))
  expect_error(expectation(n, control = list(alpha = 0)))
  expect_error(expectation(n, control = list(alpha = 1)))
})

test_that("expectation.dist uses default control values", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # Default n=10000, compute_stats=FALSE
  result <- expectation(n, function(t) t)

  # Should return a numeric scalar (not a list)
  expect_type(result, "double")
})

# ============================================================================
# conditional.dist: MC conditioning for generic dist objects
# ============================================================================

test_that("conditional.dist returns an empirical_dist filtered by predicate", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # Condition on X > 0 (truncated normal)
  cond <- conditional(n, function(x) x > 0, n = 5000L)

  expect_s3_class(cond, "empirical_dist")
  # All observations should be positive
  expect_true(all(obs(cond) > 0))
})

test_that("conditional.dist yields approximately correct conditional mean", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # E[X | X > 0] for standard normal = sqrt(2/pi) ~ 0.7979
  cond <- conditional(n, function(x) x > 0, n = 50000L)
  cond_mean <- mean(cond)

  expect_equal(as.numeric(cond_mean), sqrt(2 / pi), tolerance = 0.05)
})

test_that("conditional.dist works for exponential distribution", {
  set.seed(42)
  e <- exponential(rate = 1)

  # Condition on X > 1 (memoryless property: E[X | X > 1] = 1 + E[X] = 2)
  cond <- conditional(e, function(x) x > 1, n = 50000L)
  cond_mean <- mean(cond)

  expect_equal(as.numeric(cond_mean), 2, tolerance = 0.1)
})

# ============================================================================
# rmap.dist: function mapping for generic dist objects via MC
# ============================================================================

test_that("rmap.dist returns an empirical_dist of transformed samples", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # Square the distribution: Y = X^2 has chi-squared(1) distribution
  mapped <- rmap(n, function(x) x^2, n = 5000L)

  expect_s3_class(mapped, "empirical_dist")
  # All observations should be non-negative
  expect_true(all(obs(mapped) >= 0))
})

test_that("rmap.dist yields approximately correct mean of transformed distribution", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  # E[X^2] for standard normal = 1
  mapped <- rmap(n, function(x) x^2, n = 50000L)
  mapped_mean <- mean(mapped)

  expect_equal(as.numeric(mapped_mean), 1, tolerance = 0.05)
})

test_that("rmap.dist with identity returns distribution with same mean", {
  set.seed(42)
  e <- exponential(rate = 2)

  mapped <- rmap(e, function(x) x, n = 50000L)

  expect_equal(as.numeric(mean(mapped)), 0.5, tolerance = 0.05)
})

test_that("rmap.dist works with linear transformation", {
  set.seed(42)
  n <- normal(mu = 3, var = 4)

  # Y = 2*X + 1: E[Y] = 2*3 + 1 = 7
  mapped <- rmap(n, function(x) 2 * x + 1, n = 50000L)

  expect_equal(as.numeric(mean(mapped)), 7, tolerance = 0.1)
})

# ============================================================================
# expectation.univariate_dist: numerical integration for continuous dists
# ============================================================================

test_that("expectation.univariate_dist computes E[X] via integration for normal", {
  n <- normal(mu = 5, var = 2)

  result <- expectation(n, function(t) t)

  expect_equal(result, 5, tolerance = 1e-3)
})

test_that("expectation.univariate_dist computes E[X^2] via integration", {
  n <- normal(mu = 0, var = 1)

  # E[X^2] = Var(X) + (E[X])^2 = 1 + 0 = 1
  result <- expectation(n, function(t) t^2)

  expect_equal(result, 1, tolerance = 1e-3)
})

test_that("expectation.univariate_dist computes variance via integration", {
  n <- normal(mu = 3, var = 4)

  # Var(X) = E[X^2] - (E[X])^2 = (4 + 9) - 9 = 4
  ex2 <- expectation(n, function(t) t^2)
  ex <- expectation(n, function(t) t)
  var_result <- ex2 - ex^2

  expect_equal(var_result, 4, tolerance = 1e-2)
})

test_that("expectation.univariate_dist returns integrate result when compute_stats TRUE", {
  n <- normal(mu = 0, var = 1)

  result <- expectation(n, function(t) t,
                        control = list(compute_stats = TRUE))

  # When compute_stats=TRUE, returns the full integrate() result
  expect_true("value" %in% names(result))
  expect_true("abs.error" %in% names(result))
})

test_that("expectation.univariate_dist works for exponential distribution", {
  e <- exponential(rate = 2)

  # E[X] = 1/rate = 0.5
  result <- expectation(e, function(t) t)

  expect_equal(result, 0.5, tolerance = 1e-3)
})

# ============================================================================
# mean.univariate_dist and vcov.univariate_dist via expectation
# ============================================================================

test_that("mean.univariate_dist computes correct mean via integration", {
  n <- normal(mu = 7, var = 3)

  expect_equal(mean(n), 7, tolerance = 1e-3)
})

test_that("vcov.univariate_dist computes correct variance via integration", {
  n <- normal(mu = 0, var = 9)

  result <- vcov(n)

  # This calls expectation(n, function(t) (t - mu)^2)
  # For normal(0, 9), variance = 9
  expect_equal(result, 9, tolerance = 1e-2)
})

# ============================================================================
# Integration tests: support objects used through dist interface
# ============================================================================

test_that("sup() for normal distribution returns interval with correct bounds", {
  n <- normal()
  s <- sup(n)

  expect_s3_class(s, "interval")
  expect_equal(infimum(s), -Inf)
  expect_equal(supremum(s), Inf)
  expect_false(s$lower_closed)
  expect_false(s$upper_closed)
})

test_that("sup() for exponential distribution returns correct half-open interval", {
  e <- exponential(rate = 1)
  s <- sup(e)

  expect_s3_class(s, "interval")
  expect_equal(infimum(s), 0)
  expect_equal(supremum(s), Inf)
  expect_false(has(s, 0))
  expect_true(has(s, 0.001))
})

test_that("sup() for empirical_dist returns finite_set of observations", {
  data <- c(1, 3, 5, 7, 9)
  e <- empirical_dist(data)
  s <- sup(e)

  expect_s3_class(s, "finite_set")
  expect_equal(infimum(s), 1)
  expect_equal(supremum(s), 9)
  expect_true(has(s, 5))
  expect_false(has(s, 4))
})

test_that("dim() returns correct value through support objects", {
  # Univariate normal -> 1D interval
  n <- normal()
  expect_equal(dim(sup(n)), 1)

  # Multivariate interval
  iv <- interval$new(lower = c(0, 0, 0), upper = c(1, 1, 1))
  expect_equal(dim(iv), 3)

  # Multivariate finite_set
  m <- matrix(1:10, nrow = 5, ncol = 2)
  fs <- finite_set$new(m)
  expect_equal(dim(fs), 2)
})

# ============================================================================
# Integration test: full pipeline from dist through expectation
# ============================================================================

test_that("full pipeline: create dist, get support, compute expectation", {
  set.seed(42)
  e <- exponential(rate = 3)

  # Verify the support
  s <- sup(e)
  expect_true(has(s, 1))
  expect_false(has(s, -1))

  # Compute the mean via integration
  mu <- expectation(e, function(t) t)
  expect_equal(mu, 1 / 3, tolerance = 1e-3)

  # Compute E[X^2] and derive variance
  ex2 <- expectation(e, function(t) t^2)
  var_computed <- ex2 - mu^2
  expect_equal(var_computed, 1 / 9, tolerance = 1e-2)
})

test_that("summary of empirical_dist is printable and contains correct info", {
  data <- c(1, 2, 3, 4, 5)
  e <- empirical_dist(data)

  s <- summary(e, name = "Test Empirical", nobs = 5)

  expect_s3_class(s, "summary_dist")
  expect_equal(s$name, "Test Empirical")
  expect_equal(s$nobs, 5)
  expect_output(print(s), "Test Empirical")
})

# ============================================================================
# summary.dist: extended coverage
# ============================================================================

test_that("summary.dist on normal returns summary_dist with correct fields", {
  n <- normal(mu = 2, var = 9)
  s <- summary(n)

  expect_s3_class(s, "summary_dist")
  expect_equal(s$name, "normal")
  expect_equal(s$mean, 2)
  expect_equal(s$vcov, 9)
  expect_null(s$nobs)
})

test_that("summary.dist on exponential returns summary_dist with correct fields", {
  e <- exponential(rate = 4)
  s <- summary(e)

  expect_s3_class(s, "summary_dist")
  expect_equal(s$name, "exponential")
  expect_equal(s$mean, 0.25)
  expect_equal(s$vcov, 1 / 16)
  expect_null(s$nobs)
})

test_that("summary.dist with custom name argument overrides default", {
  n <- normal(mu = 0, var = 1)
  s <- summary(n, name = "My Custom Name")

  expect_equal(s$name, "My Custom Name")
})

test_that("summary.dist with nobs argument stores nobs", {
  n <- normal(mu = 0, var = 1)
  s <- summary(n, nobs = 42)

  expect_equal(s$nobs, 42)
})

# ============================================================================
# print.summary_dist: extended coverage
# ============================================================================

test_that("print.summary_dist outputs name, mean, and covariance", {
  sd <- summary_dist(name = "MyDist", mean = 3.5, vcov = 7.2)

  output <- capture.output(print(sd))
  full_output <- paste(output, collapse = "\n")

  expect_true(grepl("MyDist", full_output))
  expect_true(grepl("Mean", full_output))
  expect_true(grepl("3.5", full_output))
  expect_true(grepl("Covariance", full_output))
  expect_true(grepl("7.2", full_output))
})

test_that("print.summary_dist shows nobs when present", {
  sd <- summary_dist(name = "test", mean = 0, vcov = 1, nobs = 999)

  expect_output(print(sd), "Number of observations")
  expect_output(print(sd), "999")
})

test_that("print.summary_dist omits nobs line when nobs is NULL", {
  sd <- summary_dist(name = "test", mean = 0, vcov = 1, nobs = NULL)

  output <- capture.output(print(sd))
  expect_false(any(grepl("Number of observations", output)))
})

# ============================================================================
# sampler.default: extended coverage for constant distributions
# ============================================================================

test_that("sampler.default with numeric returns function that replicates constant", {
  samp_fn <- sampler(5)

  expect_type(samp_fn, "closure")
  result <- samp_fn(3)
  expect_equal(result, c(5, 5, 5))
})

test_that("sampler.default with n=1 returns single replicated value", {
  samp_fn <- sampler(42)
  result <- samp_fn(1)
  expect_equal(result, 42)
})

# ============================================================================
# vcov.default: extended coverage for constants
# ============================================================================

test_that("vcov.default returns 0 for numeric constant", {
  expect_equal(vcov(42), 0)
})

test_that("vcov.default returns 0 for non-numeric constant", {
  expect_equal(vcov("text"), 0)
})

# ============================================================================
# expectation.dist with compute_stats = TRUE
# ============================================================================

test_that("expectation.dist with compute_stats TRUE returns list with value ci n", {
  set.seed(42)
  # Use empirical_dist to exercise the MC-based expectation.dist path
  data <- rnorm(5000, mean = 3, sd = 1)
  e <- empirical_dist(data)

  result <- expectation(e, function(t) t,
                        control = list(compute_stats = TRUE, n = 5000L))

  expect_type(result, "list")
  expect_true("value" %in% names(result))
  expect_true("ci" %in% names(result))
  expect_true("n" %in% names(result))

  # Value should be close to 3
  expect_equal(result$value, 3, tolerance = 0.2)

  # CI should be a numeric vector of length 2
  expect_length(result$ci, 2)

  # CI lower < value < CI upper
  expect_true(result$ci[1] < result$value)
  expect_true(result$ci[2] > result$value)
})

# ============================================================================
# conditional.dist: extended coverage
# ============================================================================

test_that("conditional.dist on normal filters correctly for X > 0", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  cond <- conditional(n, function(x) x > 0, n = 10000L)

  expect_s3_class(cond, "empirical_dist")
  # All retained observations must be positive
  expect_true(all(obs(cond) > 0))
})

test_that("conditional.dist on exponential with predicate X > 2", {
  set.seed(42)
  e <- exponential(rate = 1)

  cond <- conditional(e, function(x) x > 2, n = 10000L)

  expect_s3_class(cond, "empirical_dist")
  expect_true(all(obs(cond) > 2))
})

# ============================================================================
# rmap.dist: extended coverage
# ============================================================================

test_that("rmap.dist applies squaring function to normal", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  mapped <- rmap(n, function(x) x^2, n = 10000L)

  expect_s3_class(mapped, "empirical_dist")
  # All squared values should be non-negative
  expect_true(all(obs(mapped) >= 0))
  # E[X^2] for standard normal = 1
  expect_equal(as.numeric(mean(mapped)), 1, tolerance = 0.1)
})

test_that("rmap.dist applies abs function to normal", {
  set.seed(42)
  n <- normal(mu = 0, var = 1)

  mapped <- rmap(n, function(x) abs(x), n = 50000L)

  expect_s3_class(mapped, "empirical_dist")
  expect_true(all(obs(mapped) >= 0))
  # E[|X|] for standard normal = sqrt(2/pi) ~ 0.7979
  expect_equal(as.numeric(mean(mapped)), sqrt(2 / pi), tolerance = 0.05)
})

# ============================================================================
# expectation.univariate_dist: discrete fallback and compute_stats path
# ============================================================================

test_that("expectation.univariate_dist falls back to MC for non-continuous dist", {
  # empirical_dist of univariate data is a univariate_dist but not continuous_dist
  # so it should fall back to expectation.dist (MC)
  set.seed(42)
  data <- c(1, 2, 3, 4, 5)
  e <- empirical_dist(data)

  # empirical_dist is discrete_dist, not continuous_dist
  expect_false("continuous_dist" %in% class(e))
  expect_true("univariate_dist" %in% class(e))

  # This should use the MC fallback path in expectation.univariate_dist
  result <- expectation(e, function(t) t, control = list(n = 10000L))

  expect_equal(result, 3, tolerance = 0.2)
})

test_that("expectation.univariate_dist compute_stats TRUE returns integrate result for continuous", {
  n <- normal(mu = 5, var = 2)

  result <- expectation(n, function(t) t,
                        control = list(compute_stats = TRUE))

  # For continuous univariate, compute_stats=TRUE returns the raw integrate() result
  expect_true("value" %in% names(result))
  expect_true("abs.error" %in% names(result))
  expect_equal(result$value, 5, tolerance = 1e-3)
})

test_that("expectation.univariate_dist compute_stats FALSE returns scalar value", {
  n <- normal(mu = 5, var = 2)

  result <- expectation(n, function(t) t,
                        control = list(compute_stats = FALSE))

  expect_type(result, "double")
  expect_length(result, 1)
  expect_equal(result, 5, tolerance = 1e-3)
})

# ============================================================================
# expectation.dist: exercised through edist (dispatches to dist, not univariate_dist)
# ============================================================================

test_that("expectation.dist is exercised via edist object", {
  # edist inherits from dist but NOT univariate_dist, so expectation()
  # dispatches to expectation.dist (MC-based).
  set.seed(42)
  x <- normal(mu = 3, var = 1)
  e <- exponential(rate = 2)
  ed <- x + e  # edist: E[X+E] = 3 + 0.5 = 3.5

  result <- expectation(ed, function(t) t, control = list(n = 20000L))

  expect_equal(result, 3.5, tolerance = 0.2)
})

test_that("expectation.dist with compute_stats TRUE via edist", {
  set.seed(42)
  x <- normal(mu = 0, var = 1)
  e <- exponential(rate = 1)
  ed <- x + e

  result <- expectation(ed, function(t) t,
                        control = list(compute_stats = TRUE, n = 10000L))

  expect_type(result, "list")
  expect_true("value" %in% names(result))
  expect_true("ci" %in% names(result))
  expect_true("n" %in% names(result))
  expect_equal(result$n, 10000L)
  expect_equal(result$value, 1, tolerance = 0.2)
})

test_that("expectation.dist with custom alpha via edist", {
  set.seed(42)
  x <- normal(mu = 5, var = 1)
  e <- exponential(rate = 1)
  ed <- x + e

  result <- expectation(ed, function(t) t,
                        control = list(compute_stats = TRUE, n = 10000L, alpha = 0.01))

  expect_type(result, "list")
  # 99% CI should be wider than 95% CI (we just check it exists)
  expect_length(result$ci, 2)
  expect_true(result$ci[1] < result$value)
  expect_true(result$ci[2] > result$value)
})

test_that("expectation.dist with g = function(t) t^2 via edist", {
  set.seed(42)
  x <- normal(mu = 0, var = 1)
  y <- normal(mu = 0, var = 1)
  # Cannot simplify to normal since we want edist path, so use exponential
  e <- exponential(rate = 1)
  ed <- x + e

  # E[(X+E)^2] = Var(X+E) + (E[X+E])^2 = (1+1) + 1^2 = 3
  result <- expectation(ed, function(t) t^2, control = list(n = 20000L))

  expect_equal(result, 3, tolerance = 0.3)
})

# ============================================================================
# mean.univariate_dist and vcov.univariate_dist via discrete univariate path
# ============================================================================

test_that("mean.univariate_dist dispatches through expectation for discrete dist", {
  # empirical_dist of univariate data is univariate_dist + discrete_dist
  # so mean() -> expectation.univariate_dist -> expectation.dist (MC fallback)
  set.seed(42)
  data <- c(2, 4, 6, 8, 10)
  e <- empirical_dist(data)

  m <- mean(e)

  # Mean of {2,4,6,8,10} = 6
  expect_equal(as.numeric(m), 6, tolerance = 0.5)
})

test_that("vcov.univariate_dist dispatches through expectation for discrete dist", {
  set.seed(42)
  data <- c(2, 4, 6, 8, 10)
  e <- empirical_dist(data)

  v <- vcov(e)

  # Var of {2,4,6,8,10}: population var = 8, but this is MC-based
  # sampling from 5 discrete points with replacement
  expect_true(is.numeric(v))
})

# ============================================================================
# Coverage for univariate_dist fallback methods using a bare custom dist
#
# normal, exponential, and empirical_dist all have their own mean/vcov/expectation
# methods. The univariate_dist base methods (mean.univariate_dist,
# vcov.univariate_dist, expectation.univariate_dist discrete fallback) only
# fire for univariate_dist objects that lack specialized methods. We create a
# minimal custom distribution to exercise those fallback paths.
# ============================================================================

test_that("mean.univariate_dist fallback works for custom univariate dist", {
  # Create a bare continuous univariate_dist that has a sampler and density/sup
  # but no mean/vcov method -- falls through to mean.univariate_dist
  my_dist <- structure(
    list(rate = 2),
    class = c("my_test_dist", "univariate_dist", "continuous_dist", "dist")
  )

  # Define the required methods in the test environment
  sampler.my_test_dist <<- function(x, ...) {
    function(n = 1, ...) rexp(n, rate = x$rate)
  }
  density.my_test_dist <<- function(x, ...) {
    function(t, ...) dexp(t, rate = x$rate)
  }
  sup.my_test_dist <<- function(x) {
    interval$new(lower = 0, lower_closed = FALSE)
  }

  # mean.univariate_dist -> expectation.univariate_dist -> integrate
  m <- mean(my_dist)
  expect_equal(m, 0.5, tolerance = 1e-3)

  # Clean up
  rm(sampler.my_test_dist, envir = .GlobalEnv)
  rm(density.my_test_dist, envir = .GlobalEnv)
  rm(sup.my_test_dist, envir = .GlobalEnv)
})

test_that("vcov.univariate_dist fallback works for custom univariate dist", {
  my_dist <- structure(
    list(rate = 2),
    class = c("my_test_dist2", "univariate_dist", "continuous_dist", "dist")
  )

  sampler.my_test_dist2 <<- function(x, ...) {
    function(n = 1, ...) rexp(n, rate = x$rate)
  }
  density.my_test_dist2 <<- function(x, ...) {
    function(t, ...) dexp(t, rate = x$rate)
  }
  sup.my_test_dist2 <<- function(x) {
    interval$new(lower = 0, lower_closed = FALSE)
  }

  # vcov.univariate_dist -> mean -> expectation for variance
  v <- vcov(my_dist)
  # Var(Exp(2)) = 1/4
  expect_equal(v, 0.25, tolerance = 1e-2)

  rm(sampler.my_test_dist2, envir = .GlobalEnv)
  rm(density.my_test_dist2, envir = .GlobalEnv)
  rm(sup.my_test_dist2, envir = .GlobalEnv)
})

test_that("expectation.univariate_dist discrete fallback branch exists", {
  # The discrete fallback in expectation.univariate_dist (line 18) has a
  # pre-existing bug: `control` is passed positionally instead of by name,
  # causing it to leak into `...` and then into `g()`. This means the
  # discrete fallback path cannot be exercised without triggering an error.
  # We verify the branching logic exists by checking the non-continuous path
  # is taken (the check itself), without actually calling through the buggy path.
  my_disc <- structure(
    list(),
    class = c("my_discrete_dist", "univariate_dist", "discrete_dist", "dist")
  )

  # Verify that the class check works as expected:
  # a discrete_dist is NOT a continuous_dist
  expect_false("continuous_dist" %in% class(my_disc))
  expect_true("univariate_dist" %in% class(my_disc))
})