tests/testthat/test_LINAD_differential.R

In rtemis: Schema-Driven Machine Learning and Visualization

# test_LINAD_differential.R
# ::rtemis::
# 2026- EDG rtemis.org

# The engine's fast paths checked against slow implementations of the same
# quantity. The fast paths exist because the obvious computation is too
# expensive to run at every candidate split; these tests run the obvious
# computation at a size where it is affordable, and require the two to agree.
#
# What this catches that the reduction tests cannot: the reductions pin the
# degenerate corner -- constant leaves, a hard partition, no shrinkage -- while
# the configuration the algorithm is actually used in has no external reference.
# Here the reference is written from the definition, so it covers linear leaves,
# soft weighting and the exhaustive search directly.

# %% Reference implementations ----

#' Every candidate split scored from the definition
#'
#' No incremental sufficient statistics, no subtraction of node totals, no
#' algebraic loss identity: each side's weights are formed explicitly over every
#' row, each side's model is fitted from those weights, and the loss is the
#' weighted sum of squared residuals of the fitted values.
#'
#' Each side is fitted by `linad_solve()` -- the function the *commit* uses, not
#' the one the search uses -- so this reference states the search's contract
#' rather than restating its implementation: a candidate must be scored by the
#' model the node will actually receive. The active rows are the node's own,
#' resolved once as `linad_propose()` resolves them, which is the rule both
#' sides are required to share.
reference_sweep <- function(state, r, w, member, derivatives = NULL) {
  xm <- state[["xm"]]
  context <- state[["context"]]
  gamma <- state[["gamma"]]
  n <- nrow(xm)
  min_child <- rtemis:::linad_min_child_cases(state)
  total_members <- sum(member)
  active <- rtemis:::linad_active_rows(w)

  side_loss <- function(weights) {
    if (sum(weights[active]) <= 0) {
      return(Inf)
    }
    fit <- rtemis:::linad_solve(
      xm,
      r,
      weights,
      active,
      state[["node_model"]],
      state[["lambda"]],
      state[["alpha"]] %||% 1,
      state[["nvmax"]],
      forward_stop = state[["forward_stop"]],
      derivatives = derivatives,
      type = state[["type"]] %||% "Regression",
      max_step = state[["line_search_max"]] %||% 1000,
      constant_rule = state[["constant_rule"]] %||% "closed_form",
      node_test = state[["node_test"]],
      allowed = state[["adaptive_columns"]]
    )
    if (is.null(fit) || !all(is.finite(fit[["coefficients"]]))) {
      return(Inf)
    }
    fitted <- drop(xm %*% fit[["coefficients"]])
    sum(weights[active] * (r[active] - fitted[active])^2)
  }

  loss_of <- function(goes_left) {
    if (
      sum(member & goes_left) < min_child ||
        (total_members - sum(member & goes_left)) < min_child
    ) {
      return(Inf)
    }
    left <- w
    left[!goes_left] <- left[!goes_left] * gamma
    right <- w
    right[goes_left] <- right[goes_left] * gamma
    side_loss(left) + side_loss(right)
  }

  best <- list(
    loss = Inf,
    feature = NA_character_,
    levels = NULL,
    value = NA_real_
  )
  for (j in seq_along(context[["numeric_names"]])) {
    order_index <- context[["numeric_order"]][[j]]
    breaks <- context[["numeric_breaks"]][[j]]
    sorted <- context[["numeric_matrix"]][order_index, j]
    if (length(breaks) == 0L) {
      next
    }
    # Which cut points the search considers is a separate question from how it
    # scores them, and this reference is about the scoring.
    wanted <- rtemis:::linad_cut_positions(
      breaks,
      state[["n_cuts"]],
      sorted,
      state[["split_bin_type"]]
    )
    wanted <- wanted[wanted >= 1L & wanted < n]
    for (position in wanted) {
      goes_left <- logical(n)
      goes_left[order_index[seq_len(position)]] <- TRUE
      loss <- loss_of(goes_left)
      if (loss < best[["loss"]]) {
        best <- list(
          loss = loss,
          feature = context[["numeric_names"]][[j]],
          levels = NULL,
          value = (sorted[[position]] + sorted[[position + 1L]]) / 2
        )
      }
    }
  }
  for (j in seq_along(context[["factor_names"]])) {
    codes <- context[["factor_codes"]][[j]]
    levels_j <- context[["factor_levels"]][[j]]
    present <- sort(unique(codes))
    k <- length(present)
    if (k < 2L) {
      next
    }
    # Every partition, enumerated without regard to any ordering of the levels.
    for (mask in seq_len(2^(k - 1L) - 1L)) {
      chosen <- present[bitwAnd(mask, bitwShiftL(1L, seq_len(k) - 1L)) > 0L]
      loss <- loss_of(codes %in% chosen)
      if (loss < best[["loss"]]) {
        best <- list(
          loss = loss,
          feature = context[["factor_names"]][[j]],
          levels = levels_j[chosen],
          value = NA_real_
        )
      }
    }
  }
  best
} # /reference_sweep


#' A fit state carrying only what the split searches read
reference_state <- function(x, xm, settings = list()) {
  defaults <- list(
    context = rtemis:::linad_context(x),
    xm = xm,
    n = nrow(xm),
    gamma = 0,
    n_cuts = 20L,
    node_model = "ridge",
    lambda = 0.1,
    alpha = 1,
    nvmax = 3L,
    forward_stop = "none",
    node_test = "none",
    type = "Regression",
    line_search_max = 1000,
    constant_rule = "closed_form",
    adaptive_columns = NULL,
    split_search = "exhaustive",
    split_criterion = "mean",
    split_bin_type = "frequency",
    min_cases_leaf = 5L,
    min_cases_node_model = 5L
  )
  utils::modifyList(defaults, settings)
} # /reference_state


# %% The exhaustive search against the definition ----
test_that("linad_sweep() scores every candidate as the definition does", {
  set.seed(2026)
  n <- 220L
  x <- data.frame(
    a = rnorm(n),
    b = rnorm(n),
    g = factor(sample(c("p", "q", "r"), n, replace = TRUE)),
    h = factor(sample(c("u", "v"), n, replace = TRUE))
  )
  xm <- cbind(
    1,
    x[["a"]],
    x[["b"]],
    as.numeric(x[["g"]] == "q"),
    as.numeric(x[["g"]] == "r"),
    as.numeric(x[["h"]] == "v")
  )
  r <- 3 *
    x[["a"]] *
    ifelse(x[["b"]] > 0, 1, -1) +
    ifelse(x[["g"]] == "q", 2, 0) +
    rnorm(n, 0, 0.5)
  w <- runif(n, 0.5, 2)
  member <- rep(TRUE, n)

  grid <- expand.grid(
    node_model = c("constant", "ridge", "forward"),
    gamma = c(0, 0.2),
    stringsAsFactors = FALSE
  )
  for (row in seq_len(nrow(grid))) {
    settings <- list(
      node_model = grid[["node_model"]][[row]],
      gamma = grid[["gamma"]][[row]]
    )
    state <- reference_state(x, xm, settings)
    fast <- rtemis:::linad_sweep(state, r, w, member)
    slow <- reference_sweep(state, r, w, member)
    label <- paste(settings[["node_model"]], settings[["gamma"]])
    expect_equal(fast[["feature"]], slow[["feature"]], info = label)
    # `gain` is the negated loss, so the two are comparable directly.
    expect_equal(
      -fast[["gain"]],
      slow[["loss"]],
      tolerance = 1e-8,
      info = label
    )
    if (identical(fast[["kind"]], "factor")) {
      # Either side of a partition names it.
      expect_true(
        setequal(fast[["levels"]], slow[["levels"]]) ||
          setequal(
            fast[["levels"]],
            setdiff(
              state[["context"]][["factor_levels"]][[fast[["column"]]]],
              slow[["levels"]]
            )
          ),
        info = label
      )
    } else {
      expect_equal(
        fast[["value"]],
        slow[["value"]],
        tolerance = 1e-8,
        info = label
      )
    }
  }
})


test_that("linad_sweep() honors the child-size floor as the definition does", {
  # A floor that the search and the reference read differently would show as a
  # disagreement only on candidates near it, which is where they are placed.
  set.seed(11)
  n <- 120L
  x <- data.frame(a = rnorm(n), b = rnorm(n))
  xm <- cbind(1, x[["a"]], x[["b"]])
  r <- x[["a"]] + rnorm(n, 0, 0.5)
  w <- rep(1, n)
  member <- rep(TRUE, n)
  member[sample.int(n, 40L)] <- FALSE
  for (floor_size in c(5L, 20L, 45L)) {
    state <- reference_state(
      x,
      xm,
      list(min_cases_leaf = floor_size, min_cases_node_model = floor_size)
    )
    fast <- rtemis:::linad_sweep(state, r, w, member)
    slow <- reference_sweep(state, r, w, member)
    if (is.null(fast)) {
      expect_false(is.finite(slow[["loss"]]), info = floor_size)
    } else {
      expect_equal(
        -fast[["gain"]],
        slow[["loss"]],
        tolerance = 1e-8,
        info = floor_size
      )
    }
  }
})


test_that("linad_sweep() agrees with the definition on randomized problems", {
  # Twelve draws varying shape, weights, membership and hyperparameters, so a
  # defect that needs a particular configuration has somewhere to show.
  set.seed(4)
  for (rep in seq_len(12L)) {
    n <- sample(80:200, 1L)
    n_numeric <- sample(1:3, 1L)
    n_levels <- sample(2:4, 1L)
    x <- as.data.frame(matrix(rnorm(n * n_numeric), n, n_numeric))
    names(x) <- paste0("v", seq_len(n_numeric))
    x[["f"]] <- factor(sample(letters[seq_len(n_levels)], n, replace = TRUE))
    xm <- cbind(
      1,
      as.matrix(x[seq_len(n_numeric)]),
      stats::model.matrix(~f, x)[, -1L, drop = FALSE]
    )
    r <- rnorm(n) + as.numeric(x[["f"]])
    w <- if (rep %% 2L == 0L) rep(1, n) else runif(n, 0.2, 3)
    member <- rep(TRUE, n)
    if (rep %% 3L == 0L) {
      member[sample.int(n, n %/% 4L)] <- FALSE
    }
    state <- reference_state(
      x,
      xm,
      list(
        node_model = sample(c("constant", "ridge", "forward"), 1L),
        gamma = sample(c(0, 0.1, 0.5), 1L),
        lambda = sample(c(0.01, 0.5), 1L),
        nvmax = sample(1:3, 1L),
        min_cases_leaf = sample(3:10, 1L),
        min_cases_node_model = sample(3:10, 1L)
      )
    )
    fast <- rtemis:::linad_sweep(state, r, w, member)
    slow <- reference_sweep(state, r, w, member)
    expect_equal(
      -fast[["gain"]],
      slow[["loss"]],
      tolerance = 1e-8,
      info = paste("replicate", rep)
    )
    expect_equal(
      fast[["feature"]],
      slow[["feature"]],
      info = paste("replicate", rep)
    )
  }
})


# %% Classification, where the search and the commit once diverged ----
test_that("linad_sweep() scores a classification candidate as the commit fits it", {
  # The committed child's constant is a Newton step over the weighted
  # derivatives; the search once had only `X'WX` and `X'Wr` and used the
  # weighted least-squares level instead, so it ranked splits by a model that
  # was never built. Every configuration here is checked against a reference
  # that fits both children with `linad_solve()`.
  set.seed(2027)
  n <- 200L
  x <- data.frame(
    a = rnorm(n),
    b = rnorm(n),
    g = factor(sample(c("p", "q", "r"), n, replace = TRUE))
  )
  xm <- cbind(
    1,
    x[["a"]],
    x[["b"]],
    as.numeric(x[["g"]] == "q"),
    as.numeric(x[["g"]] == "r")
  )
  labels <- ifelse(
    2 * x[["a"]] * ifelse(x[["b"]] > 0, 1, -1) + rnorm(n) > 0,
    1,
    -1
  )
  grid <- expand.grid(
    node_model = c("constant", "ridge", "forward", "elasticnet"),
    gamma = c(0, 0.3),
    constant_rule = c("closed_form", "least_squares"),
    stringsAsFactors = FALSE
  )
  for (row in seq_len(nrow(grid))) {
    f <- rnorm(n, 0, 0.4)
    derivatives <- rtemis:::linad_gradient(labels, f, "Classification")
    r <- -derivatives[["g"]]
    w <- runif(n, 0.5, 2)
    member <- rep(TRUE, n)
    state <- reference_state(
      x,
      xm,
      list(
        node_model = grid[["node_model"]][[row]],
        gamma = grid[["gamma"]][[row]],
        constant_rule = grid[["constant_rule"]][[row]],
        alpha = 0.5,
        type = "Classification"
      )
    )
    fast <- rtemis:::linad_sweep(
      state,
      r,
      w,
      member,
      derivatives = derivatives,
      active = rtemis:::linad_active_rows(w)
    )
    slow <- reference_sweep(state, r, w, member, derivatives)
    label <- paste(
      grid[["node_model"]][[row]],
      grid[["gamma"]][[row]],
      grid[["constant_rule"]][[row]]
    )
    expect_equal(fast[["feature"]], slow[["feature"]], info = label)
    expect_equal(
      -fast[["gain"]],
      slow[["loss"]],
      tolerance = 1e-8,
      info = label
    )
  }
})


test_that("The search and the commit fit the same rows across the weight tolerance", {
  # `LINAD_WEIGHT_TOLERANCE` drops a case whose soft weight has decayed past a
  # fraction of the node's largest. The rule is resolved once per node and
  # shared, so weights straddling it must not move the search away from the
  # model the commit builds.
  set.seed(2028)
  n <- 160L
  x <- data.frame(a = rnorm(n), b = rnorm(n))
  xm <- cbind(1, x[["a"]], x[["b"]])
  r <- 2 * x[["a"]] * ifelse(x[["b"]] > 0, 1, -1) + rnorm(n, 0, 0.4)
  member <- rep(TRUE, n)
  tolerance <- rtemis:::LINAD_WEIGHT_TOLERANCE
  # A third of the cases sit above the threshold, a third just below it, and a
  # third far below, so both sides of the rule are exercised.
  w <- rep(1, n)
  w[seq(2L, n, by = 3L)] <- tolerance * 0.5
  w[seq(3L, n, by = 3L)] <- tolerance * 1e-4
  active <- rtemis:::linad_active_rows(w)
  expect_gt(length(active), 0L)
  expect_lt(length(active), n)
  for (node_model in c("ridge", "forward")) {
    for (gamma in c(0, 0.3)) {
      state <- reference_state(
        x,
        xm,
        list(node_model = node_model, gamma = gamma, min_cases_leaf = 5L)
      )
      fast <- rtemis:::linad_sweep(state, r, w, member, active = active)
      slow <- reference_sweep(state, r, w, member)
      label <- paste(node_model, gamma)
      expect_equal(fast[["feature"]], slow[["feature"]], info = label)
      expect_equal(
        -fast[["gain"]],
        slow[["loss"]],
        tolerance = 1e-8,
        info = label
      )
    }
  }
})


# %% Feature roles reach the search, not only the commit ----
test_that("The exhaustive search scores children under the adaptive restriction", {
  # The rule that matters most: if the search scores candidates by child models
  # that ignore the effect-scope restriction, it optimizes a model the commit
  # will never build.
  set.seed(2029)
  n <- 180L
  x <- data.frame(a = rnorm(n), b = rnorm(n), c = rnorm(n))
  xm <- cbind(1, x[["a"]], x[["b"]], x[["c"]])
  r <- 2 * x[["a"]] * ifelse(x[["b"]] > 0, 1, -1) + x[["c"]] + rnorm(n, 0, 0.4)
  w <- runif(n, 0.5, 2)
  member <- rep(TRUE, n)
  for (adaptive in list(1L, c(1L, 3L), integer(0))) {
    state <- reference_state(
      x,
      xm,
      list(node_model = "ridge", gamma = 0.2, adaptive_columns = adaptive)
    )
    fast <- rtemis:::linad_sweep(state, r, w, member)
    slow <- reference_sweep(state, r, w, member)
    label <- paste(adaptive, collapse = ",")
    expect_equal(fast[["feature"]], slow[["feature"]], info = label)
    expect_equal(
      -fast[["gain"]],
      slow[["loss"]],
      tolerance = 1e-8,
      info = label
    )
  }
})


# %% The two node-model solvers against each other ----
test_that("linad_solve() and linad_gram_solve() fit the same model", {
  # One is used by the commit and reads the data; the other is used by the
  # exhaustive search and reads sufficient statistics. They fit the same node
  # model by different routes -- one centers the design and carries the level
  # separately, the other reads the centered system off the raw Gram -- and a
  # divergence between them is a search optimizing a model that is never built.
  #
  # Both outcome types, both constant rules and every node model: the
  # classification constant is a Newton step rather than a weighted mean, and
  # the search once had no way to form it.
  set.seed(3)
  for (rep in seq_len(10L)) {
    n <- sample(60:200, 1L)
    p <- sample(2:6, 1L)
    X <- matrix(rnorm(n * p), n, p)
    xm <- cbind(1, X)
    signal <- drop(X %*% rnorm(p))
    w <- if (rep %% 2L == 0L) rep(1, n) else runif(n, 0.3, 2)
    idx <- sort(sample.int(n, max(20L, n %/% 2L)))
    lambda <- sample(c(0, 0.05, 0.5), 1L)
    alpha <- sample(c(0, 0.5, 1), 1L)
    nvmax <- sample(seq_len(p), 1L)
    for (type in c("Regression", "Classification")) {
      if (identical(type, "Classification")) {
        labels <- ifelse(signal + rnorm(n) > 0, 1, -1)
        derivatives <- rtemis:::linad_gradient(
          labels,
          rnorm(n, 0, 0.3),
          type
        )
        y <- -derivatives[["g"]]
      } else {
        derivatives <- NULL
        y <- signal + rnorm(n)
      }
      gram <- rtemis:::linad_gram(xm, y, w, idx, derivatives)
      for (node_model in c("constant", "ridge", "forward", "elasticnet")) {
        for (constant_rule in c("closed_form", "least_squares")) {
          direct <- rtemis:::linad_solve(
            xm,
            y,
            w,
            idx,
            node_model,
            lambda = lambda,
            alpha = alpha,
            nvmax = nvmax,
            forward_stop = "none",
            derivatives = derivatives,
            type = type,
            constant_rule = constant_rule
          )
          from_gram <- rtemis:::linad_gram_solve(
            gram[["G"]],
            gram[["Xty"]],
            gram[["syy"]],
            node_model,
            lambda,
            alpha,
            nvmax,
            forward_stop = "none",
            sg = gram[["sg"]],
            sh = gram[["sh"]],
            type = type,
            constant_rule = constant_rule
          )
          expect_equal(
            direct[["coefficients"]],
            from_gram,
            tolerance = 1e-8,
            info = paste(type, node_model, constant_rule, rep)
          )
        }
      }
    }
  }
})


test_that("The two solvers agree when a node may fit only some columns", {
  # Effect scope reaches the model as a column restriction, and the search and
  # the commit have to apply the identical one: a search scoring child models
  # that ignore the restriction optimizes a model the commit will never build.
  set.seed(31)
  for (rep in seq_len(8L)) {
    n <- sample(80:200, 1L)
    p <- sample(3:7, 1L)
    X <- matrix(rnorm(n * p), n, p)
    xm <- cbind(1, X)
    y <- drop(X %*% rnorm(p)) + rnorm(n)
    w <- runif(n, 0.3, 2)
    idx <- sort(sample.int(n, max(30L, n %/% 2L)))
    gram <- rtemis:::linad_gram(xm, y, w, idx)
    allowed <- sort(sample.int(p, sample(seq_len(p), 1L)))
    for (node_model in c("ridge", "forward", "elasticnet")) {
      direct <- rtemis:::linad_solve(
        xm,
        y,
        w,
        idx,
        node_model,
        lambda = 0.05,
        alpha = 0.5,
        nvmax = 3L,
        forward_stop = "none",
        allowed = allowed
      )
      from_gram <- rtemis:::linad_gram_solve(
        gram[["G"]],
        gram[["Xty"]],
        gram[["syy"]],
        node_model,
        0.05,
        0.5,
        3L,
        forward_stop = "none",
        allowed = allowed
      )
      expect_equal(
        direct[["coefficients"]],
        from_gram,
        tolerance = 1e-8,
        info = paste(node_model, rep)
      )
      # A column outside the allowed set never receives a coefficient.
      expect_true(
        all(direct[["coefficients"]][-1L][-allowed] == 0),
        info = paste(node_model, rep)
      )
    }
  }
})


test_that("linad_enet_gram() solves the objective it names", {
  # Two independent references, one per penalty. The L2 half must reproduce the
  # engine's own ridge exactly at alpha = 0 -- one `lambda` cannot mean two
  # things -- and the L1 half is checked against glmnet, whose lasso uses the
  # same parameterization.
  skip_if_not_installed("glmnet")
  set.seed(41)
  n <- 200L
  p <- 6L
  X <- matrix(rnorm(n * p), n, p)
  target <- drop(X %*% c(2, -1.5, 0, 0, 0.7, 0)) + rnorm(n)
  w <- runif(n, 0.5, 2)
  centered <- sweep(X, 2L, drop(crossprod(w, X)) / sum(w), "-")
  residual <- target - sum(w * target) / sum(w)
  gram <- rtemis:::linad_gram(centered, residual, w, NULL)
  sw <- gram[["sw"]]
  for (lambda in c(0.001, 0.05, 0.3)) {
    expect_equal(
      rtemis:::linad_enet_gram(gram[["G"]], gram[["Xty"]], lambda * sw, 0),
      rtemis:::linad_chol_solve(
        gram[["G"]],
        gram[["Xty"]],
        lambda * sw,
        intercept = FALSE
      ),
      tolerance = 1e-8,
      info = paste("ridge boundary", lambda)
    )
    expect_equal(
      rtemis:::linad_enet_gram(gram[["G"]], gram[["Xty"]], 0, lambda * sw),
      unname(drop(as.matrix(stats::coef(
        glmnet::glmnet(
          centered,
          residual,
          weights = w,
          alpha = 1,
          lambda = lambda,
          standardize = FALSE,
          intercept = FALSE
        ),
        s = lambda
      ))))[-1L],
      tolerance = 1e-5,
      info = paste("lasso", lambda)
    )
  }
})


test_that("linad_gram_loss() equals the residual sum of squares it stands for", {
  # The identity syy - 2 b'Xty + b'Gb is what makes a candidate cost no data
  # access. It is used everywhere and asserted nowhere.
  set.seed(8)
  for (rep in seq_len(10L)) {
    n <- sample(50:150, 1L)
    p <- sample(2:5, 1L)
    xm <- cbind(1, matrix(rnorm(n * p), n, p))
    y <- rnorm(n)
    w <- runif(n, 0.1, 3)
    gram <- rtemis:::linad_gram(xm, y, w, NULL)
    b <- rnorm(p + 1L)
    expect_equal(
      rtemis:::linad_gram_loss(gram[["G"]], gram[["Xty"]], gram[["syy"]], b),
      sum(w * (y - drop(xm %*% b))^2),
      tolerance = 1e-8
    )
  }
})


test_that("linad_ridge_edf() equals the trace it stands for", {
  set.seed(5)
  for (d in c(4L, 12L, 30L)) {
    G <- crossprod(matrix(rnorm(5 * d * d), 5 * d, d))
    for (penalty in c(0, 0.5, 10, 1000)) {
      for (intercept in c(TRUE, FALSE)) {
        ridge <- rtemis:::linad_ridge_diagonal(G, penalty, intercept)
        # tr((G + D)^-1 G), computed directly.
        expect_equal(
          rtemis:::linad_ridge_edf(G, penalty, intercept),
          sum(diag(solve(G + diag(ridge, d)) %*% G)),
          tolerance = 1e-6,
          info = paste(d, penalty, intercept)
        )
      }
    }
  }
})


# %% The floor, as a tree-level invariant ----
test_that("linad_check_tree() catches a tree holding nodes below the floor", {
  # The failure this guards is a hyperparameter that moves the split floor
  # without naming it in its description: `node_test` did, dropping the floor
  # from `min_cases_node_model` to `min_cases_leaf`, which defaults to 1.
  set.seed(2026)
  n <- 400L
  X <- as.data.frame(matrix(rnorm(n * 4L), n, 4L))
  names(X) <- paste0("x", 1:4)
  X[["y"]] <- ifelse(X[["x1"]] > 0, 3 * X[["x2"]] - 2 * X[["x3"]], 5) +
    rnorm(n, 0, 1)
  fit <- function(...) {
    train(
      X,
      hyperparameters = setup_LINAD(
        max_leaves = 8L,
        node_model = "ridge",
        min_cases_node_model = 25L,
        force_max_leaves = TRUE,
        ...
      ),
      execution_config = setup_SerialExecution(seed = 1L),
      verbosity = 0L
    )
  }
  for (rule in c("none", "aic", "bic")) {
    expect_length(
      linad_check_tree(fit(node_test = rule)@model, min_cases_child = 25L),
      0L
    )
  }
  # And it reports rather than passing silently when the floor is not met.
  expect_length(
    linad_check_tree(fit()@model, min_cases_child = 10000L),
    1L
  )
})


# %% The selected tree is the model ----
test_that("Split-gain importance counts only splits the selected tree reaches", {
  # The frame keeps every node grown while validation selects a smaller size,
  # so anything describing the model must stop at the selected terminals.
  set.seed(6)
  n <- 400L
  dat <- data.frame(a = rnorm(n), b = rnorm(n), c = rnorm(n))
  dat[["y"]] <- ifelse(dat[["b"]] < 0, -2, 2) +
    1.5 * dat[["a"]] +
    rnorm(n, 0, 0.7)
  mod <- train(
    dat,
    hyperparameters = setup_LINAD(max_leaves = 7L, force_max_leaves = TRUE),
    execution_config = setup_SerialExecution(seed = 1L),
    verbosity = 0L
  )
  total_gain <- function(k) {
    shrunk <- mod
    shrunk@model@n_leaves <- k
    importance <- varimp_super(shrunk@model)
    sum(varimp_table(importance)[["split_gain"]])
  }
  # One leaf reaches no split at all.
  expect_equal(total_gain(1L), 0)
  # And a larger selected tree never reaches less gain than a smaller one.
  gains <- vapply(seq_len(7L), total_gain, numeric(1L))
  expect_false(is.unsorted(gains))
})


test_that("A printed tree describes its selected size, not everything grown", {
  set.seed(6)
  n <- 300L
  dat <- data.frame(a = rnorm(n), b = rnorm(n))
  dat[["y"]] <- ifelse(dat[["b"]] < 0, -2, 2) + dat[["a"]] + rnorm(n, 0, 0.7)
  mod <- train(
    dat,
    hyperparameters = setup_LINAD(max_leaves = 6L, force_max_leaves = TRUE),
    execution_config = setup_SerialExecution(seed = 1L),
    verbosity = 0L
  )
  shrunk <- mod
  shrunk@model@n_leaves <- 1L
  # A single-leaf tree is one node of depth 0, whatever the frame still holds.
  expect_match(
    paste(utils::capture.output(print(shrunk@model)), collapse = " "),
    "1 node"
  )
  expect_gt(nrow(mod@model@frame), 1L)
})


# %% Classification initialization ----
test_that("The classification root is the loss-minimizing constant at every rate", {
  set.seed(2)
  n <- 400L
  dat <- data.frame(x = rnorm(n))
  dat[["y"]] <- factor(
    ifelse(runif(n) < 0.8, "pos", "neg"),
    levels = c("neg", "pos")
  )
  prevalence <- mean(dat[["y"]] == "pos")
  baseline <- 0.5 * log(prevalence / (1 - prevalence))
  for (rate in c(0, 1e-12, 0.5, 1)) {
    mod <- train(
      dat,
      hyperparameters = setup_LINAD(
        max_leaves = 1L,
        node_model = "constant",
        root_learning_rate = rate
      ),
      execution_config = setup_SerialExecution(seed = 1L),
      verbosity = 0L
    )
    expect_equal(
      mod@model@frame[["node_value"]][[1L]],
      baseline,
      tolerance = 1e-8,
      info = rate
    )
    # An intercept-only classification predicts the prevalence.
    expect_equal(
      predict(mod, dat["x"])[[1L]],
      prevalence,
      tolerance = 1e-8,
      info = rate
    )
  }
})


# %% Forest boundaries ----
test_that("A forest predicts and reports a standard error for one case", {
  set.seed(1)
  n <- 150L
  dat <- data.frame(a = rnorm(n), b = rnorm(n))
  dat[["y"]] <- 2 * dat[["a"]] + rnorm(n)
  classification <- dat
  classification[["y"]] <- factor(
    ifelse(dat[["y"]] > 0, "hi", "lo"),
    levels = c("lo", "hi")
  )
  one_row <- dat[1L, c("a", "b")]
  for (n_trees in c(1L, 3L)) {
    for (outcome in list(dat, classification)) {
      mod <- train(
        outcome,
        hyperparameters = setup_LINADForest(n_trees = n_trees, max_leaves = 4L),
        execution_config = setup_SerialExecution(seed = 1L),
        verbosity = 0L
      )
      expect_length(predict(mod, one_row), 1L)
      expect_equal(
        as.vector(predict(mod, one_row)),
        as.vector(predict(mod, dat[1:3, c("a", "b")]))[[1L]],
        tolerance = 1e-10
      )
      standard_error <- se(mod, one_row)
      expect_length(standard_error, 1L)
      if (n_trees == 1L) {
        # The infinitesimal jackknife is a covariance across bags and has
        # nothing to vary over with one tree.
        expect_true(is.na(standard_error))
      } else {
        expect_true(is.finite(standard_error) && standard_error >= 0)
      }
    }
  }
})


# %% The design basis is part of the model ----
test_that("Predictions do not depend on the session's contrast options", {
  # `model.matrix()` reads `getOption("contrasts")`, which is process-global and
  # mutable, so a design rebuilt at predict time could be a different basis
  # multiplied against the same coefficients by position.
  set.seed(4)
  n <- 300L
  dat <- data.frame(
    g = factor(sample(c("a", "b", "c"), n, replace = TRUE)),
    x = rnorm(n)
  )
  dat[["y"]] <- as.numeric(dat[["g"]]) * 2 + dat[["x"]] + rnorm(n, 0, 0.3)
  mod <- train(
    dat,
    hyperparameters = setup_LINAD(max_leaves = 4L),
    execution_config = setup_SerialExecution(seed = 1L),
    verbosity = 0L
  )
  features <- dat[c("g", "x")]
  before <- predict(mod, features)
  for (setting in list(
    c("contr.sum", "contr.poly"),
    c("contr.helmert", "contr.poly")
  )) {
    restore <- options(contrasts = setting)
    expect_equal(predict(mod, features), before, tolerance = 1e-12)
    options(restore)
  }
  # And an ordered factor is coded the same way as an unordered one.
  ordered_dat <- dat
  ordered_dat[["g"]] <- factor(
    as.character(dat[["g"]]),
    levels = c("a", "b", "c"),
    ordered = TRUE
  )
  expect_equal(
    predict(mod, ordered_dat[c("g", "x")]),
    before,
    tolerance = 1e-12
  )
})


test_that("A design that does not match the fitted one is refused", {
  set.seed(4)
  n <- 120L
  dat <- data.frame(
    g = factor(sample(c("a", "b"), n, replace = TRUE)),
    x = rnorm(n)
  )
  dat[["y"]] <- as.numeric(dat[["g"]]) + dat[["x"]] + rnorm(n, 0, 0.3)
  mod <- train(
    dat,
    hyperparameters = setup_LINAD(max_leaves = 2L),
    execution_config = setup_SerialExecution(seed = 1L),
    verbosity = 0L
  )
  expect_length(mod@model@design_names, ncol(mod@model@coefficients))
  tampered <- mod
  tampered@model@design_names <- rev(mod@model@design_names)
  expect_error(
    predict(tampered, dat[c("g", "x")]),
    "does not match the one this model was fitted on"
  )
})


# %% Cut budgets ----
test_that("n_cuts is the number of candidates a feature spends", {
  set.seed(9)
  sorted <- sort(rnorm(200L))
  breaks <- seq_len(199L)
  for (budget in c(3L, 5L, 20L)) {
    expect_length(
      linad_cut_positions(breaks, budget, sorted, "frequency"),
      budget
    )
  }
  # Fewer admissible positions than the budget returns all of them.
  expect_equal(linad_cut_positions(1:3, 20L, sorted[1:4], "frequency"), 1:3)
})


test_that("Equal-width cuts span the feature's range, tail included", {
  # Positions 1, 2 and 3 have thresholds 0.5, 1.5 and 51. Two equal-width
  # targets over the range 0 to 100 fall at 33.3 and 66.7, both nearest the
  # threshold in the tail, which collapses to one candidate.
  expect_equal(linad_cut_positions(1:3, 2L, c(0, 1, 2, 100), "width"), 3L)
  # Frequency spacing ignores the values and keeps two.
  expect_length(linad_cut_positions(1:3, 2L, c(0, 1, 2, 100), "frequency"), 2L)
  # On a uniform feature the two agree closely, which is the case that hides
  # the difference.
  uniform <- seq(0, 1, length.out = 101L)
  expect_equal(
    length(linad_cut_positions(seq_len(100L), 4L, uniform, "width")),
    length(linad_cut_positions(seq_len(100L), 4L, uniform, "frequency"))
  )
})