tests/testthat/test-laplace-gp-dispatch.R

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

# test-laplace-gp-dispatch.R
# Pin the dispatch contract for tulpa_laplace(spatial = spatial_gp(...)).
# Builds a tiny pre-validated GP spec by hand so the test does not depend
# on the full validate_gp(data) chain or any mesh package.

# Helper: build a synthetic NNGP neighbor structure for a small problem.
make_synthetic_gp_spec <- function(n_obs = 50, k = 5L, seed = 11L) {
  set.seed(seed)
  coords <- cbind(runif(n_obs), runif(n_obs))
  ord <- order(coords[, 1], coords[, 2])
  coords_ord <- coords[ord, , drop = FALSE]

  nn_idx  <- matrix(0L, n_obs, k)
  nn_dist <- matrix(0,  n_obs, k)
  for (i in 2:n_obs) {
    d <- sqrt((coords_ord[1:(i - 1), 1] - coords_ord[i, 1])^2 +
              (coords_ord[1:(i - 1), 2] - coords_ord[i, 2])^2)
    nc <- min(length(d), k)
    sel <- order(d)[seq_len(nc)]
    nn_idx[i, seq_len(nc)]  <- sel
    nn_dist[i, seq_len(nc)] <- d[sel]
  }

  spec <- spatial_gp(coords = ~ x1 + x2, cov = "exponential", nn = k)
  # Manually populate the fields validate_gp() would fill.
  spec$n_obs        <- n_obs
  spec$n_spatial    <- n_obs
  spec$n_unique     <- n_obs
  spec$obs_to_loc   <- seq_len(n_obs)
  spec$unique_coords <- coords_ord
  spec$coords_matrix <- coords_ord
  spec$nn           <- k
  spec$neighbor_info <- list(
    nn_idx       = nn_idx,
    nn_dist      = nn_dist,
    nn_order     = ord,
    nn_order_inv = order(ord)
  )
  list(spec = spec, n_obs = n_obs, coords = coords_ord, ord = ord)
}

test_that("laplace_gp_at refuses an unvalidated spec with a clear error", {
  raw_spec <- spatial_gp(coords = ~ x1 + x2)
  expect_error(
    laplace_gp_at(
      y = rbinom(10, 1, 0.5), n_trials = rep(1L, 10),
      X = matrix(1, 10, 1), spatial = raw_spec
    ),
    "neighbor_info is NULL"
  )
})

skip_on_cran_budget()

test_that("dispatch_laplace_spatial routes GP specs (no longer errors)", {
  skip_on_cran_budget()
  ss <- make_synthetic_gp_spec()
  X  <- matrix(1.0, ss$n_obs, 1)
  y  <- rbinom(ss$n_obs, 1, 0.5)

  res <- tryCatch(
    dispatch_laplace_spatial(
      y = y, n_trials = rep(1L, ss$n_obs), X = X,
      re_idx = NULL, n_re_groups = 0L, sigma_re = 1.0,
      spatial = ss$spec, family = "binomial", phi = 1.0,
      max_iter = 5L, tol = 1e-3, n_threads = 1L
    ),
    error = function(e) e
  )
  msg <- if (inherits(res, "error")) conditionMessage(res) else NA_character_

  # Asserted on both branches: the dispatcher may reach the kernel and return,
  # or the kernel may raise its own error, but the legacy dispatch gap
  # ("not yet supported") must never come back either way.
  expect_false(isTRUE(grepl("not yet supported", msg, fixed = TRUE)))
  if (is.na(msg)) expect_true(is.list(res))
})

test_that("gp_cov_type maps cov + nu correctly", {
  skip_on_cran_budget()
  expect_equal(gp_cov_type(list(cov = "exponential")), 0L)
  expect_equal(gp_cov_type(list(cov = "matern", nu = 1.5)), 1L)
  # 4, not 2: 2 is the Gaussian kernel in tulpa::CovType, and handing a
  # Matern-5/2 request that code is what fitted it as a squared exponential on
  # every sampler mode (gcol33/tulpa#686).
  expect_equal(gp_cov_type(list(cov = "matern", nu = 2.5)), 4L)

  # Reject unsupported covariances with an actionable error.
  expect_error(gp_cov_type(list(cov = "matern", nu = 0.5)),
               "Matern with nu in")
  # The one code names the same kernel on the Laplace scatter, the PG sweep and
  # the sampler, so a fit and a prediction from it are the same model.
  expect_equal(cpp_test_compute_cov(0.7, 2.0, 1.3, 4L),
               2.0 * (1 + sqrt(5) * 0.7 / 1.3 + (sqrt(5) * 0.7 / 1.3)^2 / 3) *
                 exp(-sqrt(5) * 0.7 / 1.3))
  expect_error(gp_cov_type(list(cov = "gaussian")),
               "exponential.*matern")
  expect_error(gp_cov_type(list(cov = "spherical")),
               "exponential.*matern")
})

skip_on_cran()

test_that("NNGP Laplace threads an additional iid RE block (issue #74)", {
  skip_on_cran()
  ss <- make_synthetic_gp_spec()
  X  <- matrix(1.0, ss$n_obs, 1)
  set.seed(7)
  y  <- rbinom(ss$n_obs, 1, 0.5)
  n_re <- 3L
  re_idx <- sample.int(n_re, ss$n_obs, replace = TRUE)

  fit <- dispatch_laplace_spatial(
    y = y, n_trials = rep(1L, ss$n_obs), X = X,
    re_idx = re_idx, n_re_groups = n_re, sigma_re = 1.0,
    spatial = ss$spec, family = "binomial", phi = 1.0,
    max_iter = 50L, tol = 1e-6, n_threads = 1L
  )

  # Latent layout [beta (p), re (n_re), w_gp (n_spatial)].
  expect_length(fit$mode, ncol(X) + n_re + ss$spec$n_spatial)
  expect_true(is.finite(fit$log_marginal))
})

test_that("tulpa_laplace(spatial = spatial_gp(...)) runs end-to-end", {
  skip_on_cran()
  ss <- make_synthetic_gp_spec(n_obs = 80)
  set.seed(101)
  X  <- cbind(1, rnorm(ss$n_obs))
  beta_true <- c(-0.3, 0.5)
  eta <- as.numeric(X %*% beta_true)
  y  <- rbinom(ss$n_obs, 1, plogis(eta))

  fit <- tulpa_laplace(
    y = y, n_trials = rep(1L, ss$n_obs), X = X,
    re_list = list(),
    family = "binomial",
    spatial = ss$spec,
    max_iter = 50L, tol = 1e-6, n_threads = 1L,
    return_hessian = TRUE
  )

  # mode = c(beta, spatial_effects) — spatial_effects has length n_spatial.
  expect_length(fit$mode, ncol(X) + ss$spec$n_spatial)
  expect_true(is.finite(fit$log_marginal))
  # Since issue #16: H_beta is the Schur'd marginal precision for NNGP via
  # Q_nngp = (I-A)' D^{-1} (I-A). SE quality is exercised in
  # dev_notes/test_issue16_nngp_se.R against tulpaObs end-to-end; here the
  # property is that it is produced at all. This was guarded behind
  # `if (!is.null(fit$H_beta))`, which asserts nothing in exactly the case
  # worth catching: the builder degrades to a warning and a NULL, so the
  # conditional turned a silent loss of every fixed-effect SE into a pass.
  expect_false(is.null(fit$H_beta))
  expect_equal(dim(fit$H_beta), c(ncol(X), ncol(X)))
  expect_true(all(is.finite(fit$H_beta)))
})

test_that("tulpa_laplace(spatial = spatial_gp(...)) converges on the sparse path", {
  skip_on_cran()
  # n_x = p + n_spatial must clear SPARSE_THRESHOLD (200) so we route through
  # the CHOLMOD sparse Newton path. Earlier the NNGP precision scatter
  # double-counted off-diagonals between neighbours (storage is lower-triangle
  # only via SparseHessianBuilder, so (k, kp) and (kp, k) hit the same slot
  # and both got added), which made every iteration's H indefinite. CHOLMOD
  # then warned 'matrix not positive definite' on every Newton step and the
  # mode stayed at zero. This test pins the fixed behaviour: convergence on a
  # sparse-path problem with non-trivial β and a spatial field that moves.
  ss <- make_synthetic_gp_spec(n_obs = 250L, k = 5L, seed = 11L)
  set.seed(101)
  X  <- cbind(1, rnorm(ss$n_obs))
  beta_true <- c(-0.3, 0.5)
  eta <- as.numeric(X %*% beta_true)
  y  <- rbinom(ss$n_obs, 1, plogis(eta))

  fit <- expect_warning(
    tulpa_laplace(
      y = y, n_trials = rep(1L, ss$n_obs), X = X,
      re_list = list(),
      family = "binomial",
      spatial = ss$spec,
      max_iter = 50L, tol = 1e-6, n_threads = 1L
    ),
    regexp = NA
  )

  expect_true(isTRUE(fit$converged))
  expect_length(fit$mode, ncol(X) + ss$spec$n_spatial)
  expect_true(is.finite(fit$log_marginal))
  expect_gt(max(abs(fit$mode[-seq_len(ncol(X))])), 1e-6)
  expect_lt(max(abs(fit$mode[seq_len(ncol(X))] - beta_true)), 0.6)
})

test_that("GP Laplace recovers the field with repeated coordinates (obs_to_loc)", {
  skip_on_cran()
  # Multiple observations share each unique location, so n_spatial < N and the
  # single-point GP kernel must attach each obs to its field node via
  # obs_to_loc. With the old identity map (obs i -> node i) later observations
  # were dropped and the earlier ones mis-attached, so the field did not track
  # truth.
  set.seed(1)
  n_loc <- 25L; reps <- 4L
  coords_u <- cbind(runif(n_loc), runif(n_loc))
  D <- as.matrix(dist(coords_u)); K <- exp(-D / 0.35)
  f_u <- as.numeric(t(chol(K + 1e-6 * diag(n_loc))) %*% rnorm(n_loc))
  f_u <- f_u - mean(f_u)
  loc <- rep(seq_len(n_loc), each = reps)
  coords <- coords_u[loc, , drop = FALSE]
  y <- rpois(length(loc), exp(0.5 + f_u[loc]))
  df <- data.frame(y = y, lon = coords[, 1], lat = coords[, 2])

  fit <- tulpa(y ~ 1, data = df, family = "poisson",
               spatial = spatial_gp(~ lon + lat), mode = "laplace")

  # Field at the unique locations = trailing n_loc latent entries.
  fhat <- tail(fit$mode, n_loc)
  expect_length(fhat, n_loc)
  expect_gt(stats::cor(fhat, f_u), 0.6)          # field tracks truth
  expect_lt(abs(fit$mode[1] - 0.5), 0.3)          # intercept recovered
})