tests/testthat/test_SplineProjection.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

library(testthat)
library(DistributionalDDA)

## Regression test for Issue 3 (TODO_orthonormality_issue.md):
## .SplineCoefsFromLogDensity must return the L2 projection coefficients
## c = solve(S, <f, zeta>), not the raw inner products <f, zeta>.
## ZB-splines are NOT orthonormal, so the Gram matrix S = getMetricTensor(basis)
## is not the identity and the solve(S, .) step is essential. Before the fix the
## function returned S c instead of c (reconstruction errors of 50-750%).
##
## Design: build a clr/log-density that lies EXACTLY in the basis span,
## f = sum_i c_true_i * zeta_i = coef2function(c_true, basis). Then by
## construction <f, zeta_j> = (S c_true)_j, so the correct projection
## solve(S, <f,zeta>) recovers c_true, whereas the raw inner products (the
## pre-fix return value) equal S c_true != c_true.
##
## The reference measure carries a SEEDED, fine Monte-Carlo sample so that
## CovarianceOperator() -> cov() is deterministic; its quadrature error scales
## like 1/sqrt(nsim) and is ~0.014 at nsim = 2e5, so tolerance 0.1 has ample
## margin while the pre-fix value (~10.9) misses it by ~100x. We must NOT use the
## packaged default basis here: its measure sample is drawn at load time without
## a seed and is therefore not reproducible across sessions.

test_that(".SplineCoefsFromLogDensity recovers in-span coefficients (Issue 3, solve(S, .))", {
  set.seed(123)
  measure <- ReferenceMeasureUniform(0, 1, sample = stats::runif(2e5, 0, 1))
  basis   <- BayesSpaceBasis(BayesSpace(referenceMeasure = measure),
                             nbasis = 7, order = 4)

  S <- getMetricTensor(basis)
  expect_gt(max(abs(S - diag(nrow(S)))), 1)          # sanity: ZB-splines, S != I

  p      <- nrow(S)
  c_true <- sin(seq_len(p))
  f      <- coef2function(c_true, basis)             # f = sum_i c_true_i * zeta_i

  ## the raw inner products (== the pre-fix return value) are S c_true, far from
  ## c_true; this invariant documents exactly what the solve(S, .) step corrects
  b <- as.vector(CovarianceOperator(measure, f, as.function(basis)))
  expect_gt(sqrt(sum((b - c_true)^2)), 1)

  ## the fixed projection recovers c_true (pre-fix err ~10.9; post-fix ~0.014)
  c_hat <- as.vector(.SplineCoefsFromLogDensity(f, basis))
  expect_lt(sqrt(sum((c_hat - c_true)^2)), 0.1)
})

test_that("density/log-density -> splines is shape-safe for multiple observations and basis changes (Issue 3)", {
  ## The projection applies solve(S, .) ROW-WISE: CovarianceOperator returns an
  ## n x p inner-product matrix (one row per observation) on the multi-observation
  ## / basis-change reproject path (as.BDataSplines onto a different basis;
  ## scalarProduct across two bases). A single-observation-only implementation
  ## (e.g. as.vector(solve(S, .))) flattens that matrix and either crashes or
  ## returns a malformed object -- these assertions guard against that regression.
  set.seed(7)
  p5 <- getBayesSpaceBasisDimension(B2Unif01Basis5)
  p7 <- getBayesSpaceBasisDimension(B2Unif01Basis7)

  ## multi-row splines on basis5, reprojected onto basis7: shape must be preserved
  X5 <- BDataSplines(matrix(rnorm(3 * p5), nrow = 3), B2Unif01Basis5)
  Y7 <- as.BDataSplines(X5, basis = B2Unif01Basis7)
  expect_equal(numberOfCases(Y7), 3)
  expect_equal(ncol(unclass(Y7)), p7)

  ## a single-row basis change must stay exactly one case (not collapse to a vector)
  X1 <- BDataSplines(matrix(rnorm(p5), nrow = 1), B2Unif01Basis5)
  expect_equal(numberOfCases(as.BDataSplines(X1, basis = B2Unif01Basis7)), 1)

  ## scalarProduct between two proper spline objects on different bases reprojects
  ## one onto the other internally (n > 1); it must run row-wise without crashing
  ## (pre-fix the as.vector() flattening crashed solve() here)
  Y7p <- BDataSplines(matrix(rnorm(3 * p7), nrow = 3), B2Unif01Basis7)
  sp  <- scalarProduct(X5, Y7p)
  expect_equal(dim(sp), c(3L, 3L))
  expect_true(all(is.finite(sp)))
})