# ---- .adghNodes1 -------------------------------------------------------------

test_that(".adghNodes1(1) returns trivial node", {
  g <- admixr2:::.adghNodes1(1L)
  expect_equal(g$x, 0)
  expect_equal(g$w, 1)
})

test_that(".adghNodes1(m): weights sum to 1 (probabilists' form)", {
  for (m in c(2L, 3L, 5L, 7L)) {
    g <- admixr2:::.adghNodes1(m)
    expect_equal(sum(g$w), 1, tolerance = 1e-12,
                 label = sprintf("sum(w), m=%d", m))
  }
})

test_that(".adghNodes1(m): weights are non-negative", {
  for (m in c(2L, 3L, 5L, 7L)) {
    g <- admixr2:::.adghNodes1(m)
    expect_true(all(g$w >= 0), label = sprintf("all(w>=0), m=%d", m))
  }
})

test_that(".adghNodes1(m): E[x^2] = 1 under N(0,1) (probabilists' GH)", {
  for (m in c(2L, 3L, 5L, 7L)) {
    g <- admixr2:::.adghNodes1(m)
    expect_equal(sum(g$w * g$x^2), 1, tolerance = 1e-10,
                 label = sprintf("E[x^2]=1, m=%d", m))
  }
})

test_that(".adghNodes1(m): E[x] = 0 (symmetry)", {
  for (m in c(2L, 3L, 5L, 7L)) {
    g <- admixr2:::.adghNodes1(m)
    expect_equal(sum(g$w * g$x), 0, tolerance = 1e-12,
                 label = sprintf("E[x]=0, m=%d", m))
  }
})

test_that(".adghNodes1(m): returns m nodes", {
  for (m in c(1L, 3L, 5L)) {
    g <- admixr2:::.adghNodes1(m)
    expect_length(g$x, m)
    expect_length(g$w, m)
  }
})

test_that(".adghNodes1(m): m < 1 errors", {
  expect_error(admixr2:::.adghNodes1(0L))
})

# ---- .adghNodeGrid -----------------------------------------------------------

test_that(".adghNodeGrid(m, 0): degenerate grid (no etas)", {
  g <- admixr2:::.adghNodeGrid(5L, 0L)
  expect_equal(nrow(g$X), 1L)
  expect_equal(ncol(g$X), 0L)
  expect_equal(g$W, 1)
})

test_that(".adghNodeGrid(m, 1): m nodes, weights sum to 1", {
  g <- admixr2:::.adghNodeGrid(5L, 1L)
  expect_equal(nrow(g$X), 5L)
  expect_equal(ncol(g$X), 1L)
  expect_equal(sum(g$W), 1, tolerance = 1e-12)
})

test_that(".adghNodeGrid(m, 2): m^2 nodes, weights sum to 1", {
  g <- admixr2:::.adghNodeGrid(3L, 2L)
  expect_equal(nrow(g$X), 9L)
  expect_equal(ncol(g$X), 2L)
  expect_equal(sum(g$W), 1, tolerance = 1e-12)
})

test_that(".adghNodeGrid(m, n): E[x_j^2] = 1 for each dimension", {
  g <- admixr2:::.adghNodeGrid(5L, 2L)
  for (j in 1:2)
    expect_equal(sum(g$W * g$X[, j]^2), 1, tolerance = 1e-10,
                 label = sprintf("E[x_%d^2]=1", j))
})

test_that(".adghNodeGrid(m, n): marginal weights are correct (colSums via product)", {
  m <- 3L; d <- 2L
  g <- admixr2:::.adghNodeGrid(m, d)
  g1 <- admixr2:::.adghNodes1(m)
  # marginal weight for first dim: sum over all second-dim combinations
  # = product of marginal weights * sum of second-dim weights
  expect_equal(length(g$W), m^d)
  expect_true(all(g$W > 0))
})

# ---- adghControl -------------------------------------------------------------

test_that("adghControl() returns correct class and key defaults", {
  ctl <- adghControl()
  expect_s3_class(ctl, "adghControl")
  expect_equal(ctl$n_nodes,    5L)
  expect_equal(ctl$grad,       "analytical")
  expect_equal(ctl$maxeval,    500L)
  expect_equal(ctl$seed,       12345L)
  expect_equal(ctl$cores,      as.integer(rxode2::rxCores()))
  expect_equal(ctl$n_restarts, 1L)
  expect_equal(ctl$workers,    1L)
  # the DEFAULT: "r,s" reduces to "r" under correct specification, so it is the
  # conservative choice rather than the aggressive one, and it degrades to "r"
  # (reporting "r") wherever it cannot be built.
  expect_equal(ctl$covMethod,  "r,s")
  expect_equal(ctl$returnAdmr, FALSE)
})

test_that("adghControl(): default grad='analytical' coerces BOBYQA to LBFGS", {
  ctl <- adghControl()
  expect_equal(ctl$algorithm, "NLOPT_LD_LBFGS")
})

test_that("adghControl(): grad='fd' coerces BOBYQA to LBFGS", {
  ctl <- adghControl(grad = "fd")
  expect_equal(ctl$algorithm, "NLOPT_LD_LBFGS")
})

test_that("adghControl(): grad='fd' coerces BOBYQA to LBFGS", {
  ctl <- adghControl(grad = "fd")
  expect_equal(ctl$algorithm, "NLOPT_LD_LBFGS")
})

test_that("adghControl(): grad='none' keeps BOBYQA", {
  ctl <- adghControl(grad = "none")
  expect_equal(ctl$algorithm, "NLOPT_LN_BOBYQA")
})

test_that("adghControl(): grad != 'none' + explicit non-BOBYQA algorithm kept", {
  ctl <- adghControl(grad = "fd", algorithm = "NLOPT_LD_SLSQP")
  expect_equal(ctl$algorithm, "NLOPT_LD_SLSQP")
})

test_that("adghControl(): an EXPLICIT derivative-free algorithm turns the gradient off", {
  # The combination that changed in 0.4.1, and the one this file did not pin.
  #
  # adgh used to special-case exactly "NLOPT_LN_BOBYQA" and upgrade it to LBFGS
  # whenever a gradient was requested -- necessarily, because BOBYQA was
  # adghControl's own DEFAULT, so naming it was indistinguishable from leaving it
  # alone. With the default now NULL that ambiguity is gone, and naming a
  # derivative-free algorithm means what it says, exactly as it does for the
  # other three controls.
  #
  # Pinned because it is a silent change for anyone who explicitly restated the
  # old default: they get a derivative-free fit where they had a quasi-Newton one.
  for (g in c("analytical", "fd")) {
    ctl <- suppressMessages(adghControl(grad = g, algorithm = "NLOPT_LN_BOBYQA"))
    expect_equal(ctl$algorithm, "NLOPT_LN_BOBYQA", info = g)
    expect_equal(ctl$grad, "none", info = g)
  }
  # ... and it is not silent: the reconciliation says so.
  expect_message(adghControl(grad = "analytical", algorithm = "NLOPT_LN_BOBYQA"))
  # The documented escape hatch -- drop the argument and NULL picks LBFGS.
  expect_equal(adghControl(grad = "analytical")$algorithm, "NLOPT_LD_LBFGS")
  expect_equal(adghControl(grad = "analytical")$grad, "analytical")
})

test_that("adghControl(): internal n_sim=1L for .admRunRestarts() compat", {
  ctl <- adghControl()
  expect_type(ctl$n_sim, "integer")
  expect_equal(ctl$n_sim, 1L)
})

test_that("adghControl(): internal sampling='sobol' for .admRunRestarts() compat", {
  ctl <- adghControl()
  expect_equal(ctl$sampling, "sobol")
})

test_that("adghControl(): n_nodes stored as integer", {
  ctl <- adghControl(n_nodes = 3)
  expect_type(ctl$n_nodes, "integer")
  expect_equal(ctl$n_nodes, 3L)
})

test_that("adghControl(): maxeval stored as integer", {
  ctl <- adghControl(maxeval = 200)
  expect_type(ctl$maxeval, "integer")
  expect_equal(ctl$maxeval, 200L)
})

test_that("adghControl(): cov_h_outer default is eps^(1/4)", {
  ctl <- adghControl()
  expect_equal(ctl$cov_h_outer, .Machine$double.eps^(1/4), tolerance = 1e-15)
})

test_that("adghControl(): covMethod 'none' accepted", {
  ctl <- adghControl(covMethod = "none")
  expect_equal(ctl$covMethod, "none")
})

test_that("adghControl(): returnAdmr stored correctly", {
  ctl <- adghControl(returnAdmr = TRUE)
  expect_true(ctl$returnAdmr)
})

test_that("adghControl(): unknown argument errors", {
  expect_error(adghControl(nonexistent_arg = 42), regexp = "unused argument")
})

test_that("adghControl(): n_nodes=0 errors via checkmate", {
  expect_error(adghControl(n_nodes = 0L), regexp = "n_nodes")
})

test_that("adghControl(): cores=0 errors via checkmate", {
  expect_error(adghControl(cores = 0L), regexp = "cores")
})

test_that("adghControl(): maxeval=0 errors via checkmate", {
  expect_error(adghControl(maxeval = 0L), regexp = "maxeval")
})

test_that("adghControl(): invalid grad errors", {
  expect_error(adghControl(grad = "newton"))
})

test_that("adghControl(): invalid covMethod errors", {
  expect_error(adghControl(covMethod = "hessian"))
})

test_that("adghControl(): grad_h, grad_bounds, cov_h stored", {
  ctl <- adghControl(grad_h = 1e-5, grad_bounds = 3, cov_h = 5e-4)
  expect_equal(ctl$grad_h,      1e-5)
  expect_equal(ctl$grad_bounds, 3)
  expect_equal(ctl$cov_h,       5e-4)
})

# ---- S3 hooks ----------------------------------------------------------------

test_that("getValidNlmixrCtl.adgh passes through adghControl", {
  ctl <- adghControl(maxeval = 99L)
  out <- admixr2:::getValidNlmixrCtl.adgh(ctl)
  expect_s3_class(out, "adghControl")
  expect_equal(out$maxeval, 99L)
})

test_that("getValidNlmixrCtl.adgh returns default for NULL", {
  out <- admixr2:::getValidNlmixrCtl.adgh(NULL)
  expect_s3_class(out, "adghControl")
})

test_that("getValidNlmixrCtl.adgh converts plain list with n_nodes", {
  out <- admixr2:::getValidNlmixrCtl.adgh(list(n_nodes = 3L, maxeval = 77L))
  expect_s3_class(out, "adghControl")
  expect_equal(out$n_nodes, 3L)
})

test_that("getValidNlmixrCtl.adgh extracts nested adghControl", {
  ctl <- adghControl(maxeval = 66L)
  out <- admixr2:::getValidNlmixrCtl.adgh(list(ctl))
  expect_s3_class(out, "adghControl")
  expect_equal(out$maxeval, 66L)
})

test_that("getValidNlmixrCtl.adgh builds from nested list with studies key", {
  out <- admixr2:::getValidNlmixrCtl.adgh(list(list(studies = list(), maxeval = 55L)))
  expect_s3_class(out, "adghControl")
})

test_that("nmObjHandleControlObject.adghControl: assigns control to env", {
  ctl <- adghControl(maxeval = 22L)
  e   <- new.env(parent = emptyenv())
  admixr2:::nmObjHandleControlObject.adghControl(ctl, e)
  expect_true(exists("adghControl", envir = e, inherits = FALSE))
  expect_equal(get("adghControl", envir = e)$maxeval, 22L)
})

test_that("nmObjGetControl.adgh: retrieves control from env", {
  ctl <- adghControl(maxeval = 33L)
  e   <- new.env(parent = emptyenv())
  e$adghControl <- ctl
  out <- admixr2:::nmObjGetControl.adgh(list(e))
  expect_s3_class(out, "adghControl")
  expect_equal(out$maxeval, 33L)
})

test_that("nmObjGetControl.adgh: falls back to 'control' slot", {
  ctl <- adghControl(maxeval = 44L)
  e   <- new.env(parent = emptyenv())
  e$control <- ctl
  out <- admixr2:::nmObjGetControl.adgh(list(e))
  expect_s3_class(out, "adghControl")
  expect_equal(out$maxeval, 44L)
})

test_that("nmObjGetControl.adgh: errors when no control found", {
  e <- new.env(parent = emptyenv())
  expect_error(admixr2:::nmObjGetControl.adgh(list(e)),
               regexp = "cannot find adgh control")
})

# ---- grid memo ---------------------------------------------------------------

test_that("the tensor grid memo answers by key, not by whatever it cached last", {
  a  <- admixr2:::.adghNodeGrid(4L, 2L)
  b  <- admixr2:::.adghNodeGrid(6L, 2L)
  a2 <- admixr2:::.adghNodeGrid(4L, 2L)          # after b displaced a
  expect_equal(a, a2)
  expect_equal(nrow(a$X), 16L)
  expect_equal(nrow(b$X), 36L)
})

test_that("a uniform nv is the ordinary tensor grid, not a second copy of it", {
  expect_identical(admixr2:::.admNodeGridNv(c(3L, 3L)),
                   admixr2:::.adghNodeGrid(3L, 2L))
  g <- admixr2:::.admNodeGridNv(c(3L, 5L))
  expect_equal(nrow(g$X), 15L)
  expect_equal(sum(g$W), 1, tolerance = 1e-12)
})

# ---- batched solves ----------------------------------------------------------

test_that("grouping studies does not move the objective, across doses", {
  skip_on_cran(); skip_if_not_installed("rxode2")
  # The batched path must agree with one-call-per-study EXACTLY, including when
  # the group spans different dosing schedules -- those stack with offset ids,
  # and getting that wrong is invisible in the answer's shape.
  mod <- function() {
    ini({ tcl <- log(3.2); tv <- log(21); bwt <- .6
          eta.cl ~ .09; add.err <- .6 })
    model({ cl <- exp(tcl + eta.cl)*(WT/70)^bwt
            v <- exp(tv); cp <- linCmt(); cp ~ add(add.err) }) }
  ui  <- suppressMessages(rxode2::rxode2(mod))
  ov  <- admixr2:::.admOutputVar(ui)
  pop <- admPopulation(WT = c(mean = 76, sd = 15))
  mk  <- function(n, d, tms) admStudy(model = mod, population = pop, n = n,
                                      dose = d, times = tms, strata_nodes = 3L)
  g <- suppressMessages(suppressWarnings(admixr2:::.admMaterialise(
    admStudies(a = mk(200L, 100, c(1, 4, 8)),
               b = mk(160L, 250, c(2, 6, 12)),      # different dose AND times
               c = mk(240L, 175, c(1, 4, 8))))))
  ctl <- adghControl(studies = g, print = 0L, covMethod = "none")
  p   <- admixr2:::.admDriverPinfo(ui, ctl)
  u   <- admixr2:::.admDriverUnits(g, ui, ov)
  stu <- suppressMessages(admixr2:::.admCheckCovariates(ui, p, u$studies))
  rx  <- admixr2:::.admLoadModel(ui)
  gr  <- admixr2:::.adghNodeGrid(p$n_nodes %||% 5L, p$n_eta)
  # more than one event key present, so the stacking path is the one under test
  expect_gt(length(unique(vapply(stu, function(z) z$ev_key %||% "", ""))), 1L)
  # ...and they really do land in ONE group, which is the point of stacking.
  expect_equal(length(unique(admixr2:::.adghSolveGroupKeys(stu, ov))), 1L)

  # Accumulated the way .adghNLL() accumulates: a running total, not sum() over
  # a vector, or the two differ in the last bit purely by summation order (1e-13
  # on ~920) and say nothing about the batching.
  per_study <- function(pv) {
    pars <- admixr2:::.admUnpack(pv, p)
    tot <- 0
    for (s in stu) {
      m <- admixr2:::.adghMoments(pars, p, s, rx, s$output %||% ov, gr, 1L)
      tot <- tot + if (identical(s$method, "var"))
        admixr2:::nll_var_cpp(s$E, s$v_diag, m$E, diag(m$V), s$n)
      else admixr2:::nll_cov_cpp(s$E, s$V, m$E, m$V, s$n)
    }
    tot
  }
  p0 <- admixr2:::.admBuildOptVec(p)$p0
  set.seed(7)
  for (k in 0:2) {
    pv <- if (k == 0L) p0 else p0 + rnorm(length(p0), 0, 0.05)
    expect_equal(admixr2:::.adghNLL(pv, p, stu, rx, ov, gr, 1L), per_study(pv),
                 tolerance = 0)
  }
})

test_that(".admSimulateMany reproduces .admSimulate, same events or not", {
  skip_on_cran(); skip_if_not_installed("rxode2")
  mod <- function() {
    ini({ tcl <- log(3.2); tv <- log(21); eta.cl ~ .09; add.err <- .6 })
    model({ cl <- exp(tcl + eta.cl); v <- exp(tv)
            cp <- linCmt(); cp ~ add(add.err) }) }
  ui <- suppressMessages(rxode2::rxode2(mod))
  ov <- admixr2:::.admOutputVar(ui); rx <- admixr2:::.admLoadModel(ui)
  mk <- function(d, tms) list(E = rep(1, length(tms)), V = diag(length(tms)),
                              n = 100L, times = tms, ev = rxode2::et(amt = d))
  for (lab in c("shared", "mixed")) {
    st <- if (lab == "shared")
      list(a = mk(100, c(1, 4, 8)), b = mk(100, c(1, 4, 8)))
    else
      list(a = mk(100, c(1, 4, 8)), b = mk(250, c(2, 6, 12)))
    u  <- admixr2:::.admDriverUnits(st, ui, ov)
    p  <- admixr2:::.admDriverPinfo(ui, adghControl(studies = st, print = 0L))
    n  <- 60L
    z  <- admixr2:::.admMakeZ(n, p, 1L, "sobol")
    pl <- admixr2:::.admMakeParamsList(n, p, 1L)
    e  <- z[[1L]] %*% t(diag(sqrt(.09), 1L))
    colnames(e) <- p$eta_col_names
    th <- c(tcl = log(3.2), tv = log(21))
    one <- lapply(u$studies, function(s)
      admixr2:::.admSimulate(rx, th, "add.err", e, s, ov, pl[[1L]], 1L))
    many <- admixr2:::.admSimulateMany(rx, th, "add.err",
              rep(list(e), length(u$studies)), u$studies, ov,
              rep(list(pl[[1L]]), length(u$studies)), 1L)
    expect_equal(length(many), length(one), info = lab)
    for (k in seq_along(one))
      expect_equal(many[[k]], one[[k]], tolerance = 0, info = lab)
  }
})