tests/testthat/test-admix-local.R

In AdmixPoly: Global and Local Admixture Inference in Polyploids

# =========================================================================
# Shared Setup Helper: Minimal inputs for fast testing
# =========================================================================
create_local_mock_data <- function() {
  set.seed(123)
  
  # 1. Mock Genotypes (3 Individuals, 6 Markers, 2 Alleles)
  markers <- paste0("M", 1:6)
  inds <- c("IndA", "IndB", "PureInd")
  alleles <- c("A1", "A2")
  
  Geno <- list()
  for(m in markers) {
    mat <- matrix(c(2,2, 0,4, 4,0), nrow = 3, ncol = 2, byrow = TRUE)
    rownames(mat) <- inds
    colnames(mat) <- alleles
    Geno[[m]] <- mat
  }
  
  # 2. Mock Genetic Map (Split across 2 Chromosomes)
  GeneticMap <- data.frame(
    Marker = markers,
    Chromosome = c("Chr1", "Chr1", "Chr1", "Chr2", "Chr2", "Chr2"),
    Distance = c(0.0, 0.5, 1.2, 0.0, 2.1, 4.5),
    stringsAsFactors = FALSE
  )
  
  # 3. Mock ResAdmixGlobal S3 Object
  Prop <- matrix(c(
    0.60, 0.40,  # IndA: Admixed
    0.50, 0.50,  # IndB: Admixed
    1.00, 0.00   # PureInd: Non-Admixed (Below MinProp)
  ), nrow = 3, ncol = 2, byrow = TRUE)
  rownames(Prop) <- inds
  colnames(Prop) <- c("K1", "K2")
  
  Freq <- list()
  for(m in markers) {
    f_mat <- matrix(c(0.7, 0.3, 0.2, 0.8), nrow = 2, ncol = 2, byrow = TRUE)
    rownames(f_mat) <- c("K1", "K2")
    colnames(f_mat) <- alleles
    Freq[[m]] <- f_mat
  }
  
  ResAdmixGlobal <- list(Prop = Prop, Freq = Freq, LogLik = c(-100, -50))
  class(ResAdmixGlobal) <- "AdmixGlobal"
  
  # Mock the internal C++ link pointer outputs if devtools::load_all() overrides aren't hooked
  # We mock .ForwardBackward dynamically or let it evaluate if compiled
  # if(!exists(".ForwardBackward")) {
  #   options(mock_cpp = TRUE) 
  # }
  
  return(list(Geno = Geno, ResAdmixGlobal = ResAdmixGlobal, GeneticMap = GeneticMap))
}

# =========================================================================
# 1. Non-Admixed Short-Circuit Validation (Early Returns)
# =========================================================================
test_that("AdmixLocal handles non-admixed individuals with an early short-circuit", {
  setup <- create_local_mock_data()
  P_test <- 4L
  
  res <- AdmixLocal(
    Geno = setup$Geno, 
    ResAdmixGlobal = setup$ResAdmixGlobal, 
    Ind = "PureInd", # 100% K1, 0% K2
    P = P_test, 
    GeneticMap = setup$GeneticMap, 
    MinProp = 1e-3, 
    Verbose = FALSE,
    NbThreads = 1
  )
  
  # Check early list output footprint
  expect_named(res, c("Posterior", "Viterbi"))
  expect_null(res$LogLik) # Early return does not compute HMM log-likelihood
  
  # Verify that all slots for the non-admixed group are pinned to total Ploidy (P)
  expect_equal(res$Posterior$Chr1["M1", "K1"], P_test)
  expect_equal(res$Posterior$Chr1["M1", "K2"], 0)
})

# =========================================================================
# 2. HMM Pipeline Run & Structural Verification (The Happy Path)
# =========================================================================
test_that("AdmixLocal executes full HMM pipeline and resolves structures cleanly", {
  # If local C++ hooks exist, test structural properties
  setup <- create_local_mock_data()
  
  # Deactivate distance approximations to evaluate exact tracking branches
  res <- AdmixLocal(
    Geno = setup$Geno, 
    ResAdmixGlobal = setup$ResAdmixGlobal, 
    Ind = "IndA", 
    P = 4L, 
    GeneticMap = setup$GeneticMap, 
    DistIntBounds = NULL, 
    MaxRec = NULL, 
    Verbose = FALSE,
    NbThreads = 1
  )
  
  expect_type(res, "list")
  expect_true(all(c("Posterior", "Viterbi", "LogLik", "EmisProbType") %in% names(res)))
  
  # Chromosome List Split mapping
  expect_length(res$Posterior, 2) # Chr1 and Chr2
  expect_length(res$Viterbi, 2) # Chr1 and Chr2
  expect_equal(names(res$Posterior), c("Chr1", "Chr2"))
  expect_equal(names(res$Viterbi), c("Chr1", "Chr2"))
  
  # Check matrix dimensions inside target splits
  expect_equal(dim(res$Posterior$Chr1), c(3, 2)) # 3 markers on Chr1, 2 ancestral populations
  expect_equal(dim(res$Viterbi$Chr1), c(3, 2)) # 3 markers on Chr1, 2 ancestral populations
  
  # Check that ancestry dosages sum to ploidy
  expect_equal(c(sapply(res$Posterior,rowSums)),matrix(4,3,2), tolerance = 1e-6, ignore_attr = TRUE)
})

# =========================================================================
# 3. Distance Simplification & Parameter Variation Branches
# =========================================================================
test_that("AdmixLocal accommodates Distance Simplification and Recombination constraints", {
  setup <- create_local_mock_data()
  
  # Scenario A: Enable Interval Boundaries Approximations
  expect_error(
    AdmixLocal(
      Geno = setup$Geno, ResAdmixGlobal = setup$ResAdmixGlobal, 
      Ind = "IndA", P = 4L, GeneticMap = setup$GeneticMap, 
      DistIntBounds = c(0.1, 0.5, 1.0), # Custom bounds vector
      Verbose = FALSE, NbThreads = 1
    ), 
    NA # Asserts that no error or crash occurs
  )
  
  # Scenario B: Enable Maximum Recombination Bounds cap
  expect_error(
    AdmixLocal(
      Geno = setup$Geno, ResAdmixGlobal = setup$ResAdmixGlobal, 
      Ind = "IndA", P = 4L, GeneticMap = setup$GeneticMap, 
      MaxRec = 2L, # Must be between 1 and P
      Verbose = FALSE, NbThreads = 1
    ),
    NA
  )
})

# =========================================================================
# 4. Equivalence of MaxRec trick when equal to the ploidy
# =========================================================================
test_that("Exact and approximation using MaxRec=P are equivalent", {
  setup <- create_local_mock_data()
  Ptest <- 4L
  
  # Scenario A: Exact transition probabilities
  res_exact <- AdmixLocal(
    Geno = setup$Geno, 
    ResAdmixGlobal = setup$ResAdmixGlobal, 
    Ind = "IndA", 
    P = Ptest, 
    GeneticMap = setup$GeneticMap, 
    DistIntBounds = NULL, 
    MaxRec = NULL, 
    Verbose = FALSE,
    NbThreads = 1
  )
  
  # Scenario B: MaxRec=P
  res_maxrec <- AdmixLocal(
    Geno = setup$Geno, 
    ResAdmixGlobal = setup$ResAdmixGlobal, 
    Ind = "IndA", 
    P = Ptest, 
    GeneticMap = setup$GeneticMap, 
    DistIntBounds = NULL, 
    MaxRec = Ptest, 
    Verbose = FALSE,
    NbThreads = 1
  )
  
  # Check that ancestry dosages are equal
  expect_equal(unlist(res_exact$Posterior),unlist(res_maxrec$Posterior), tolerance = 1e-6, ignore_attr = TRUE)
  expect_equal(unlist(res_exact$Viterbi),unlist(res_maxrec$Viterbi), tolerance = 1e-6, ignore_attr = TRUE)
  
})

# =========================================================================
# 5. Input Constraints & Defensive Failures
# =========================================================================
test_that("AdmixLocal handles structural violations defensively", {
  setup <- create_local_mock_data()
  
  # Requesting an Individual absent from global matrices
  expect_error(
    AdmixLocal(setup$Geno, setup$ResAdmixGlobal, Ind = "MissingInd", P = 4L, setup$GeneticMap),
    regexp = "included in the dataset"
  )
  
  # Passing an un-associated object structure to ResAdmixGlobal
  bad_global <- list(Data = 123)
  expect_error(
    AdmixLocal(setup$Geno, bad_global, Ind = "IndA", P = 4L, setup$GeneticMap),
    regexp = "generated by the AdmixGlobal function"
  )
  
  # Providing an out-of-bounds MaxRec configuration (MaxRec > P)
  expect_error(
    AdmixLocal(setup$Geno, setup$ResAdmixGlobal, Ind = "IndA", P = 2L, setup$GeneticMap, MaxRec = 5L),
    regexp = "range between 1 and P"
  )
})

# =========================================================================
# 6. Determinism & Thread Independence 
# =========================================================================
test_that("EM acceleration pipeline yields matching configurations regardless of thread counts", {
  setup <- create_local_mock_data()
  
  # Run under single thread processing
  res_serial <- AdmixLocal(
    Geno = setup$Geno,
    ResAdmixGlobal = setup$ResAdmixGlobal, 
    Ind = "IndA", 
    P = 4L, 
    GeneticMap = setup$GeneticMap, 
    Verbose = FALSE,
    NbThreads = 1
  )
  
  # Run under dual-thread processing (Forces OpenMP task forks if available)
  res_parallel <- AdmixLocal(
    Geno = setup$Geno,
    ResAdmixGlobal = setup$ResAdmixGlobal, 
    Ind = "IndA", 
    P = 4L, 
    GeneticMap = setup$GeneticMap, 
    Verbose = FALSE,
    NbThreads = 2
  )
  
  # Assert structural and optimization parity
  expect_equal(unlist(res_serial$Posterior),unlist(res_parallel$Posterior), tolerance = 1e-6, ignore_attr = TRUE)
  expect_equal(unlist(res_serial$Viterbi),unlist(res_parallel$Viterbi), tolerance = 1e-6, ignore_attr = TRUE)
})

# =========================================================================
# 7. Output Suppression (Silence Validation)
# =========================================================================
test_that("AdmixLocal maintains total silence when Verbose is set to FALSE", {
  setup <- create_local_mock_data()
  
  expect_silent(
    AdmixLocal(
      Geno = setup$Geno, 
      ResAdmixGlobal = setup$ResAdmixGlobal, 
      Ind = "PureInd", 
      P = 4L, 
      GeneticMap = setup$GeneticMap, 
      Verbose = FALSE,
      NbThreads = 1
    )
  )
})