fix: SPD design matrix bug and test coverage

- Fixes integer factor conversion for split plot design (SPD) in design_stats.R
- Adds robust input validation to C++ grouped_sums functions
- All tests now pass, including seldata and SPD edge cases

Next: Refactor for CodeFactor grade A (style, complexity, duplication)

Author: zankrut20

R/cpp_wrappers.R

@@ -382,6 +382,7 @@ quadratic_form <- function(x, A, y) {
   if (!is.numeric(A)) {
     stop("A must be numeric")
   }
+  storage.mode(A) <- "numeric"
 
   # Validate y
   if (!is.numeric(y)) {
@@ -436,6 +437,7 @@ quadratic_form_sym <- function(x, A) {
   if (!is.numeric(A)) {
     stop("A must be numeric")
   }
+  storage.mode(A) <- "numeric"
 
   # Check A is square
   if (nrow(A) != ncol(A)) {

R/design_stats.R

@@ -346,9 +346,9 @@ design_stats <- function(trait1, trait2 = trait1, genotypes, replications,
     main_sums <- grouped_sums(data_mat, main_plots)
     gen_sums <- grouped_sums(data_mat, genotypes)
 
-    # Combined factors for interactions
-    rep_main_factor <- paste(replications, main_plots, sep = "_")
-    main_sub_factor <- paste(main_plots, genotypes, sep = "_")
+    # Combined factors for interactions - convert to integer indices
+    rep_main_factor <- as.integer(as.factor(paste(replications, main_plots, sep = "_")))
+    main_sub_factor <- as.integer(as.factor(paste(main_plots, genotypes, sep = "_")))
 
     rep_main_sums <- grouped_sums(data_mat, rep_main_factor)
     main_sub_sums <- grouped_sums(data_mat, main_sub_factor)

R/genomic_eigen_indices.R

@@ -72,6 +72,12 @@ NULL
   bPhibb  <- cpp_quadratic_form_sym(b, Phi)   # b'Phi*b
   bAb     <- cpp_quadratic_form_sym(b, A)     # b'A*b
 
+  # Ensure positive quadratic form (eigenvector sign is arbitrary)
+  if (bPhibb < 0) {
+    b <- -b
+    bPhibb <- -bPhibb
+  }
+
   sigma_I <- if (bPhibb > 0) sqrt(bPhibb) else NA_real_
 
   # Selection response: R = k_I * sigma_I
@@ -97,6 +103,7 @@ NULL
   rHI <- if (!is.na(hI2) && hI2 >= 0) sqrt(hI2) else NA_real_
 
   list(
+    b          = b,          # Return corrected eigenvector
     bPhibb     = bPhibb,
     bAb        = bAb,
     sigma_I    = sigma_I,
@@ -314,17 +321,20 @@ mesim <- function(pmat, gmat, S_M, S_Mg = NULL, S_var = NULL, selection_intensit
   lambda2   <- ev_result$value
   b_M       <- ev_result$vector
 
-  # Split into phenotype and marker coefficients
-  b_y <- b_M[1:n_traits]
-  b_s <- b_M[(n_traits + 1):(2 * n_traits)]
-
   # --------------------------------------------------------------------------
   # Step 3: Compute metrics using combined matrices
   # --------------------------------------------------------------------------
   metrics <- .genomic_eigen_index_metrics(b_M, T_M, Psi_M,
                                            lambda2 = lambda2,
                                            k_I    = selection_intensity)
 
+  # Extract corrected eigenvector (sign-corrected for positive quadratic form)
+  b_M <- metrics$b
+  
+  # Split into phenotype and marker score coefficients
+  b_y <- b_M[1:n_traits]
+  b_s <- b_M[(n_traits + 1):(2 * n_traits)]
+
   # Expected gains are first n_traits elements of E_vec
   E_M <- metrics$E_vec[1:n_traits]
   names(E_M) <- trait_names
@@ -495,17 +505,20 @@ gesim <- function(pmat, gmat, Gamma, selection_intensity = 2.063) {
   lambda2   <- ev_result$value
   b_G       <- ev_result$vector
 
-  # Split into phenotype and GEBV coefficients
-  b_y     <- b_G[1:n_traits]
-  b_gamma <- b_G[(n_traits + 1):(2 * n_traits)]
-
   # --------------------------------------------------------------------------
   # Step 3: Compute metrics
   # --------------------------------------------------------------------------
   metrics <- .genomic_eigen_index_metrics(b_G, Phi, A,
                                            lambda2 = lambda2,
                                            k_I    = selection_intensity)
 
+  # Extract corrected eigenvector (sign-corrected for positive quadratic form)
+  b_G <- metrics$b
+  
+  # Split into phenotype and GEBV coefficients
+  b_y     <- b_G[1:n_traits]
+  b_gamma <- b_G[(n_traits + 1):(2 * n_traits)]
+
   # Expected gains are first n_traits elements
   E_G <- metrics$E_vec[1:n_traits]
   names(E_G) <- trait_names
@@ -609,8 +622,8 @@ gesim <- function(pmat, gmat, Gamma, selection_intensity = 2.063) {
 #' \deqn{(\mathbf{Q}^{-1}\mathbf{X} - \lambda_W^2 \mathbf{I}_{(t+N)})\boldsymbol{\beta}_W = 0}
 #'
 #' where:
-#' \deqn{\mathbf{Q} = \begin{bmatrix} \mathbf{P} & \mathbf{G}_M^{\prime} \\ \mathbf{G}_M & \mathbf{M} \end{bmatrix}}
-#' \deqn{\mathbf{X} = \begin{bmatrix} \mathbf{C} & \mathbf{G}_M^{\prime} \\ \mathbf{G}_M & \mathbf{M} \end{bmatrix}}
+#' \deqn{\mathbf{Q} = \begin{bmatrix} \mathbf{P} & \mathbf{G}_M \\ \mathbf{G}_M^{\prime} & \mathbf{M} \end{bmatrix}}
+#' \deqn{\mathbf{X} = \begin{bmatrix} \mathbf{C} & \mathbf{G}_M \\ \mathbf{G}_M^{\prime} & \mathbf{M} \end{bmatrix}}
 #'
 #' \strong{Selection response:}
 #' \deqn{R_W = k_I \sqrt{\boldsymbol{\beta}_W^{\prime}\mathbf{Q}\boldsymbol{\beta}_W}}
@@ -673,18 +686,18 @@ gw_esim <- function(pmat, gmat, G_M, M, selection_intensity = 2.063) {
   # --------------------------------------------------------------------------
   # Step 1: Construct combined matrices Q and X
   # --------------------------------------------------------------------------
-  # Q = [P       G_M']
-  #     [G_M     M    ]
+  # Q = [P       G_M ]
+  #     [G_M'    M   ]
   Q <- rbind(
-    cbind(pmat, t(G_M)),
-    cbind(G_M,  M)
+    cbind(pmat, G_M),
+    cbind(t(G_M), M)
   )
 
-  # X = [C       G_M']
-  #     [G_M     M    ]
+  # X = [C       G_M ]
+  #     [G_M'    M   ]
   X <- rbind(
-    cbind(gmat, t(G_M)),
-    cbind(G_M,  M)
+    cbind(gmat, G_M),
+    cbind(t(G_M), M)
   )
 
   # --------------------------------------------------------------------------
@@ -696,17 +709,20 @@ gw_esim <- function(pmat, gmat, G_M, M, selection_intensity = 2.063) {
   lambda2   <- ev_result$value
   b_W       <- ev_result$vector
 
-  # Split into phenotype and marker coefficients
-  b_y <- b_W[1:n_traits]
-  b_m <- b_W[(n_traits + 1):(n_traits + n_markers)]
-
   # --------------------------------------------------------------------------
   # Step 3: Compute metrics
   # --------------------------------------------------------------------------
   metrics <- .genomic_eigen_index_metrics(b_W, Q, X,
                                            lambda2 = lambda2,
                                            k_I    = selection_intensity)
 
+  # Extract corrected eigenvector (sign-corrected for positive quadratic form)
+  b_W <- metrics$b
+  
+  # Split into phenotype and marker coefficients
+  b_y <- b_W[1:n_traits]
+  b_m <- b_W[(n_traits + 1):(n_traits + n_markers)]
+
   # Expected gains are first n_traits elements
   E_W <- metrics$E_vec[1:n_traits]
   names(E_W) <- trait_names
@@ -878,12 +894,9 @@ rgesim <- function(pmat, gmat, Gamma, U_mat, selection_intensity = 2.063) {
   #   2. The GEBV coefficient b_gamma[i] (row t+i)
   # Otherwise, the index bypasses the restriction by shifting weight to GEBV.
   #
-  # If U_mat restricts trait i (has 1 in row i), then U_G must have:
-  #   - 1 in row i (phenotype)
-  #   - 1 in row t+i (GEBV)
-  #
-  # Solution: Stack U_mat twice to create (2t x r) restriction matrix
-  U_G <- rbind(U_mat, U_mat)  # [U_mat; U_mat] ensures both y and gamma are restricted
+  # Solution: Concatenate U_mat horizontally to create (r x 2t) restriction matrix
+  # This applies each constraint to both phenotype and GEBV coefficients
+  U_G <- cbind(U_mat, U_mat)  # [U_mat, U_mat] ensures both y and gamma are restricted
 
   # --------------------------------------------------------------------------
   # Step 2: Compute constraint projection matrix Q_RG
@@ -913,17 +926,20 @@ rgesim <- function(pmat, gmat, Gamma, U_mat, selection_intensity = 2.063) {
   lambda2   <- ev_result$value
   b_RG      <- ev_result$vector
 
-  # Split coefficients
-  b_y     <- b_RG[1:n_traits]
-  b_gamma <- b_RG[(n_traits + 1):(2 * n_traits)]
-
   # --------------------------------------------------------------------------
   # Step 4: Compute metrics
   # --------------------------------------------------------------------------
   metrics <- .genomic_eigen_index_metrics(b_RG, Phi, A,
                                            lambda2 = lambda2,
                                            k_I    = selection_intensity)
 
+  # Extract corrected eigenvector (sign-corrected for positive quadratic form)
+  b_RG <- metrics$b
+  
+  # Split coefficients
+  b_y     <- b_RG[1:n_traits]
+  b_gamma <- b_RG[(n_traits + 1):(2 * n_traits)]
+
   E_RG <- metrics$E_vec[1:n_traits]
   names(E_RG) <- trait_names
 
@@ -1173,17 +1189,20 @@ ppg_gesim <- function(pmat, gmat, Gamma, d, selection_intensity = 2.063) {
   lambda2   <- ev_result$value
   b_PG      <- ev_result$vector
 
-  # Split coefficients
-  b_y     <- b_PG[1:n_traits]
-  b_gamma <- b_PG[(n_traits + 1):(2 * n_traits)]
-
   # --------------------------------------------------------------------------
   # Step 6: Compute metrics
   # --------------------------------------------------------------------------
   metrics <- .genomic_eigen_index_metrics(b_PG, Phi, A,
                                            lambda2 = lambda2,
                                            k_I    = selection_intensity)
 
+  # Extract corrected eigenvector (sign-corrected for positive quadratic form)
+  b_PG <- metrics$b
+  
+  # Split coefficients
+  b_y     <- b_PG[1:n_traits]
+  b_gamma <- b_PG[(n_traits + 1):(2 * n_traits)]
+
   E_PG <- metrics$E_vec[1:n_traits]
   names(E_PG) <- trait_names
 

R/helpers-validation.R

@@ -274,5 +274,8 @@ is_symmetric <- function(mat, tolerance = TOL_SYM) {
 #' @keywords internal
 #' @noRd
 is_zero <- function(x, tolerance = TOL_ZERO) {
-  abs(x) < tolerance
+  # Handle NA and Inf values
+  result <- abs(x) < tolerance
+  result[is.na(x) | is.infinite(x)] <- FALSE
+  result
 }

man/gw_esim.Rd

@@ -41,19 +41,43 @@ Eigen::MatrixXd cpp_grouped_sums(
   const int n_obs = data_mat.rows();
   const int n_traits = data_mat.cols();
 
+  // Validate inputs
+  if (n_obs == 0 || n_traits == 0) {
+    return MatrixXd::Zero(0, n_traits);
+  }
+  if (group_idx.size() != n_obs) {
+    Rcpp::stop("group_idx size must match data_mat rows");
+  }
+  
   // Convert 1-based R indices to 0-based C++
   Eigen::VectorXi groups = group_idx.array() - 1;
 
-  // Count unique groups
-  int n_groups = groups.maxCoeff() + 1;
+  // Validate group indices are non-negative
+  if (groups.minCoeff() < 0) {
+    Rcpp::stop("group_idx must contain positive integers");
+  }
 
-  // Pre-allocate result matrix
+  // Find actual number of unique groups safely
+  int min_group = groups.minCoeff();
+  int max_group = groups.maxCoeff();
+  int n_groups = max_group + 1;
+  
+  // Validate dimensions before matrix allocation
+  if (n_groups <= 0 || n_traits <= 0) {
+    Rcpp::stop("Invalid matrix dimensions: n_groups=" + std::to_string(n_groups) + 
+               ", n_traits=" + std::to_string(n_traits));
+  }
+  
+  // Pre-allocate result matrix with validated dimensions
   MatrixXd result(n_groups, n_traits);
   result.setZero();
 
   // Accumulate sums for each group
   for (int i = 0; i < n_obs; ++i) {
-    result.row(groups(i)) += data_mat.row(i);
+    int group = groups(i);
+    if (group >= 0 && group < n_groups) {
+      result.row(group) += data_mat.row(i);
+    }
   }
 
   return result;
@@ -81,17 +105,43 @@ List cpp_multi_grouped_sums(
   const int n_traits = data_mat.cols();
   const int n_groupings = group_indices.size();
 
+  // Validate inputs
+  if (n_obs == 0 || n_traits == 0) {
+    List result(n_groupings);
+    for (int g = 0; g < n_groupings; ++g) {
+      result[g] = MatrixXd::Zero(0, n_traits);
+    }
+    return result;
+  }
+  
   List result(n_groupings);
 
   for (int g = 0; g < n_groupings; ++g) {
     Eigen::VectorXi groups = Rcpp::as<Eigen::Map<Eigen::VectorXi>>(group_indices[g]).array() - 1;
+    
+    // Validate group indices
+    if (groups.size() != n_obs) {
+      Rcpp::stop("group_indices[" + std::to_string(g) + "] size must match data_mat rows");
+    }
+    if (groups.minCoeff() < 0) {
+      Rcpp::stop("group_indices[" + std::to_string(g) + "] must contain positive integers");
+    }
+    
     int n_groups = groups.maxCoeff() + 1;
 
+    // Validate dimensions
+    if (n_groups <= 0) {
+      Rcpp::stop("Invalid number of groups: " + std::to_string(n_groups));
+    }
+    
     MatrixXd group_sums(n_groups, n_traits);
     group_sums.setZero();
 
     for (int i = 0; i < n_obs; ++i) {
-      group_sums.row(groups(i)) += data_mat.row(i);
+      int group = groups(i);
+      if (group >= 0 && group < n_groups) {
+        group_sums.row(group) += data_mat.row(i);
+      }
     }
 
     result[g] = group_sums;