test_that("aedl validates basic inputs", {
  expect_error(aedl(c(1, 2), c(1), nsim = 1, nsim_chi = 1, n_tau = 1),
               "same length")
  expect_error(aedl(c(1, 2), c(1, 0), nsim = 1, nsim_chi = 1, n_tau = 1),
               "positive")
  expect_error(aedl(c(1, NA), c(1, 1), nsim = 1, nsim_chi = 1, n_tau = 1),
               "missing")
})

test_that("plugin method returns a valid aedl object", {
  set.seed(1)

  res <- aedl(
    yi = c(0.10, 0.25, 0.05, 0.30),
    vi = c(0.04, 0.05, 0.03, 0.06),
    method = "plugin",
    t_grid = seq(-6, 6, by = 0.2),
    nsim = 100,
    nsim_chi = 500,
    n_tau = 50
  )

  expect_s3_class(res, "aedl")
  expect_equal(res$method, "plugin")
  expect_true(is.finite(res$estimate))
  expect_true(is.finite(res$se))
  expect_true(res$ci.lb < res$ci.ub)
  expect_true(res$pval >= 0 && res$pval <= 1)
  expect_true(res$tau2 >= 0)
  expect_true(is.na(res$I2c))
})

test_that("i2c method returns a valid aedl object", {
  set.seed(2)

  res <- aedl(
    yi = c(0.10, 0.25, 0.05, 0.30),
    vi = c(0.04, 0.05, 0.03, 0.06),
    method = "i2c",
    a = 0,
    b = 1,
    t_grid = seq(-6, 6, by = 0.2),
    I2_grid = seq(0, 0.95, by = 0.05),
    nsim = 100,
    nsim_chi = 500,
    n_tau = 50
  )

  expect_s3_class(res, "aedl")
  expect_equal(res$method, "i2c")
  expect_true(is.finite(res$estimate))
  expect_true(is.finite(res$se))
  expect_true(is.finite(res$ci.lb))
  expect_true(is.finite(res$ci.ub))
  expect_true(is.finite(res$pval))
  expect_true(res$ci.lb < res$ci.ub)
  expect_true(res$pval >= 0 && res$pval <= 1)
  expect_true(res$tau2 >= 0)
  expect_true(res$I2hat >= 0 && res$I2hat < 1)
  expect_true(is.finite(res$I2c))
  expect_true(res$I2c >= 0 && res$I2c < 1)
  expect_equal(res$tau2_for_dist, res$tau2c)
})

test_that("print and confint methods work", {
  res <- aedl(
    yi = c(0.10, 0.25, 0.05, 0.30),
    vi = c(0.04, 0.05, 0.03, 0.06),
    method = "plugin",
    t_grid = seq(-6, 6, by = 0.2),
    nsim = 100,
    nsim_chi = 500,
    n_tau = 50,
    seed = 3
  )

  expect_output(print(res), "Almost-exact")
  ci <- confint(res)
  expect_true(is.matrix(ci))
  expect_equal(nrow(ci), 1)
  expect_equal(ncol(ci), 2)
  expect_true(ci[1, 1] < ci[1, 2])
})

test_that("README-style lightweight example works", {
  yi <- c(0.10, 0.25, 0.05, 0.30)
  vi <- c(0.04, 0.05, 0.03, 0.06)

  res <- aedl(
    yi,
    vi,
    t_grid = seq(-6, 6, by = 0.2),
    nsim = 100,
    nsim_chi = 500,
    n_tau = 50,
    seed = 1
  )

  expect_s3_class(res, "aedl")
  expect_true(is.finite(res$pval))
  expect_true(res$ci.lb < res$ci.ub)
})

test_that("density helper functions return finite density values", {
  vi <- c(0.04, 0.05, 0.03, 0.06)

  set.seed(4)
  tau_density <- aedl_tau2_density(
    vi = vi,
    tau2 = 0.01,
    nsim_chi = 500,
    n_tau = 50
  )

  expect_type(tau_density, "list")
  expect_equal(length(tau_density$x), length(tau_density$y))
  expect_true(all(is.finite(tau_density$x)))
  expect_true(all(is.finite(tau_density$y)))

  set.seed(5)
  stat_density <- aedl_stat_density(
    t_grid = seq(-4, 4, by = 0.5),
    vi = vi,
    tau2 = 0.01,
    nsim = 100,
    nsim_chi = 500,
    n_tau = 50
  )

  expect_type(stat_density, "list")
  expect_equal(length(stat_density$x), length(stat_density$y))
  expect_true(all(is.finite(stat_density$x)))
  expect_true(all(is.finite(stat_density$y)))
})