tests/testthat/test_FunctionBasisGram.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

library(testthat)
library(DistributionalDDA)

## Regression test for B8 (TODO_orthonormality_issue.md):
## BayesSpaceBasis.function must store the GRAM matrix
## S = A %*% scalarProductOfLogDensities(basis, space) %*% t(A), NOT its inverse.
## Every consumer (scalarProduct, mvar, princomp via chol(S)) evaluates c^T S c,
## so storing solve(...) made inner products on function-constructed bases come
## out with the wrong magnitude and frequently the wrong sign. This mirrors the
## .basisfd convention already pinned in test_GramMatrix.R. The test FAILS on the
## pre-fix code (solve() present) and PASSES after dropping the solve().

## deterministic, strongly non-orthogonal 2-function basis on Unif[0,1]
## (Legendre-type ramps; the cov-based Gram is close to [[.5,.5],[.5,1]])
.rawfun <- function(x) {
  f1 <- sqrt(0.5) * sqrt(12)  * (x - 0.5)
  e  <- sqrt(0.5) * sqrt(180) * ((x - 0.5)^2 - 1 / 12)
  cbind(f1, f1 + e)
}

## trapezoidal rule on an equispaced grid (matches tests/testthat/test_GramMatrix.R)
.trapz <- function(y, dx) dx * (sum(y) - (y[1] + y[length(y)]) / 2)

test_that("BayesSpaceBasis.function dispatches to the .function method", {
  b <- BayesSpaceBasis(.rawfun)
  expect_s3_class(b, "BayesSpaceBasis")
  expect_false(b$centered)                  # .function sets FALSE; .basisfd sets TRUE
})

test_that("getMetricTensor of a function-basis is the Gram, not its inverse (B8)", {
  b <- BayesSpaceBasis(.rawfun)
  A <- b$A
  ## (i) equals the constructor's own Gram formula (same in-session sample)
  G <- A %*% scalarProductOfLogDensities(.rawfun, space = B2Unif01) %*% t(A)
  expect_equal(unname(getMetricTensor(b)), unname(G), tolerance = 1e-8)
  ## (ii) and an INDEPENDENT fine-grid quadrature of the Gram (A = I here, so
  ## zeta = .rawfun), so the check in (i) is not merely the implementation's own
  ## formula echoed back
  tt <- seq(0, 1, length.out = 20001); dx <- tt[2] - tt[1]
  Z  <- .rawfun(tt)
  Gq <- matrix(0, 2, 2)
  for (i in 1:2) for (j in 1:2) Gq[i, j] <- .trapz(Z[, i] * Z[, j], dx)
  expect_equal(unname(getMetricTensor(b)), unname(Gq), tolerance = 0.1)
  ## the real regression guard: the pre-fix inverse Gram is off by O(1), far
  ## outside any quadrature tolerance, so this expect_false flips on a revert
  expect_false(isTRUE(all.equal(unname(getMetricTensor(b)), unname(solve(G)),
                                tolerance = 1e-6)))
})

test_that("scalarProduct on a function-basis equals the Gram inner product c1^T G c2 (B8)", {
  b  <- BayesSpaceBasis(.rawfun)
  A  <- b$A
  G  <- A %*% scalarProductOfLogDensities(.rawfun, space = B2Unif01) %*% t(A)
  c1 <- c(-3, 1); c2 <- c(-1, -2)           # deterministic, mixed sign
  x  <- BDataSplines(t(c1), b)
  y  <- BDataSplines(t(c2), b)
  ip_gram    <- as.numeric(t(c1) %*% G        %*% c2)
  ip_inverse <- as.numeric(t(c1) %*% solve(G) %*% c2)
  expect_equal(as.numeric(scalarProduct(x, y)), ip_gram, tolerance = 1e-8)
  ## the pre-fix code returned ip_inverse, which has the opposite sign here
  expect_gt(ip_gram, 0)
  expect_lt(ip_inverse, 0)
})

test_that("default ZB-spline bases (.basisfd) are unaffected by the B8 fix", {
  expect_true(B2Unif01Basis7$centered)
  S7  <- getMetricTensor(B2Unif01Basis7)
  bsp <- fda::create.bspline.basis(c(0, 1), nbasis = 8, norder = 4)
  expect_equal(unname(S7),
               unname(B2Unif01Basis7$A %*% fda::inprod(bsp, bsp) %*% t(B2Unif01Basis7$A)),
               tolerance = 1e-8)
})