# =======================================================================
# aiDIF tests
#
# 0.1.0 tested structure only. These add numerical checks, regression
# tests for the bugs fixed in 0.2.0, and coverage of the new arguments.
# =======================================================================

make_mle <- function(d_shift_g2 = 0, n_items = 4L, se_a = 0.06, se_d = 0.07,
                     seed = 1) {
  set.seed(seed)
  a <- runif(n_items, 0.8, 1.5)
  d <- rnorm(n_items, 0, 0.5)
  make_df <- function(av, dv) {
    df <- data.frame(a1 = av, d1 = dv)
    rownames(df) <- paste0("item", seq_len(n_items))
    df
  }
  make_vc <- function() {
    v <- diag(c(rbind(se_a^2, se_d^2)), n_items * 2)
    for (i in seq_len(n_items)) {
      ia <- (i - 1) * 2 + 1; id <- ia + 1
      v[ia, id] <- v[id, ia] <- -0.3 * se_a * se_d
    }
    v
  }
  pn <- c(rbind(paste0("item", seq_len(n_items), ".a1"),
                paste0("item", seq_len(n_items), ".d1")))
  d2 <- d - 0.5 + d_shift_g2 * c(1, rep(0, n_items - 1))
  list(par.names = list(internal = pn, original = pn),
       est     = list(group.1 = make_df(a, d), group.2 = make_df(a, d2)),
       var.cov = list(group.1 = make_vc(), group.2 = make_vc()))
}

# ---------------------------------------------------------------- validation

test_that("check_aidif_mle rejects non-lists and missing components", {
  expect_error(check_aidif_mle("nope"), "must be a list")
  expect_error(check_aidif_mle(list(est = list(g1 = data.frame(a1 = 1)))),
               "var.cov")
})

test_that("validation catches missing a1/d1 columns (new in 0.2.0)", {
  m <- make_mle()
  names(m$est$group.1) <- c("slope", "intercept")
  expect_error(check_aidif_mle(m), "no slope column")
})

test_that("validation catches non-consecutive intercept columns", {
  m <- make_mle()
  m$est$group.1 <- cbind(m$est$group.1, d3 = 0)
  m$est$group.2 <- cbind(m$est$group.2, d3 = 0)
  expect_error(check_aidif_mle(m), "consecutively|d1\\.\\.d2")
})

test_that("validation catches a covariance matrix of the wrong size", {
  m <- make_mle()
  m$var.cov$group.1 <- diag(2)
  expect_error(check_aidif_mle(m), "parameters")
})

test_that("compatible_mles compares item names, not just counts", {
  h <- make_mle(); a <- make_mle()
  rownames(a$est$group.1) <- rev(rownames(a$est$group.1))
  rownames(a$est$group.2) <- rownames(a$est$group.1)
  expect_error(check_compatible_mles(h, a), "Item names differ")
})

test_that("read_ai_scored returns aidif_data and fit_aidif accepts it", {
  m <- make_mle()
  dd <- read_ai_scored(m, m)
  expect_s3_class(dd, "aidif_data")
  expect_s3_class(fit_aidif(dd), "aidif")
  expect_error(fit_aidif(dd, m), "leave `ai_mle` as NULL")
})

# ------------------------------------------------------------ scoring bias

test_that("scoring_bias_test returns the documented structure", {
  out <- scoring_bias_test(make_mle(), make_mle(d_shift_g2 = 0.4))
  expect_s3_class(out, "dasb")
  expect_true(all(c("shift_g1", "shift_g2", "DASB", "se", "lower", "upper",
                    "z", "p_val") %in% names(out)))
  expect_equal(nrow(out), 4L)
})

test_that("results are NOT rounded at storage time", {
  out <- scoring_bias_test(make_mle(), make_mle(d_shift_g2 = 0.4))
  # Note: `se` is exactly 2 * se_d here (four equal variances), so it is
  # representable at 4dp and cannot discriminate. `z` and `p_val` can.
  expect_false(isTRUE(all.equal(out$z,     round(out$z, 4))))
  expect_false(isTRUE(all.equal(out$p_val, round(out$p_val, 4))))
})

test_that("uniform AI drift produces zero DASB to machine precision", {
  h <- make_mle()
  a <- h
  a$est$group.1$d1 <- h$est$group.1$d1 + 0.1
  a$est$group.2$d1 <- h$est$group.2$d1 + 0.1
  out <- scoring_bias_test(h, a)
  expect_true(all(abs(out$DASB) < 1e-12))
})

test_that("planted DASB is recovered in the point estimate", {
  out <- scoring_bias_test(make_mle(), make_mle(d_shift_g2 = 0.5))
  expect_equal(unname(out["item1", "DASB"]), 0.5, tolerance = 1e-8)
  expect_true(all(abs(out[-1, "DASB"]) < 1e-12))
  expect_lt(out["item1", "p_val"], 0.05)
})

test_that("p-values do not underflow to exactly zero (0.2.0 fix)", {
  # |z| ~ 10.7: 1 - pnorm(|z|) rounds to exactly 1 in double precision, so the
  # old 2 * (1 - pnorm(|z|)) returned 0, while 2 * pnorm(-|z|) is ~1e-26.
  # Beyond |z| ~ 37 even pnorm(-|z|) underflows, so a huge shift tests nothing.
  out <- scoring_bias_test(make_mle(), make_mle(d_shift_g2 = 1.5))
  expect_gt(abs(out["item1", "z"]), 8.5)
  expect_gt(out["item1", "p_val"], 0)
  expect_lt(out["item1", "p_val"], 1e-20)
  # Confirm the old formula really would have failed here.
  expect_identical(2 * (1 - pnorm(abs(out["item1", "z"]))), 0)
})

test_that("independent variance equals the sum of four marginals", {
  h <- make_mle(); a <- make_mle(d_shift_g2 = 0.4)
  out <- scoring_bias_test(h, a)
  v <- diag(h$var.cov$group.1)[2] + diag(h$var.cov$group.2)[2] +
       diag(a$var.cov$group.1)[2] + diag(a$var.cov$group.2)[2]
  expect_equal(unname(out["item1", "se"]), sqrt(v), tolerance = 1e-10)
})

test_that("paired design shrinks the SE relative to independent", {
  h <- make_mle(); a <- make_mle(d_shift_g2 = 0.4)
  cc <- list(group.1 = 0.5 * h$var.cov$group.1,
             group.2 = 0.5 * h$var.cov$group.2)
  ind <- scoring_bias_test(h, a)
  par <- scoring_bias_test(h, a, design = "paired", cross_cov = cc)
  expect_true(all(par$se < ind$se))
})

test_that("paired design requires and validates cross_cov", {
  h <- make_mle(); a <- make_mle(d_shift_g2 = 0.4)
  expect_error(scoring_bias_test(h, a, design = "paired"), "requires `cross_cov`")
  expect_error(
    scoring_bias_test(h, a, design = "paired",
                      cross_cov = list(group.1 = diag(2), group.2 = diag(2))),
    "must match")
})

test_that("linked metric removes an item-common offset", {
  h <- make_mle(n_items = 8)
  a <- make_mle(n_items = 8)
  # Uniform group-by-condition offset of 0.3 on every item, plus a real
  # effect of 0.6 on item 1 only.
  a$est$group.2$d1 <- a$est$group.2$d1 + 0.3
  a$est$group.2$d1[1] <- a$est$group.2$d1[1] + 0.6
  raw    <- scoring_bias_test(h, a, metric = "common")
  linked <- scoring_bias_test(h, a, metric = "linked")
  expect_equal(attr(linked, "offset"), 0.3, tolerance = 1e-6)
  expect_equal(unname(linked[1, "DASB"]), 0.6, tolerance = 1e-6)
  expect_true(all(abs(linked[-1, "DASB"]) < 1e-6))
  # The raw version wrongly attributes the offset to every item.
  expect_true(all(abs(raw[-1, "DASB"] - 0.3) < 1e-6))
})

test_that("multiplicity adjustment is applied and is never anti-conservative", {
  out <- scoring_bias_test(make_mle(), make_mle(d_shift_g2 = 0.4),
                           adjust = "BH")
  expect_true("p_adj" %in% names(out))
  expect_true(all(out$p_adj >= out$p_val - 1e-12))
  expect_error(scoring_bias_test(make_mle(), make_mle(), adjust = "nonsense"),
               "must be one of")
})

test_that("the removed `fun` argument now errors rather than being ignored", {
  expect_error(scoring_bias_test(make_mle(), make_mle(), fun = "d_fun3"))
})

# ------------------------------------------------------------------- engine

test_that("estimate_robust_scale returns class rdif with a rho curve", {
  fit <- estimate_robust_scale(make_mle(d_shift_g2 = 0.5))
  expect_s3_class(fit, "rdif")
  expect_true(is.data.frame(fit$rho.plot))
  expect_true(all(c("theta", "rho") %in% names(fit$rho.plot)))
  expect_true(any(is.finite(fit$rho.plot$rho)))
  expect_type(fit$multiple.solutions, "logical")
})

test_that("the rho grid brackets the estimate", {
  fit <- estimate_robust_scale(make_mle(d_shift_g2 = 0.5))
  expect_gte(fit$est, min(fit$rho.plot$theta))
  expect_lte(fit$est, max(fit$rho.plot$theta))
})

test_that("grid_rho_search survives values entirely outside [-2, 2]", {
  expect_silent(grid_rho_search(c(3.1, 3.2, 3.3), function(th) rep(0.01, 3), 1.96))
  expect_silent(grid_rho_search(c(-5, -5.1, -4.9), function(th) rep(0.01, 3), 1.96))
})

test_that("DIF and DASB tables use the same row names", {
  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  expect_identical(rownames(mod$dif_human), rownames(mod$scoring_bias))
  expect_identical(rownames(mod$dif_human), rownames(mod$ai_effect))
})

# ---------------------------------------------------------------- fit_aidif

test_that("fit_aidif works with and without an AI condition", {
  mod_h <- fit_aidif(make_mle(d_shift_g2 = 0.5))
  expect_null(mod_h$ai_fit)
  expect_null(mod_h$scoring_bias)

  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  expect_s3_class(mod$ai_fit, "rdif")
  expect_s3_class(mod$scoring_bias, "dasb")
})

test_that("fit_aidif rejects incompatible inputs and bad arguments", {
  expect_error(fit_aidif(make_mle(n_items = 4), make_mle(n_items = 6)), "items")
  expect_error(fit_aidif(make_mle(), make_mle(), alpha = 1.5), "alpha")
  expect_error(fit_aidif(make_mle(), make_mle(), metric = "nope"))
})

# ------------------------------------------------------------- S3 behaviour

test_that("summary returns an object rather than NULL (0.2.0 change)", {
  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  s <- summary(mod)
  expect_s3_class(s, "summary.aidif")
  expect_false(is.null(s))
  expect_equal(s$n_items, 4L)
  expect_false(is.null(s$human$n_iter))
  expect_false(is.null(s$human$multiple_solutions))
})

test_that("summary can recompute the multiplicity adjustment", {
  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  s <- summary(mod, adjust = "holm")
  expect_true("p_adj" %in% names(s$dasb))
  expect_true(all(s$dasb$p_adj >= s$dasb$p_val - 1e-12))
})

test_that("summary prints the lines 0.1.0 silently dropped", {
  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  txt <- paste(capture.output(print(summary(mod))), collapse = "\n")
  expect_match(txt, "Items: 4")
  expect_match(txt, "Iterations:")
  expect_match(txt, "Multiple solutions:")
})

test_that("all four plot types run, including rho", {
  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  pdf(NULL)
  on.exit(dev.off())
  for (ty in c("dif_forest", "dasb", "weights", "rho")) {
    expect_silent(plot(mod, type = ty))
  }
})

test_that("anchor_weights and ai_effect_summary accept the fitted object", {
  mod <- fit_aidif(make_mle(d_shift_g2 = 0.5), make_mle(d_shift_g2 = 0.8))
  w <- anchor_weights(mod)
  expect_named(w, c("human_weight", "ai_weight"))
  expect_true(all(w$human_weight >= 0 & w$human_weight <= 1))

  eff <- ai_effect_summary(mod)
  expect_true(all(eff$status %in% c("stable_clean", "stable_dif", "introduced",
                                    "masked", "new_direction")))
})

# --------------------------------------------------------------- simulation

test_that("simulate_aidif_data output feeds fit_aidif", {
  dat <- simulate_aidif_data(n_items = 5, dif_items = 1, dasb_items = 3,
                             seed = 7)
  expect_named(dat, c("human", "ai"))
  expect_no_error(fit_aidif(dat$human, dat$ai))
})

test_that("make_aidif_eg recovers its planted DASB at item 3", {
  eg  <- make_aidif_eg()
  sb  <- scoring_bias_test(eg$human, eg$ai, metric = "linked")
  expect_equal(which.max(abs(sb$DASB)), 3L)
})

# ----------------------------------------------------------------- adapters

test_that("as_aidif rejects objects it cannot convert", {
  expect_error(as_aidif(list()), "No as_aidif\\(\\) method")
})

# The mirt adapter is the piece most likely to break on a mirt update, and a
# wrong covariance slice would not raise an error -- it would just report the
# wrong standard errors. These tests check the values against mirt's own
# output rather than only checking that the call succeeds.

test_that("as_aidif matches mirt's estimates and SEs (polytomous)", {
  skip_if_not_installed("mirt")
  grp <- rep(c("a", "b"), each = 196)
  m <- mirt::multipleGroup(mirt::Science, 1, group = grp, SE = TRUE,
                           verbose = FALSE)
  x <- as_aidif(m, groups = c("a", "b"))

  expect_named(x, c("par.names", "est", "var.cov"))
  expect_named(x$est$group.1, c("a1", "d1", "d2", "d3"))
  expect_equal(nrow(x$est$group.1), 4L)
  expect_equal(dim(x$var.cov$group.1), c(16L, 16L))

  # Estimates must equal mirt's coefficients.
  cf <- mirt::coef(m, simplify = TRUE)$a$items[, c("a1", "d1", "d2", "d3")]
  expect_equal(as.matrix(x$est$group.1), as.matrix(cf),
               ignore_attr = TRUE, tolerance = 1e-8)

  # SEs must equal mirt's, which proves the covariance slice is at the right
  # positions rather than merely the right size.
  pv <- mirt::mod2values(m)
  er <- pv[pv$est, ]
  er <- er[order(er$parnum), ]
  se_mirt <- sqrt(diag(mirt::extract.mirt(m, "vcov")))
  key <- paste(er$group, er$item, er$name, sep = ".")
  want <- paste("a", rep(rownames(cf), each = 4),
                rep(c("a1", "d1", "d2", "d3"), times = 4), sep = ".")
  expect_equal(unname(sqrt(diag(x$var.cov$group.1))),
               unname(se_mirt[match(want, key)]), tolerance = 1e-10)
})

test_that("as_aidif renames the 2PL intercept d to d1 (dichotomous)", {
  skip_if_not_installed("mirt")
  d <- mirt::expand.table(mirt::LSAT7)
  g <- rep(c("a", "b"), length.out = nrow(d))
  m <- mirt::multipleGroup(d, 1, group = g, itemtype = "2PL", SE = TRUE,
                           verbose = FALSE)
  x <- as_aidif(m, groups = c("a", "b"))

  expect_named(x$est$group.1, c("a1", "d1"))
  expect_equal(dim(x$var.cov$group.1), c(2L * nrow(x$est$group.1)) * c(1L, 1L))
  expect_silent(check_aidif_mle(x, label = "x"))
})

test_that("as_aidif output feeds fit_aidif end to end", {
  skip_if_not_installed("mirt")
  d <- mirt::expand.table(mirt::LSAT7)
  g <- rep(c("a", "b"), length.out = nrow(d))
  m <- mirt::multipleGroup(d, 1, group = g, itemtype = "2PL", SE = TRUE,
                           verbose = FALSE)
  x <- as_aidif(m, groups = c("a", "b"))

  # Same calibration in both roles: every DASB must be exactly zero.
  fit <- fit_aidif(x, ai_mle = NULL)
  expect_s3_class(fit, "aidif")

  fit2 <- fit_aidif(x, x)
  expect_true(all(abs(fit2$scoring_bias$DASB) < 1e-12))
})

test_that("as_aidif errors on unknown groups and on too many", {
  skip_if_not_installed("mirt")
  grp <- rep(c("a", "b"), each = 196)
  m <- mirt::multipleGroup(mirt::Science, 1, group = grp, SE = TRUE,
                           verbose = FALSE)
  expect_error(as_aidif(m, groups = c("a", "zzz")), "not found")
  expect_error(as_aidif(m, groups = "a"), "exactly two")
})