tests/testthat/test-loglik_math_grad_parity.R

In xnicher: Estimates Ecological Niche Models Using Ellipses

test_that("hybrid analytic gradient matches central FD", {
  occ <- as.matrix(example_env_occ_2d)
  M   <- as.matrix(example_env_m_2d)

  set.seed(3)
  n_m     <- nrow(M)
  den_idx <- sample.int(n_m, min(500L, n_m))
  kde_idx <- sample.int(n_m, min(500L, n_m))
  M_kde   <- M[kde_idx, , drop = FALSE]
  w_occ   <- as.numeric(kde_gaussian(occ, M_kde))
  w_den   <- as.numeric(kde_gaussian(M[den_idx, , drop = FALSE], M_kde))
  M_den   <- M[den_idx, , drop = FALSE]

  fd_grad <- function(theta) {
    n <- length(theta)
    g <- numeric(n)
    for (i in seq_len(n)) {
      dx <- abs(theta[i]) * 1e-6 + 1e-8
      th_p <- theta
      th_p[i] <- theta[i] + dx
      th_m <- theta
      th_m[i] <- theta[i] - dx
      fp <- xnicher:::loglik_niche_math_ip_weighted_cpp(
        th_p, occ, M_den, w_occ, w_den, eta = 1)
      fm <- xnicher:::loglik_niche_math_ip_weighted_cpp(
        th_m, occ, M_den, w_occ, w_den, eta = 1)
      g[i] <- (fp - fm) / (2 * dx)
    }
    g
  }

  for (k in seq_len(5L)) {
    theta <- start_theta(occ) + stats::rnorm(length(start_theta(occ)),
                                             sd = 0.05)
    res <- xnicher:::loglik_niche_math_ip_weighted_grad_cpp(
      theta, occ, M_den, w_occ, w_den, eta = 1)
    g_fd <- fd_grad(theta)

    # Tolerance dominated by FD precision (~1e-6 relative)
    expect_equal(res$gradient, g_fd, tolerance = 1e-5)
  }
})

test_that("penalized weighted hybrid analytic gradient matches central FD", {
  occ <- as.matrix(example_env_occ_2d)
  M   <- as.matrix(example_env_m_2d)

  set.seed(31)
  n_m     <- nrow(M)
  den_idx <- sample.int(n_m, min(400L, n_m))
  kde_idx <- sample.int(n_m, min(400L, n_m))
  M_kde   <- M[kde_idx, , drop = FALSE]
  w_occ   <- as.numeric(kde_gaussian(occ, M_kde))
  w_den   <- as.numeric(kde_gaussian(M[den_idx, , drop = FALSE], M_kde))
  M_den   <- M[den_idx, , drop = FALSE]

  log_sd <- log(apply(occ, 2L, stats::sd))
  lam    <- 1.5

  fd_grad_pen <- function(theta) {
    n <- length(theta)
    g <- numeric(n)
    for (i in seq_len(n)) {
      dx <- abs(theta[i]) * 1e-6 + 1e-8
      th_p <- theta
      th_p[i] <- theta[i] + dx
      th_m <- theta
      th_m[i] <- theta[i] - dx
      fp <- xnicher:::loglik_niche_math_weighted_cpp(
        th_p, occ, M_den, w_occ, w_den,
        prior_log_sigma_center = log_sd, prior_log_sigma_lambda = lam,
        eta = 1)
      fm <- xnicher:::loglik_niche_math_weighted_cpp(
        th_m, occ, M_den, w_occ, w_den,
        prior_log_sigma_center = log_sd, prior_log_sigma_lambda = lam,
        eta = 1)
      g[i] <- (fp - fm) / (2 * dx)
    }
    g
  }

  for (k in seq_len(5L)) {
    theta <- start_theta(occ) + stats::rnorm(length(start_theta(occ)),
                                             sd = 0.05)
    res <- xnicher:::loglik_niche_math_weighted_grad_cpp(
      theta, occ, M_den, w_occ, w_den,
      prior_log_sigma_center = log_sd, prior_log_sigma_lambda = lam,
      eta = 1)
    g_fd <- fd_grad_pen(theta)

    expect_equal(res$gradient, g_fd, tolerance = 1e-5)
  }
})