tests/testthat/test_harmonics.R

In sphunif: Uniformity Tests on the Circle, Sphere, and Hypersphere

set.seed(12314)
int_prod <- function(x, i1, i2, k1, k2) {
  g_i_k(x = x, k = k1, i = i1) * g_i_k(x = x, k = k2, i = i2)
}

test_that("d_p_k for p = 2 and k = 0,1", {

  expect_equal(d_p_k(p = 2, k = 0:3), c(1, rep(2, 3)))
  expect_equal(d_p_k(p = 2, k = 0:3, log = TRUE), c(0, rep(log(2), 3)))
  expect_equal(d_p_k(p = 30, k = 0:1), c(1, 30))

})

test_that("d_p_k = (1 + (2 * k) / (p - 2)) * C_k^{(p - 2) / 2}(1)", {

  k <- 0:5
  for (p in 2:5) {
    expect_equal(d_p_k(p = p, k = k),
                 drop(Gegen_polyn(theta = 0, k = k, p = p)) *
                   switch((p == 2) + 1, 1 + (2 * k) / (p - 2), 2 - (k == 0)))
  }

})

test_that("omega{p - 2} / C_k^{(p - 2) / 2}(1) =
           1 / c_p_k * (1 + (2 * k) / (p - 2))^(-1) * omega_{p - 1}", {

  k <- 0:5
  for (p in 2:5) {
   expect_equal(rotasym::w_p(p = p - 1) /
                  drop(Gegen_polyn(theta = 0, k = k, p = p)),
                1 / drop(Gegen_coefs(only_const = TRUE, p = p, k = k)) *
                switch((p == 2) + 1, (1 + (2 * k) / (p - 2)),
                       2 - (k == 0))^(-1) * rotasym::w_p(p = p))
  }

})

test_that("c_p_k = (1 + (2 * k) / (p - 2))^(-2) * d_p_k *
           omega_{p - 1} / omega{p - 2}", {

  k <- 0:5
  for (p in 2:5) {
    expect_equal(drop(Gegen_coefs(only_const = TRUE, p = p, k = k)),
                 switch((p == 2) + 1, (1 + (2 * k) / (p - 2)),
                        2 - (k == 0))^(-2) * d_p_k(p = p, k = k) *
                   rotasym::w_p(p = p) / rotasym::w_p(p = p - 1))
  }

})

test_that("angles_to_sphere vs. sphere_to_angles", {

  n <- 10
  for (p in c(2:4, 11)) {
    X <- r_unif_sph(n = n, p = p, M = 1)[, , 1]
    th <- t(matrix(runif(n * (p - 1)), nrow = p - 1, ncol = n) *
              c(rep(pi, p - 2), 2 * pi))
    expect_equal(angles_to_sphere(sphere_to_angles(x = X)), X)
    expect_equal(sphere_to_angles(angles_to_sphere(theta = th)), th)
    expect_equal(angles_to_sphere(sphere_to_angles(x = X[1, ])),
                 X[1, , drop = FALSE])
    expect_equal(sphere_to_angles(angles_to_sphere(theta = th[1, ])),
                 th[1, , drop = FALSE])

  }

})

test_that("Harmonics with (i, k) vs. m", {

  skip_on_cran()
  for (p in 2:5) {
    repeat {
      m <- c(rpois(n = p - 1, lambda = 2), sample(c(0, 1), size = 1))
      if (sum(m) > 0) break
    }
    k <- sum(m)
    ms <- as.matrix(do.call(expand.grid, c(rep(list(0:k), p - 1), list(0:1))))
    ms <- ms[rowSums(ms) == k, ]
    i <- which(apply(ms, 1, function(x) all(x == m)))
    X <- r_unif_sph(n = 1, p = p, M = 1)[, , 1]
    expect_message(expect_equal(
      g_i_k(x = X, i = i, k = k, show_m = TRUE), g_i_k(x = X, m = m)))
  }

})

test_that("Orthonormality of harmonics among i's with common k", {

  skip_on_cran()
  set.seed(323231)
  M <- 1e3
  for (p in 2:5) {
    for (k in 1:2) {

      dpk <- d_p_k(p = p, k = k)
      intprods <- diag(rep(1, dpk))
      for (i1 in 1:dpk) {
        for (i2 in i1:dpk) {
          intprods[i1, i2] <- int_sph_MC(f = int_prod, i1 = i1, i2 = i2,
                                         k1 = k, k2 = k, p = p, M = M,
                                         cores = 1) /
            rotasym::w_p(p = p)
          image(abs(intprods), zlim = c(0, 1.5))
          expect_lt(max(abs(intprods - diag(rep(1, dpk)))), 0.1)
        }
      }

    }
  }

})

test_that("Orthonormality of harmonics among i's with different k's = 1,2", {

  skip_on_cran()
  set.seed(323231)
  M <- 1e3
  for (p in 2:5) {

    dpk1 <- d_p_k(p = p, k = 1)
    dpk2 <- d_p_k(p = p, k = 2)
    intprods <- matrix(0, nrow = dpk1, ncol = dpk2)
    for (i1 in 1:dpk1) {
      for (i2 in 1:dpk2) {
        intprods[i1, i2] <- int_sph_MC(f = int_prod, i1 = i1, i2 = i2,
                                       k1 = 1, k2 = 2, p = p, M = M,
                                       cores = 1) /
          rotasym::w_p(p = p)
        image(abs(intprods), zlim = c(0, 1.5))
        expect_lt(max(abs(intprods)), 0.1)
      }
    }

  }

})

test_that("Orthonormality of harmonics among i's with different k's = 0,1", {

  skip_on_cran()
  set.seed(323231)
  M <- 1e3
  for (p in 2:5) {

    dpk0 <- d_p_k(p = p, k = 0)
    dpk1 <- d_p_k(p = p, k = 1)
    intprods <- matrix(0, nrow = dpk0, ncol = dpk1)
    for (i0 in 1:dpk0) {
      for (i1 in 1:dpk1) {
        intprods[i0, i1] <- int_sph_MC(f = int_prod, i1 = i0, i2 = i1,
                                       k1 = 0, k2 = 1, p = p, M = M,
                                       cores = 1) /
          rotasym::w_p(p = p)
        image(abs(intprods), zlim = c(0, 1.5))
        expect_lt(max(abs(intprods)), 0.1)
      }
    }

  }

})

test_that("Specific cases of harmonics for x = e_1", {

  for (p in 2:5) {
    for (k in 1:3) {
      for (i in 1:d_p_k(p = p, k = k)) {
        expect_equal(g_i_k(x = c(1, rep(0, p - 1)), i = i, k = k),
                     (i == 1) * sqrt((gamma(k + p - 2) * (2 * k + p - 2)) /
                                       (gamma(k + 1) * gamma(p - 1))))
      }
    }
  }

})

test_that("Funk-Hecke formula in Dai and Xu (2013)", {

  skip_on_cran()
  set.seed(323231)
  M <- 1e3
  kappa <- 2
  i0 <- 1
  for (p in 2:5) {
    for (k0 in 1:3) {
      mu <- c(1, rep(0, p - 1))
      psi <- function(x) rotasym::g_vMF(t = x, p = p, kappa = kappa)
      expect_equal(
        mean(g_i_k(x = rotasym::r_vMF(n = M, mu = mu, kappa = kappa),
                     i = i0, k = k0)) / rotasym::w_p(p = p),
        rotasym::w_p(p = p - 1) / rotasym::w_p(p = p) *
          integrate(function(t) psi(x = t) *
                      Gegen_polyn(theta = acos(t), k = k0, p = p) *
                      (1 - t^2)^((p - 3) / 2), lower = -1, upper = 1,
                    stop.on.error = FALSE)$value /
          drop(Gegen_polyn(theta = 0, k = k0, p = p)) *
          drop(g_i_k(x = mu, i = i0, k = k0)),
        tolerance = 0.05
      )
    }
  }

})

test_that("Alternative version Funk-Hecke formula", {

  skip_on_cran()
  set.seed(323231)
  M <- 1e3
  kappa <- 2
  i0 <- 1
  for (p in 2:5) {
    for (k0 in 1:3) {
    mu <- c(1, rep(0, p - 1))
    phi <- function(x) rotasym::g_vMF(t = x, p = p, kappa = kappa)
    expect_equal(
      mean(g_i_k(x = rotasym::r_vMF(n = M, mu = mu, kappa = kappa),
                 i = i0, k = k0)) / rotasym::w_p(p = p),
      switch((p == 2) + 1, (1 + (2 * k0) / (p - 2)),
             2 - (k0 == 0))^(-1) *
        drop(Gegen_coefs(k = k0, p = p, psi = function(th) phi(x = cos(th)))) *
        drop(g_i_k(x = mu, i = i0, k = k0)),
      tolerance = 0.05)
    }
  }

})