tests/testthat/test_ConversionOps.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

# Regression tests for three fixes (2026-07-02):
# (1) as.BDataLogPDFs.BDataLogDensities used Sub instead of Add (sign of log
#     reference density flipped; invisible on B2Unif01 where the reference pdf
#     is 1, wrong on any non-uniform reference).
# (2) Add.BDataLogPDFs wrapped its (correct) log-scale result in BDataDensities
#     instead of BDataLogPDFs.
# (3) dist.BDataSplines computed Euclidean distances on raw ZB-spline
#     coefficients, ignoring the Gram matrix S of the (non-orthonormal) basis
#     (finding B7); it now maps coefficients to orthonormal coordinates first.

test_that("as.BDataLogPDFs adds (not subtracts) the log reference density", {
  # density identically 1 w.r.t. N(0,1) reference => the distribution is N(0,1)
  x <- BDataLogDensities(list(function(z) rep(0, length(z))), basis = B2N01Basis5)
  lp <- as.BDataLogPDFs(x)
  z <- c(-1.5, -0.5, 0, 0.5, 1.5)
  expect_equal(lp[[1]](z), dnorm(z, log = TRUE), tolerance = 1e-10)
  # consistency with the exp route
  p <- as.BDataPDFs(x)
  expect_equal(lp[[1]](z), log(p[[1]](z)), tolerance = 1e-10)
})

test_that("Add on BDataLogPDFs returns a BDataLogPDFs with correct values", {
  lp <- BDataLogPDFs(list(function(z) dnorm(z, log = TRUE)), basis = B2N01Basis5)
  s <- lp + lp
  expect_s3_class(s, "BDataLogPDFs")
  # perturbation of N(0,1) by itself: log pdf doubles up to the reference part;
  # check against the density route which is known to be correct
  d <- as.BDataDensities(lp) + as.BDataDensities(lp)
  z <- c(-1, 0, 1)
  expect_equal(exp(s[[1]](z)), as.BDataPDFs(d)[[1]](z), tolerance = 1e-10)
})

test_that("dist on BDataSplines equals the Gram-matrix (Bayes-space) distance", {
  basis <- B2Unif01Basis7
  S <- getMetricTensor(basis)
  # non-orthonormal ZB-spline basis: S must differ from the identity
  expect_false(isTRUE(all.equal(S, diag(nrow(S)))))
  set.seed(1)
  coefs <- matrix(rnorm(3 * nrow(S)), nrow = 3)
  x <- BDataSplines(coefs, basis = basis)
  D <- as.matrix(dist(x))
  for (i in 1:2) for (j in (i + 1):3) {
    delta <- coefs[i, ] - coefs[j, ]
    expect_equal(D[i, j], sqrt(drop(t(delta) %*% S %*% delta)), tolerance = 1e-10)
  }
  # consistency with scalarProduct (which is Gram-aware)
  G <- scalarProduct(x, x)
  d12 <- sqrt(G[1, 1] + G[2, 2] - 2 * G[1, 2])
  # scalarProduct centers via cov-like machinery only for function classes;
  # for splines it is the plain quadratic form, so this must match dist exactly
  expect_equal(D[1, 2], d12, tolerance = 1e-10)
})

test_that("Add/Sub work rowwise on multi-row BDataSplines (addDatasets branch)", {
  # addDatasets() previously fell through to NULL when both operands were
  # matrices with more than one row, breaking x + y / x - y for n > 1 and,
  # downstream, spatialDepth(), median() and boxplot() on spline data
  basis <- B2Unif01Basis5
  set.seed(2)
  a <- BDataSplines(matrix(rnorm(10), nrow = 2), basis = basis)
  b <- BDataSplines(matrix(rnorm(10), nrow = 2), basis = basis)
  s <- a - b
  expect_s3_class(s, "BDataSplines")
  expect_equal(unclass(s), unclass(a) - unclass(b), ignore_attr = TRUE)
  expect_error(a - BDataSplines(matrix(rnorm(15), nrow = 3), basis = basis),
               "different number of rows")
  # depth toolbox runs end-to-end on multi-row spline data
  d <- spatialDepth(BDataSplines(matrix(rnorm(20), nrow = 4), basis = basis))
  expect_length(d, 4)
  expect_true(all(d >= 0 & d <= 1))
})

test_that("Mul/Div on BDataSplines do not recurse (is.numeric(x) guard)", {
  # Mul.BDataSplines used to test is.numeric(x); since a BDataSplines IS a
  # numeric matrix, Mul/Div (and hence * and /) entered infinite mutual
  # recursion with Mul.numeric.BData and died with a C stack overflow
  set.seed(4)
  x <- BDataSplines(matrix(rnorm(10), nrow = 2), basis = B2Unif01Basis5)
  m <- x * 2
  expect_s3_class(m, "BDataSplines")
  expect_equal(unclass(m), unclass(x) * 2, ignore_attr = TRUE)
  d <- x / c(1, 2)
  expect_s3_class(d, "BDataSplines")
  expect_equal(unclass(d)[2, ], unclass(x)[2, ] / 2, ignore_attr = TRUE)
  expect_error(x * "a", "works only with numeric")
})

test_that("norm.BDataSplines matches the Monte-Carlo clr-L2 norm", {
  set.seed(9)
  xs <- BDataSplines(matrix(rnorm(10), nrow = 2), basis = B2Unif01Basis5)
  fast <- norm(xs)
  slow <- norm(as.BDataLogDensities(xs))  # dispatches to the MC integral path
  expect_equal(fast, slow, tolerance = 0.02)
  # exact value against the quadratic form with the metric tensor
  S <- getMetricTensor(B2Unif01Basis5)
  expect_equal(fast[1], sqrt(drop(t(unclass(xs)[1, ]) %*% S %*% unclass(xs)[1, ])),
               tolerance = 1e-12)
})

test_that("background is passed through by the histogram converters", {
  # as.BDataDensities.BDataHistograms (and the two FrequencyHistograms
  # methods) hard-coded background = 0.0 instead of forwarding the argument
  h <- BDataFrequencyHistograms(matrix(c(0, 5, 5, 0), nrow = 1),
                                breaks = seq(0, 1, by = 0.25),
                                basis = B2Unif01Basis5)
  d0 <- as.BDataDensities(h)
  db <- as.BDataDensities(h, background = 1)
  expect_equal(d0[[1]](0.1), 0)          # empty bin stays zero without background
  expect_gt(db[[1]](0.1), 0)             # background lifts the empty bin
})

test_that("msd dispatches to the BData method (sqrt of metric variance)", {
  # msd.BData lacked @export/@method tags, so msd() on BData objects silently
  # fell through to compositions::msd on the raw coefficient matrix
  set.seed(5)
  x <- BDataSplines(matrix(rnorm(20), nrow = 4), basis = B2Unif01Basis5)
  expect_equal(msd(x), sqrt(mvar(x)), tolerance = 1e-12)
})

test_that("dist on an orthonormal-basis BDataSplines is plain Euclidean", {
  # identity Gram matrix => fast path, no transformation
  bs <- BayesSpaceBasis(function(x) cbind(sin(2 * pi * x), cos(2 * pi * x)),
                        bayesSpace = B2Unif01, S = diag(2))
  coefs <- matrix(c(0, 0, 3, 4), nrow = 2, byrow = TRUE)
  x <- BDataSplines(coefs, basis = bs)
  expect_equal(as.matrix(dist(x))[1, 2], 5)
})