tests/testthat/test-admix-global.R

In AdmixPoly: Global and Local Admixture Inference in Polyploids

# =========================================================================
# Shared Mock Data Generator for Scoped Tests
# =========================================================================
create_mock_geno <- function(N=4, M=6, L=3) {
  set.seed(42)
  ind_names <- paste0("Ind", 1:N)
  allele_names <- paste0("Allele", 1:L)
  
  geno_list <- list()
  for (m in 1:M) {
    # Generate mock dosage matrices summing close to an explicit ploidy configuration
    mat <- matrix(runif(N * L, min = 0, max = 2), nrow = N, ncol = L)
    # Ensure rows sum up to a constant ploidy level (e.g., tetraploid = 4)
    mat <- t(apply(mat, 1, function(r) r * 4 / sum(r)))
    rownames(mat) <- ind_names
    colnames(mat) <- allele_names
    geno_list[[paste0("Marker", m)]] <- mat
  }
  return(geno_list)
}

# =========================================================================
# 1. Structural Sanity & Convergence Tests
# =========================================================================
test_that("AdmixGlobal converges cleanly and outputs expected attributes", {
  mock_geno <- create_mock_geno()
  K_test <- 2L
  
  # Run a standard inference phase with 1 thread for stable tracing
  res <- AdmixGlobal(Geno = mock_geno, K = K_test, MaxIter = 5L, MinIter = 2L, Verbose = FALSE, NbThreads = 1)
  
  # Structural Inspections
  expect_s3_class(res, "AdmixGlobal")
  expect_named(res, c("Prop", "Freq", "LogLik"))
  
  # Dimension Integrity Mapping
  expect_equal(dim(res$Prop), c(nrow(mock_geno[[1]]), K_test))
  expect_length(res$Freq, length(mock_geno))
  expect_equal(dim(res$Freq[[1]]), c(K_test, ncol(mock_geno[[1]])))
  
  # Mathematical Soundness checks
  expect_equal(rowSums(res$Prop), rep(1.0, nrow(mock_geno[[1]])), tolerance = 1e-6, ignore_attr = TRUE)
  expect_equal(rowSums(res$Freq[[1]]), rep(1.0, K_test), tolerance = 1e-6, ignore_attr = TRUE)
  
  # Verify log-likelihood behavior (Monotonically increasing or converging)
  expect_true(length(res$LogLik) >= 2)
})

# =========================================================================
# 2. Initialization Overrides (Testing Custom Injections)
# =========================================================================
test_that("AdmixGlobal correctly integrates custom Prop and Freq matrices", {
  mock_geno <- create_mock_geno(N = 3, M = 4, L = 2)
  K_test <- 2L
  
  # Construct Valid Pre-baked Initial Arrays
  custom_prop <- matrix(c(0.7, 0.3, 0.5, 0.5, 0.1, 0.9), nrow = 3, ncol = 2, byrow = TRUE)
  rownames(custom_prop) <- rownames(mock_geno[[1]])
  colnames(custom_prop) <- c("K1", "K2")
  
  custom_freq <- list(
    matrix(c(0.6, 0.4, 0.2, 0.8), nrow = 2, ncol = 2, byrow = TRUE),
    matrix(c(0.5, 0.5, 0.5, 0.5), nrow = 2, ncol = 2, byrow = TRUE),
    matrix(c(0.9, 0.1, 0.1, 0.9), nrow = 2, ncol = 2, byrow = TRUE),
    matrix(c(0.3, 0.7, 0.7, 0.3), nrow = 2, ncol = 2, byrow = TRUE)
  )
  
  # Scenario A: Fix Freq, Update Only Prop
  res_prop_only <- AdmixGlobal(
    Geno = mock_geno, K = K_test, 
    PropInit = custom_prop, FreqInit = custom_freq,
    ParamToUpdate = "Prop", MaxIter = 3L, MinIter = 2L, Verbose = FALSE, NbThreads = 1
  )
  # Freq must remain exactly identical to custom input values
  expect_equal(res_prop_only$Freq[[1]], custom_freq[[1]], tolerance = 1e-7, ignore_attr = TRUE)
  
  # Scenario B: Fix Prop, Update Only Freq
  res_freq_only <- AdmixGlobal(
    Geno = mock_geno, K = K_test, 
    PropInit = custom_prop, FreqInit = custom_freq,
    ParamToUpdate = "Freq", MaxIter = 3L, MinIter = 2L, Verbose = FALSE, NbThreads = 1
  )
  # Prop must remain exactly identical to custom input values
  expect_equal(res_freq_only$Prop, custom_prop, tolerance = 1e-7, ignore_attr = TRUE)
  
  # Scenario C: Fix ploidy
  P_scalar <- 6
  res_p_scalar <- AdmixGlobal(
    Geno = mock_geno, K = K_test, P = P_scalar,
    MaxIter = 3L, MinIter = 2L, Verbose = FALSE, NbThreads = 1
  )
  expect_equal(rowSums(res_p_scalar$Prop), rep(1.0, nrow(mock_geno[[1]])), tolerance = 1e-6, ignore_attr = TRUE)
  expect_equal(rowSums(res_p_scalar$Freq[[1]]), rep(1.0, K_test), tolerance = 1e-6, ignore_attr = TRUE)
  P_vect <- 2:4
  res_p_vect <- AdmixGlobal(
    Geno = mock_geno, K = K_test, P = P_vect,
    MaxIter = 3L, MinIter = 2L, Verbose = FALSE, NbThreads = 1
  )
  expect_equal(rowSums(res_p_scalar$Prop), rep(1.0, nrow(mock_geno[[1]])), tolerance = 1e-6, ignore_attr = TRUE)
  expect_equal(rowSums(res_p_scalar$Freq[[1]]), rep(1.0, K_test), tolerance = 1e-6, ignore_attr = TRUE)
})

# =========================================================================
# 3. Input Validation Interceptions (Catching User Errors)
# =========================================================================
test_that("AdmixGlobal raises explicit errors when formatting properties mismatch", {
  mock_geno <- create_mock_geno(N = 3, M = 4, L = 2)
  
  # Unnamed List Inputs
  unnamed_geno <- mock_geno
  names(unnamed_geno) <- NULL
  expect_error(AdmixGlobal(unnamed_geno, K = 2), regexp = "must be a named list")
  
  # Invalid parameter updating keywords
  expect_error(AdmixGlobal(mock_geno, K = 2, ParamToUpdate = "invalid_string"), regexp = "chosen among")
  
  # Iteration Bounds conflicts (MaxIter < MinIter)
  expect_error(AdmixGlobal(mock_geno, K = 2, MaxIter = 5, MinIter = 10), regexp = "greater than or equal to MinIter")
  
  # Mismatched dimensions on incoming tracking components
  bad_prop_shape <- matrix(0.5, nrow = 10, ncol = 2) # Should match Geno N=3
  expect_error(AdmixGlobal(mock_geno, K = 2, PropInit = bad_prop_shape), regexp = "same number of rows")
})

# =========================================================================
# 4. Determinism & Thread Independence 
# =========================================================================
test_that("EM acceleration pipeline yields matching configurations regardless of thread counts", {
  mock_geno <- create_mock_geno(N = 4, M = 10, L = 3)
  K_test <- 2L
  seed_val <- 42L
  
  # Run under single thread processing
  res_serial <- AdmixGlobal(Geno = mock_geno, K = K_test, MaxIter = 6L, MinIter = 2L, Seed = seed_val, NbThreads = 1, Verbose = FALSE)
  
  # Run under dual-thread processing (Forces OpenMP task forks if available)
  res_parallel <- AdmixGlobal(Geno = mock_geno, K = K_test, MaxIter = 6L, MinIter = 2L, Seed = seed_val, NbThreads = 2, Verbose = FALSE)
  
  # Assert structural and optimization parity
  expect_equal(res_serial$Prop, res_parallel$Prop, tolerance = 1e-5)
  expect_equal(res_serial$LogLik, res_parallel$LogLik, tolerance = 1e-5)
})

# =========================================================================
# 5. Silent mode
# =========================================================================
test_that("AdmixGlobal prints absolutely nothing when Verbose is FALSE", {
  mock_geno <- create_mock_geno(N = 3, M = 4, L = 2)
  
  # expect_silent will fail if cat(), print(), message(), or warning() run
  expect_silent(
    AdmixGlobal(
      Geno = mock_geno, 
      K = 2, 
      MaxIter = 2L, 
      MinIter = 2L, 
      NbThreads = 1, 
      Verbose = FALSE # Explicitly turning it off
    )
  )
})