tests/testthat/test-constrained_genomic_indices.R

In selection.index: Analysis of Selection Index in Plant Breeding

# ==============================================================================
# Tests for Constrained Genomic Selection Indices
# ==============================================================================

# Setup: Create test data for constrained genomic indices
setup_constrained_test_data <- function() {
  set.seed(123)

  # Use actual variance-covariance matrices from example data
  gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
  pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

  # Simulate GEBVs using actual phenotypes + noise
  phen_mat <- as.matrix(seldata[1:50, 3:9])
  gebv_mat <- phen_mat * 0.7 + matrix(
    rnorm(prod(dim(phen_mat)), sd = 0.5),
    nrow(phen_mat), ncol(phen_mat)
  )
  colnames(gebv_mat) <- colnames(gmat)

  # GEBV variance-covariance matrix
  Gamma <- cov(gebv_mat)

  n_traits <- ncol(gmat)
  n_genotypes <- nrow(phen_mat)

  # Economic weights
  weights <- weight$ew

  # Create simple constraint matrix (restrict first 2 traits)
  U <- diag(n_traits)[, 1:2, drop = FALSE]

  # Desired proportional gains for PPG methods
  d <- c(2, 1, 1, 0.5, 0.5, 0, 0)

  list(
    gmat = gmat,
    pmat = pmat,
    gebv_mat = gebv_mat,
    phen_mat = phen_mat,
    Gamma = Gamma,
    weights = weights,
    n_traits = n_traits,
    n_genotypes = n_genotypes,
    U = U,
    d = d
  )
}

# ==============================================================================
# TEST: rlgsi() - Restricted Linear Genomic Selection Index
# ==============================================================================

test_that("rlgsi basic functionality with restricted_traits", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- rlgsi(
    Gamma = data$Gamma,
    wmat = data$weights,
    restricted_traits = c(1, 2)
  )

  # Check structure
  expect_s3_class(result, "rlgsi")
  expect_s3_class(result, "genomic_index")
  expect_type(result, "list")

  # Check components exist
  expect_true("b" %in% names(result))
  expect_true("E" %in% names(result))
  expect_true("R" %in% names(result))
  expect_true("GA" %in% names(result))
  expect_true("PRE" %in% names(result))
  expect_true("rHI" %in% names(result))
  expect_true("U" %in% names(result))
  expect_true("constrained_response" %in% names(result))
  expect_true("summary" %in% names(result))

  # Check dimensions
  expect_equal(length(result$b), data$n_traits)
  expect_equal(length(result$E), data$n_traits)
  expect_equal(ncol(result$U), 2) # 2 constraints
  expect_equal(length(result$constrained_response), 2)

  # Check constraints are satisfied (gains should be ~0 for restricted traits)
  # Note: numerical tolerance of 1.0 due to optimization and estimation variability
  expect_true(all(abs(result$constrained_response) < 1.0))

  # Check values are reasonable
  expect_false(is.na(result$R))
  expect_false(is.na(result$GA))
  expect_true(result$rHI >= 0 && result$rHI <= 1)
})

test_that("rlgsi works with custom U matrix", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Custom constraint: restrict last 3 traits
  U_custom <- diag(data$n_traits)[, 5:7, drop = FALSE]

  result <- rlgsi(
    Gamma = data$Gamma,
    wmat = data$weights,
    U = U_custom
  )

  expect_s3_class(result, "rlgsi")
  expect_equal(ncol(result$U), 3)
  expect_equal(length(result$constrained_response), 3)

  # Check constraints are satisfied
  expect_true(all(abs(result$constrained_response) < 1.0))
})

test_that("rlgsi handles single trait restriction", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- rlgsi(
    Gamma = data$Gamma,
    wmat = data$weights,
    restricted_traits = 3
  )

  expect_s3_class(result, "rlgsi")
  expect_equal(ncol(result$U), 1)
  expect_equal(length(result$constrained_response), 1)
  expect_true(abs(result$constrained_response[1]) < 1.0)
})

test_that("rlgsi validates inputs correctly", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Missing both restricted_traits and U
  expect_error(
    rlgsi(Gamma = data$Gamma, wmat = data$weights),
    "Either 'restricted_traits' or 'U' must be provided"
  )

  # Invalid trait indices
  expect_error(
    rlgsi(Gamma = data$Gamma, wmat = data$weights, restricted_traits = c(1, 10)),
    "restricted_traits must be valid trait indices"
  )

  # Non-square Gamma
  expect_error(
    rlgsi(Gamma = data$Gamma[1:5, ], wmat = data$weights, restricted_traits = 1),
    "Gamma must be a square matrix"
  )

  # Mismatched dimensions
  expect_error(
    rlgsi(Gamma = data$Gamma, wmat = data$weights[1:3], restricted_traits = 1),
    "wmat must have"
  )
})

test_that("rlgsi with standardization L_G", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Compute standardization constant
  w <- data$weights
  wGw <- as.numeric(t(w) %*% data$Gamma %*% w)
  L_G <- sqrt(wGw)

  result <- rlgsi(
    Gamma = data$Gamma,
    wmat = data$weights,
    restricted_traits = c(1, 2),
    L_G = L_G
  )

  expect_s3_class(result, "rlgsi")
  expect_equal(result$L_G, L_G)
  expect_false(is.na(result$R))
})

test_that("rlgsi with GAY parameter", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- rlgsi(
    Gamma = data$Gamma,
    wmat = data$weights,
    restricted_traits = c(1, 2),
    GAY = 5.0
  )

  expect_s3_class(result, "rlgsi")
  expect_false(is.na(result$PRE))
  expect_true(is.numeric(result$PRE))
})

# ==============================================================================
# TEST: ppg_lgsi() - Predetermined Proportional Gains LGSI
# ==============================================================================

test_that("ppg_lgsi basic functionality", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- ppg_lgsi(
    Gamma = data$Gamma,
    d = data$d,
    wmat = data$weights
  )

  # Check structure
  expect_s3_class(result, "ppg_lgsi")
  expect_s3_class(result, "genomic_index")
  expect_type(result, "list")

  # Check components exist
  expect_true("b" %in% names(result))
  expect_true("E" %in% names(result))
  expect_true("theta_G" %in% names(result))
  expect_true("gain_ratios" %in% names(result))
  expect_true("constrained_gains" %in% names(result))
  expect_true("desired_gains" %in% names(result))

  # Check dimensions
  expect_equal(length(result$b), data$n_traits)
  expect_equal(length(result$E), data$n_traits)
  expect_equal(length(result$desired_gains), data$n_traits)

  # Check proportionality (gain_ratios should be similar for non-zero d)
  non_zero_d <- data$d != 0
  ratios <- result$gain_ratios[non_zero_d]
  ratios <- ratios[!is.na(ratios)]
  if (length(ratios) > 1) {
    ratio_sd <- sd(ratios)
    ratio_mean <- mean(ratios)
    # Coefficient of variation should be reasonable (gains are approximately proportional)
    # Relaxed tolerance due to estimation variability with real data
    expect_true(ratio_sd / abs(ratio_mean) < 2.0)
  }

  # Check theta_G is numeric
  expect_true(is.numeric(result$theta_G))
  expect_false(is.na(result$theta_G))
})

test_that("ppg_lgsi with custom U matrix", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Apply proportional gains only to first 4 traits
  U_custom <- diag(data$n_traits)[, 1:4, drop = FALSE]
  d_custom <- c(2, 1, 1, 0.5)

  result <- ppg_lgsi(
    Gamma = data$Gamma,
    d = d_custom,
    wmat = data$weights,
    U = U_custom
  )

  expect_s3_class(result, "ppg_lgsi")
  expect_equal(length(result$desired_gains), 4)
  expect_equal(ncol(result$U), 4)
})

test_that("ppg_lgsi without weights", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- ppg_lgsi(
    Gamma = data$Gamma,
    d = data$d
  )

  expect_s3_class(result, "ppg_lgsi")
  # R is computed regardless of weights
  expect_true(is.numeric(result$R))
  # GA and PRE should be NA without weights
  expect_true(is.na(result$GA))
  expect_true(is.na(result$PRE))
})

test_that("ppg_lgsi validates inputs correctly", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Non-square Gamma
  expect_error(
    ppg_lgsi(Gamma = data$Gamma[1:5, ], d = data$d),
    "Gamma must be a square matrix"
  )

  # Mismatched d length (when U not provided)
  expect_error(
    ppg_lgsi(Gamma = data$Gamma, d = c(1, 2, 3)),
    "d must have length"
  )

  # Mismatched d and U dimensions
  U_wrong <- diag(data$n_traits)[, 1:3, drop = FALSE]
  expect_error(
    ppg_lgsi(Gamma = data$Gamma, d = data$d, U = U_wrong),
    "d must have length"
  )
})

test_that("ppg_lgsi with zero desired gains", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # All zero gains (equivalent to RLGSI)
  d_zero <- rep(0, data$n_traits)

  result <- ppg_lgsi(
    Gamma = data$Gamma,
    d = d_zero,
    wmat = data$weights
  )

  expect_s3_class(result, "ppg_lgsi")
  expect_true(is.numeric(result$theta_G))
})

# ==============================================================================
# TEST: crlgsi() - Combined Restricted LGSI
# ==============================================================================

test_that("crlgsi basic functionality with raw data", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- crlgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    wmat = data$weights,
    restricted_traits = c(1, 2),
    reliability = 0.7
  )

  # Check structure
  expect_s3_class(result, "crlgsi")
  expect_s3_class(result, "genomic_index")
  expect_type(result, "list")

  # Check components exist
  expect_true("b" %in% names(result))
  expect_true("b_y" %in% names(result))
  expect_true("b_g" %in% names(result))
  expect_true("E" %in% names(result))
  expect_true("R" %in% names(result))
  expect_true("constrained_response" %in% names(result))

  # Check dimensions
  expect_equal(length(result$b), 2 * data$n_traits) # Combined vector
  expect_equal(length(result$b_y), data$n_traits) # Phenotype coefficients
  expect_equal(length(result$b_g), data$n_traits) # GEBV coefficients
  expect_equal(length(result$E), data$n_traits) # Genetic gains per trait

  # Check constraints are satisfied
  expect_equal(length(result$constrained_response), 2) # 2 restricted traits
  # Combined indices have more numerical complexity, use relaxed tolerance
  expect_true(all(abs(result$constrained_response) < 5.0))

  # Check values are reasonable
  expect_false(is.na(result$R))
  expect_false(is.na(result$GA))
  expect_true(result$rHI >= 0 && result$rHI <= 1)
})

test_that("crlgsi with precomputed T_C and Psi_C matrices", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Compute matrices manually
  P_gebv <- cov(data$gebv_mat)
  P_yg <- cov(data$phen_mat, data$gebv_mat)

  T_C <- rbind(
    cbind(data$pmat, P_yg),
    cbind(t(P_yg), P_gebv)
  )

  reliability <- 0.7
  Gamma_gebv_g <- data$gmat * sqrt(reliability)
  # Psi_C is 2t x 2t: [G, Gamma; Gamma, Gamma]
  Psi_C <- rbind(
    cbind(data$gmat, Gamma_gebv_g),
    cbind(Gamma_gebv_g, Gamma_gebv_g)
  )

  result <- crlgsi(
    T_C = T_C,
    Psi_C = Psi_C,
    wmat = data$weights,
    restricted_traits = c(1, 2)
  )

  expect_s3_class(result, "crlgsi")
  expect_equal(length(result$b), 2 * data$n_traits)
})

test_that("crlgsi with custom U matrix", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Restrict last 3 traits: need 2r = 6 constraints (3 phenotypes + 3 GEBVs)
  # Build 2t x 6 matrix
  r <- 3
  traits_to_restrict <- 5:7
  U_custom <- matrix(0, nrow = 2 * data$n_traits, ncol = 2 * r)
  for (i in seq_along(traits_to_restrict)) {
    trait_idx <- traits_to_restrict[i]
    U_custom[trait_idx, i] <- 1
    U_custom[data$n_traits + trait_idx, r + i] <- 1
  }

  result <- crlgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    wmat = data$weights,
    U = U_custom,
    reliability = 0.7
  )

  expect_s3_class(result, "crlgsi")
  expect_equal(ncol(result$U), 6)
  expect_equal(length(result$constrained_response), 3)
})

test_that("crlgsi with vector reliability", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Different reliability per trait
  rel_vec <- seq(0.5, 0.8, length.out = data$n_traits)

  result <- crlgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    wmat = data$weights,
    restricted_traits = c(1, 2),
    reliability = rel_vec
  )

  expect_s3_class(result, "crlgsi")
  expect_false(is.na(result$R))
})

test_that("crlgsi validates inputs correctly", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  # Missing required matrices
  expect_error(
    crlgsi(wmat = data$weights, restricted_traits = 1),
    "Must provide either"
  )

  # Missing both restricted_traits and U
  expect_error(
    crlgsi(
      phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
      pmat = data$pmat, gmat = data$gmat, wmat = data$weights
    ),
    "Either 'restricted_traits' or 'U' must be provided"
  )

  # Missing gmat when using raw data
  expect_error(
    crlgsi(
      phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
      pmat = data$pmat, wmat = data$weights, restricted_traits = 1
    ),
    "gmat required"
  )
})

test_that("crlgsi summary format", {
  skip_on_cran() # heavy cross-products / TRE regex — bypass CRAN sanitizers
  data <- setup_constrained_test_data()

  result <- crlgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    wmat = data$weights,
    restricted_traits = c(1, 2),
    reliability = 0.7
  )

  # Check summary data frame
  expect_s3_class(result$summary, "data.frame")
  expect_equal(nrow(result$summary), 1)

  # Check summary has correct columns
  expect_true(any(grepl("^b_y\\.", colnames(result$summary), perl = TRUE)))
  expect_true(any(grepl("^b_g\\.", colnames(result$summary), perl = TRUE)))
  expect_true("R" %in% colnames(result$summary))
  expect_true("GA" %in% colnames(result$summary))
  expect_true("rHI" %in% colnames(result$summary))
})

# ==============================================================================
# TEST: cppg_lgsi() - Combined Predetermined Proportional Gains LGSI
# ==============================================================================

test_that("cppg_lgsi basic functionality with raw data", {
  data <- setup_constrained_test_data()

  result <- cppg_lgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    d = data$d,
    wmat = data$weights,
    reliability = 0.7
  )

  # Check structure
  expect_s3_class(result, "cppg_lgsi")
  expect_s3_class(result, "genomic_index")
  expect_type(result, "list")

  # Check components exist
  expect_true("b" %in% names(result))
  expect_true("b_y" %in% names(result))
  expect_true("b_g" %in% names(result))
  expect_true("E" %in% names(result))
  expect_true("theta_CP" %in% names(result))
  expect_true("gain_ratios" %in% names(result))

  # Check dimensions
  expect_equal(length(result$b), 2 * data$n_traits)
  expect_equal(length(result$b_y), data$n_traits)
  expect_equal(length(result$b_g), data$n_traits)
  expect_equal(length(result$desired_gains), data$n_traits)

  # Check proportionality
  non_zero_d <- data$d != 0
  ratios <- result$gain_ratios[non_zero_d]
  ratios <- ratios[!is.na(ratios)]
  if (length(ratios) > 1) {
    ratio_sd <- sd(ratios)
    ratio_mean <- mean(ratios)
    expect_true(ratio_sd / abs(ratio_mean) < 2.0)
  }

  # Check theta_CP
  expect_true(is.numeric(result$theta_CP))
  expect_false(is.na(result$theta_CP))
})

test_that("cppg_lgsi with precomputed matrices", {
  data <- setup_constrained_test_data()

  # Compute matrices manually
  P_gebv <- cov(data$gebv_mat)
  P_yg <- cov(data$phen_mat, data$gebv_mat)

  T_C <- rbind(
    cbind(data$pmat, P_yg),
    cbind(t(P_yg), P_gebv)
  )

  reliability <- 0.7
  Gamma_gebv_g <- data$gmat * sqrt(reliability)
  # Psi_C is 2t x 2t
  Psi_C <- rbind(
    cbind(data$gmat, Gamma_gebv_g),
    cbind(Gamma_gebv_g, Gamma_gebv_g)
  )

  result <- cppg_lgsi(
    T_C = T_C,
    Psi_C = Psi_C,
    d = data$d,
    wmat = data$weights
  )

  expect_s3_class(result, "cppg_lgsi")
  expect_equal(length(result$b), 2 * data$n_traits)
})

test_that("cppg_lgsi with custom U matrix", {
  data <- setup_constrained_test_data()

  # Apply proportional gains only to first 4 traits
  # Need 2t x 2r matrix (14 x 8)
  r <- 4
  traits_to_constrain <- 1:4
  U_custom <- matrix(0, nrow = 2 * data$n_traits, ncol = 2 * r)
  for (i in seq_along(traits_to_constrain)) {
    trait_idx <- traits_to_constrain[i]
    U_custom[trait_idx, i] <- 1
    U_custom[data$n_traits + trait_idx, r + i] <- 1
  }
  d_custom <- c(2, 1, 1, 0.5)

  result <- cppg_lgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    d = d_custom,
    wmat = data$weights,
    U = U_custom,
    reliability = 0.7
  )

  expect_s3_class(result, "cppg_lgsi")
  expect_equal(length(result$desired_gains), 4)
  expect_equal(ncol(result$U), 8)
})

test_that("cppg_lgsi without weights", {
  data <- setup_constrained_test_data()

  result <- cppg_lgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    d = data$d,
    reliability = 0.7
  )

  expect_s3_class(result, "cppg_lgsi")
  expect_false(is.na(result$R))
})

test_that("cppg_lgsi validates inputs correctly", {
  skip_on_cran() # error handling test or warning test
  data <- setup_constrained_test_data()

  # Missing d parameter
  expect_error(
    cppg_lgsi(
      phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
      pmat = data$pmat, gmat = data$gmat, wmat = data$weights,
      reliability = 0.7
    ),
    "argument \"d\" is missing"
  )

  # Missing required matrices
  expect_error(
    cppg_lgsi(d = data$d, wmat = data$weights),
    "Must provide either"
  )

  # Mismatched d and U dimensions
  U_wrong <- matrix(0, nrow = 2 * data$n_traits, ncol = 4) # Too few constraints
  expect_error(
    cppg_lgsi(
      phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
      pmat = data$pmat, gmat = data$gmat, d = data$d,
      wmat = data$weights, U = U_wrong, reliability = 0.7
    ),
    "d must have length"
  )
})

test_that("cppg_lgsi summary format", {
  data <- setup_constrained_test_data()

  result <- cppg_lgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    d = data$d,
    wmat = data$weights,
    reliability = 0.7
  )

  # Check summary data frame
  expect_s3_class(result$summary, "data.frame")
  expect_equal(nrow(result$summary), 1)

  # Check summary has correct columns
  expect_true(any(grepl("^b_y\\.", colnames(result$summary), perl = TRUE)))
  expect_true(any(grepl("^b_g\\.", colnames(result$summary), perl = TRUE)))
  expect_true("R" %in% colnames(result$summary))
  expect_true("theta_CP" %in% colnames(result$summary))
})

# ==============================================================================
# INTEGRATION TESTS: Compare methods
# ==============================================================================

test_that("RLGSI vs PPG-LGSI with zero gains", {
  data <- setup_constrained_test_data()

  # RLGSI with restrictions
  rlgsi_result <- rlgsi(
    Gamma = data$Gamma,
    wmat = data$weights,
    restricted_traits = c(1, 2)
  )

  # PPG-LGSI with zero desired gains (should be similar to RLGSI)
  d_restricted <- data$d
  d_restricted[1:2] <- 0

  ppg_result <- ppg_lgsi(
    Gamma = data$Gamma,
    d = d_restricted,
    wmat = data$weights
  )

  # Both should have near-zero gains for traits 1 and 2
  expect_true(all(abs(rlgsi_result$constrained_response) < 1.0))
  expect_true(all(abs(ppg_result$constrained_gains[1:2]) < 1.0))
})

test_that("CRLGSI vs CPPG-LGSI with zero gains", {
  data <- setup_constrained_test_data()

  # CRLGSI with restrictions
  crlgsi_result <- crlgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    wmat = data$weights,
    restricted_traits = c(1, 2),
    reliability = 0.7
  )

  # CPPG-LGSI with zero desired gains
  d_restricted <- data$d
  d_restricted[1:2] <- 0

  cppg_result <- cppg_lgsi(
    phen_mat = data$phen_mat,
    gebv_mat = data$gebv_mat,
    pmat = data$pmat,
    gmat = data$gmat,
    d = d_restricted,
    wmat = data$weights,
    reliability = 0.7
  )

  # Both should have reduced gains for traits 1 and 2 (constraints applied)
  # Combined indices have more numerical complexity
  expect_true(all(abs(crlgsi_result$constrained_response) < 5.0))
  expect_true(all(abs(cppg_result$constrained_gains[1:2]) < 5.0))
})

# ==============================================================================
# NUMERICAL STABILITY TESTS
# ==============================================================================

test_that("All methods handle near-singular matrices", {
  skip_on_cran() # error handling test or warning test
  # Create a nearly singular Gamma matrix
  set.seed(456)
  n_traits <- 5
  Gamma_singular <- matrix(rnorm(n_traits^2, sd = 0.01), n_traits, n_traits)
  Gamma_singular <- (Gamma_singular + t(Gamma_singular)) / 2
  diag(Gamma_singular) <- 1

  w <- rep(1, n_traits)
  d <- c(1, 1, 0, 0, 0)

  # Test RLGSI
  expect_error(
    rlgsi(Gamma = Gamma_singular, wmat = w, restricted_traits = 1),
    NA # Should not error (using ginv for robustness)
  )

  # Test PPG-LGSI
  expect_error(
    ppg_lgsi(Gamma = Gamma_singular, d = d, wmat = w),
    NA
  )
})

test_that("All methods handle small values correctly", {
  data <- setup_constrained_test_data()

  # Scale down Gamma (small variances)
  Gamma_small <- data$Gamma * 0.001

  result <- rlgsi(
    Gamma = Gamma_small,
    wmat = data$weights,
    restricted_traits = c(1, 2)
  )

  expect_s3_class(result, "rlgsi")
  expect_false(any(is.na(result$b)))
  expect_false(any(is.infinite(result$b)))
})

# ==============================================================================
# EDGE CASES AND ERROR CONDITIONS FOR 100% COVERAGE
# ==============================================================================

test_that("genomic_index_metrics handles NA/Inf bGb (line 106)", {
  data <- setup_constrained_test_data()
  T_C_zero <- matrix(0, 2 * data$n_traits, 2 * data$n_traits)
  Psi_C_zero <- matrix(0, 2 * data$n_traits, 2 * data$n_traits)

  result <- crlgsi(T_C = T_C_zero, Psi_C = Psi_C_zero, wmat = data$weights, restricted_traits = 1)
  expect_true(all(is.na(result$E)))
})

test_that("rlgsi edge cases (lines 258, 292)", {
  skip_on_cran() # error handling test or warning test
  data <- setup_constrained_test_data()

  U_bad <- matrix(0, nrow = data$n_traits - 1, ncol = 1)
  expect_error(
    rlgsi(Gamma = data$Gamma, wmat = data$weights, U = U_bad),
    "U must have .* rows"
  )

  wmat_na <- data$weights
  wmat_na[1] <- NA
  expect_error(
    rlgsi(Gamma = data$Gamma, wmat = wmat_na, restricted_traits = 1),
    "RLGSI coefficients contain NA or Inf"
  )
})

test_that("ppg_lgsi edge cases (lines 450, 462)", {
  skip_on_cran() # error handling test or warning test
  data <- setup_constrained_test_data()

  U_bad <- matrix(0, nrow = data$n_traits - 1, ncol = 1)
  expect_error(
    ppg_lgsi(Gamma = data$Gamma, d = data$d, U = U_bad),
    "U must have .* rows"
  )

  wmat_bad <- data$weights[1:(data$n_traits - 1)]
  expect_error(
    ppg_lgsi(Gamma = data$Gamma, d = data$d, wmat = wmat_bad),
    "wmat must have .* rows"
  )
})

test_that("crlgsi edge cases (lines 685, 723, 727, 733, 740, 794)", {
  skip_on_cran() # error handling test or warning test
  data <- setup_constrained_test_data()

  result <- crlgsi(
    phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
    gmat = data$gmat, wmat = data$weights, restricted_traits = 1,
    reliability = 0.7
  )
  expect_s3_class(result, "crlgsi")

  T_C_bad <- matrix(0, nrow = 2 * data$n_traits - 1, ncol = 2 * data$n_traits - 1)
  expect_error(
    crlgsi(T_C = T_C_bad, Psi_C = diag(2 * data$n_traits), wmat = data$weights, restricted_traits = 1),
    "T_C must be \\(2\\*n_traits x 2\\*n_traits\\)"
  )

  T_C_good <- matrix(0, nrow = 2 * data$n_traits, ncol = 2 * data$n_traits)
  Psi_C_bad <- matrix(0, nrow = 2 * data$n_traits - 1, ncol = 2 * data$n_traits)
  expect_error(
    crlgsi(T_C = T_C_good, Psi_C = Psi_C_bad, wmat = data$weights, restricted_traits = 1),
    "Psi_C must be \\(2\\*n_traits x 2\\*n_traits\\)"
  )

  wmat_bad <- data$weights[1:(data$n_traits - 1)]
  expect_error(
    crlgsi(T_C = T_C_good, Psi_C = T_C_good, wmat = wmat_bad, restricted_traits = 1),
    "wmat must have .* rows"
  )

  expect_error(
    crlgsi(T_C = T_C_good, Psi_C = T_C_good, wmat = data$weights, restricted_traits = data$n_traits + 1),
    "restricted_traits must be valid trait indices"
  )

  wmat_na <- data$weights
  wmat_na[1] <- NA
  expect_error(
    crlgsi(T_C = T_C_good, Psi_C = T_C_good, wmat = wmat_na, restricted_traits = 1),
    "CRLGSI coefficients contain NA or Inf"
  )
})

test_that("cppg_lgsi edge cases (lines 970, 971, 985, 1005, 1021, 1031, 1042)", {
  skip_on_cran() # error handling test or warning test
  data <- setup_constrained_test_data()

  result <- cppg_lgsi(
    phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
    gmat = data$gmat, d = data$d, wmat = data$weights
  )
  expect_s3_class(result, "cppg_lgsi")

  expect_error(
    cppg_lgsi(
      phen_mat = data$phen_mat, gebv_mat = data$gebv_mat,
      d = data$d, wmat = data$weights
    ),
    "gmat required when using raw data"
  )

  T_C_good <- diag(2 * data$n_traits)
  Psi_C_good <- diag(2 * data$n_traits)

  d_bad <- c(1, 2)
  expect_error(
    cppg_lgsi(T_C = T_C_good, Psi_C = Psi_C_good, d = d_bad, wmat = data$weights),
    "d must have length .* when U is not provided"
  )

  U_bad <- matrix(0, nrow = 2 * data$n_traits - 1, ncol = 2)
  expect_error(
    cppg_lgsi(T_C = T_C_good, Psi_C = Psi_C_good, d = data$d, U = U_bad, wmat = data$weights),
    "U must have .* rows for combined indices"
  )

  U_custom <- diag(2 * data$n_traits)
  d_C_full <- rep(1, 2 * data$n_traits)
  res_custom <- cppg_lgsi(T_C = T_C_good, Psi_C = Psi_C_good, d = d_C_full, U = U_custom, wmat = data$weights)
  expect_s3_class(res_custom, "cppg_lgsi")

  wmat_bad <- data$weights[1:(data$n_traits - 1)]
  expect_error(
    cppg_lgsi(T_C = T_C_good, Psi_C = Psi_C_good, d = data$d, wmat = wmat_bad),
    "wmat must have .* rows"
  )
})