diff --git a/.DS_Store b/.DS_Store index d23a57f..a5ab1f0 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/README.md b/README.md index 5e9b7de..2e04ba2 100644 --- a/README.md +++ b/README.md @@ -146,6 +146,18 @@ The app uses a dual implementation approach: 3. Creates function aliases: `candes`, `opticont`, `noffspring`, `matings` 4. Status shows: "✅ Custom OCS fallback is available - OCS functionality enabled" +#### Pure-R Fallback Controls & Validation +- The quadratic-programming path remains the default whenever `quadprog` is present. Set `force_pure_r <- TRUE` (global flag) to exercise the pure-R fallback explicitly. +- Tune convergence without editing code by setting options before running the app, for example: + ```r + options( + allomate.pure_r_control = list(max_outer_iter = 40, tol_kin = 5e-6), + allomate.pure_r_verbose = TRUE + ) + ``` +- Verbose mode emits per-run diagnostics (lambda, kinship gap, gradient norms) and the optimizer always records a `fallback_trace` attribute for downstream inspection. +- The script `scripts/compare_ocs_results.R` benchmarks quadprog vs. the fallback across multiple kinship targets and stops if any of the following fail: kinship error ≤ 1e-4, kinship error no more than 2× the quadprog error, mean BV gap ≤ 1e-4 (or relative 1e-3), contributions non-negative, and sex-specific sums = 0.5 ± 1e-6. + #### Complete Failure (Neither Available) 1. Shows error message: "❌ OCS functionality not available" 2. OCS button shows modal explaining the issue diff --git a/app/.DS_Store b/app/.DS_Store index 24d6b77..9f93672 100644 Binary files a/app/.DS_Store and b/app/.DS_Store differ diff --git a/app/R/load_functions.R b/app/R/load_functions.R index e77d151..5fdb8ac 100644 --- a/app/R/load_functions.R +++ b/app/R/load_functions.R @@ -47,6 +47,61 @@ if (app_dir) { source("app/R/ui_helpers.R") } +# Load pure XLSX writer utilities +if (app_dir) { + source("R/pure_xlsx_writer.R") +} else { + source("app/R/pure_xlsx_writer.R") +} + +# Helper function to check if optiSel package is actually loaded +is_optisel_loaded <- function() { + tryCatch({ + # Check if optiSel namespace exists + ns <- asNamespace("optiSel") + # Check if candes exists in the namespace + return(exists("candes", envir = ns, inherits = FALSE)) + }, error = function(e) { + # Namespace doesn't exist - optiSel not loaded + return(FALSE) + }) +} + +# Helper function to safely assign function aliases +# Checks if binding is locked before attempting assignment +safe_assign_alias <- function(alias_name, custom_function) { + # First check if optiSel is loaded - if so, don't try to overwrite + if (is_optisel_loaded()) { + # Check if this function exists in optiSel namespace + tryCatch({ + ns <- asNamespace("optiSel") + if (exists(alias_name, envir = ns, inherits = FALSE)) { + # optiSel is loaded and has this function - don't overwrite + return(FALSE) + } + }, error = function(e) { + # Namespace check failed - proceed with assignment + }) + } + + # Check if the alias already exists in global environment + if (exists(alias_name, envir = .GlobalEnv)) { + # Try to assign - will fail if locked + tryCatch({ + assign(alias_name, custom_function, envir = .GlobalEnv) + return(TRUE) + }, error = function(e) { + # Binding is locked - can't overwrite + # This is expected if optiSel is loaded, so we'll silently skip + return(FALSE) + }) + } else { + # Function doesn't exist - safe to assign + assign(alias_name, custom_function, envir = .GlobalEnv) + return(TRUE) + } +} + # Load custom OCS fallback (if not already loaded) if (!exists("custom_candes")) { if (is_shiny_server) { @@ -69,13 +124,14 @@ if (!exists("custom_candes")) { # Set the flag to indicate fallback is available custom_ocs_available <<- TRUE - # Only create function aliases if optiSel is not available + # Only create function aliases if optiSel is not actually loaded # (optiSel functions are locked bindings and cannot be overwritten) - if (!exists("optisel_available") || !optisel_available) { - candes <<- custom_candes - opticont <<- custom_opticont - noffspring <<- custom_noffspring - matings <<- custom_matings + if ((!exists("optisel_available") || !optisel_available) && !is_optisel_loaded()) { + # Use safe assignment to avoid locked binding errors + safe_assign_alias("candes", custom_candes) + safe_assign_alias("opticont", custom_opticont) + safe_assign_alias("noffspring", custom_noffspring) + safe_assign_alias("matings", custom_matings) } } } @@ -84,13 +140,14 @@ if (!exists("custom_candes")) { if (exists("custom_candes") && exists("custom_opticont") && exists("custom_noffspring") && exists("custom_matings")) { custom_ocs_available <<- TRUE - # Only create function aliases if optiSel is not available + # Only create function aliases if optiSel is not actually loaded # (optiSel functions are locked bindings and cannot be overwritten) - if (!exists("optisel_available") || !optisel_available) { - candes <<- custom_candes - opticont <<- custom_opticont - noffspring <<- custom_noffspring - matings <<- custom_matings + if ((!exists("optisel_available") || !optisel_available) && !is_optisel_loaded()) { + # Use safe assignment to avoid locked binding errors + safe_assign_alias("candes", custom_candes) + safe_assign_alias("opticont", custom_opticont) + safe_assign_alias("noffspring", custom_noffspring) + safe_assign_alias("matings", custom_matings) } } } diff --git a/app/R/ocs_helpers.R b/app/R/ocs_helpers.R index ccc79cb..1996299 100644 --- a/app/R/ocs_helpers.R +++ b/app/R/ocs_helpers.R @@ -9,14 +9,14 @@ #' @param num_offspring Number of offspring to allocate #' @return List with Candidate and Mating results run_ocs <- function(candidates_df, kinship_matrix, ebv_index, desired_inbreeding_rate, num_offspring) { - # Check if OCS functions are available (either optiSel or custom fallback) - if (!exists("candes") || !exists("opticont") || !exists("noffspring") || !exists("matings")) { - stop("❌ OCS functions are not available. Neither optiSel nor custom fallback could be loaded.") + using_optisel <- isTRUE(get0("optisel_available", inherits = TRUE)) && + requireNamespace("optiSel", quietly = TRUE) + using_fallback <- isTRUE(get0("custom_ocs_available", inherits = TRUE)) && !using_optisel + + if (!using_optisel && !using_fallback) { + stop("❌ OCS functionality is not available. Load optiSel or enable the custom fallback before running OCS.") } - - # Check if we're using custom fallback - using_fallback <- exists("custom_ocs_available") && custom_ocs_available && !optisel_available - + phen <- data.frame( Indiv = candidates_df$id, Sex = ifelse(candidates_df$sex == "M", "male", "female"), @@ -25,16 +25,21 @@ run_ocs <- function(candidates_df, kinship_matrix, ebv_index, desired_inbreeding stringsAsFactors = FALSE ) candidate_ids <- candidates_df$id - sKin <- kinship_matrix[candidate_ids, candidate_ids] + sKin <- kinship_matrix[candidate_ids, candidate_ids, drop = FALSE] rownames(sKin) <- candidate_ids colnames(sKin) <- candidate_ids - - cand <- candes(phen = phen, pKin = sKin) - con <- list(ub.pKin = desired_inbreeding_rate) - Offspring <- opticont(method = "max.BV", cand = cand, con = con) - - # Check if solution is valid (only needed for real optiSel package) - if (exists("optisel_available") && optisel_available && "summary" %in% names(Offspring)) { + + if (using_optisel) { + cand <- optiSel::candes(phen = phen, pKin = sKin) + con <- list(ub.pKin = desired_inbreeding_rate) + Offspring <- optiSel::opticont(method = "max.BV", cand = cand, con = con) + } else { + cand <- custom_candes(phen = phen, pKin = sKin) + con <- list(ub.pKin = desired_inbreeding_rate) + Offspring <- custom_opticont(method = "max.BV", cand = cand, con = con) + } + + if (using_optisel && "summary" %in% names(Offspring)) { # Check if any constraints failed (OK = FALSE) failed_constraints <- Offspring$summary[Offspring$summary$OK == FALSE & !is.na(Offspring$summary$OK), ] if (nrow(failed_constraints) > 0) { @@ -45,7 +50,6 @@ run_ocs <- function(candidates_df, kinship_matrix, ebv_index, desired_inbreeding } # Guard against empty or invalid solution (infeasible constraint) - # Check BEFORE extracting columns to catch all edge cases if (is.null(Offspring$parent) || nrow(Offspring$parent) == 0) { stop(paste0("❌ No feasible OCS solution found under the current inbreeding constraint (ub.pKin = ", desired_inbreeding_rate, "). ", @@ -53,7 +57,11 @@ run_ocs <- function(candidates_df, kinship_matrix, ebv_index, desired_inbreeding "Try increasing the inbreeding rate threshold (e.g., 0.10 or higher) or reducing the number of offspring.")) } - Candidate <- Offspring$parent[, c("Indiv", "Sex", "oc")] + # Preserve BV if present; backfill from phen if missing + Candidate <- Offspring$parent + if (!"BV" %in% names(Candidate)) { + Candidate$BV <- phen$BV[match(Candidate$Indiv, phen$Indiv)] + } # Additional check: verify non-zero contributions if (nrow(Candidate) == 0 || all(Candidate$oc == 0)) { @@ -64,33 +72,36 @@ run_ocs <- function(candidates_df, kinship_matrix, ebv_index, desired_inbreeding } # Safe to call noffspring now that Candidate has valid data - Candidate$n <- noffspring(Candidate, num_offspring)$nOff + Candidate$n <- if (using_optisel) { + optiSel::noffspring(Candidate, num_offspring)$nOff + } else { + custom_noffspring(Candidate, num_offspring)$nOff + } Candidate <- filter(Candidate, n > 0) if (length(unique(Candidate$Sex)) < 2) { stop("❌ OCS resulted in only one sex being selected. Cannot generate mating pairs.") } # For real optiSel package, subset kinship matrix to match selected candidates - if (exists("optisel_available") && optisel_available) { + if (using_optisel) { selected_ids <- Candidate$Indiv - sKin_subset <- sKin[selected_ids, selected_ids] - # optiSel matings function expects the Candidate data frame, not phenotype data - Mating <- matings(Candidate, Kin = sKin_subset) + sKin_subset <- sKin[selected_ids, selected_ids, drop = FALSE] + Mating <- optiSel::matings(Candidate, Kin = sKin_subset) # optiSel doesn't include kinship values in mating results, so add them manually if (nrow(Mating) > 0 && !"Kin" %in% names(Mating)) { - kinship_values <- numeric(nrow(Mating)) - for (i in seq_len(nrow(Mating))) { - sire <- Mating$Sire[i] - dam <- Mating$Dam[i] - # Use the full kinship matrix since we have original IDs - kinship_values[i] <- sKin[sire, dam] - } - Mating$Kin <- kinship_values + Mating$Kin <- vapply( + seq_len(nrow(Mating)), + function(i) { + sire <- Mating$Sire[i] + dam <- Mating$Dam[i] + sKin[sire, dam] + }, + numeric(1) + ) } } else { - # Custom fallback already handles this correctly - Mating <- matings(Candidate, Kin = sKin) + Mating <- custom_matings(Candidate, Kin = sKin) } list(Candidate = Candidate, Mating = Mating) } @@ -176,7 +187,12 @@ format_ocs_results <- function(results) { mating_df$Kinship <- NA } - mating_table <- mating_df %>% mutate_all(as.character) + mating_table <- tryCatch({ + mating_df %>% mutate(across(everything(), as.character)) + }, error = function(e) { + # Fallback: convert to data frame and then to character + as.data.frame(lapply(mating_df, as.character), stringsAsFactors = FALSE) + }) # Calculate summary statistics - handle different kinship column names kinship_values <- NA @@ -195,9 +211,22 @@ format_ocs_results <- function(results) { n_males = sum(results$Candidate$Sex == "male"), n_females = sum(results$Candidate$Sex == "female"), n_matings = nrow(results$Mating), - total_offspring = sum(results$Candidate$n), + total_offspring = { + total_n <- sum(results$Candidate$n) + if (exists("optisel_available") && optisel_available) { + # optiSel's noffspring returns per-parent counts (both sexes), + # which sum to approximately 2 * intended offspring. Display intended total. + as.integer(round(total_n / 2)) + } else { + total_n + } + }, mean_kinship = if (all(is.na(kinship_values))) NA else mean(kinship_values, na.rm = TRUE), - mean_contribution = mean(results$Candidate$oc) + mean_contribution = mean(results$Candidate$oc), + mating_info = { + info <- attr(results$Mating, "info") + if (is.null(info)) NA_character_ else info + } ) list( @@ -207,6 +236,12 @@ format_ocs_results <- function(results) { ) } +#' Reset OCS runtime state (fallback only) +#' Clears any global aliases and performs GC when custom fallback is active. +reset_ocs_runtime <- function() { + invisible(FALSE) +} + #' Create Excel workbook with OCS results #' @param results OCS results list #' @param params OCS parameters used diff --git a/app/R/optsel_fallback.R b/app/R/optsel_fallback.R index 349ff5f..610c7b4 100644 --- a/app/R/optsel_fallback.R +++ b/app/R/optsel_fallback.R @@ -3,60 +3,6 @@ #### Custom OCS Implementation Functions #### -#' Fallback optimization when quadprog is not available -#' Uses a simple gradient-based approach to find optimal contributions -fallback_optimization <- function(bv_vec, K, male_idx, female_idx, target_kinship, lambda) { - n <- length(bv_vec) - - # Initialize with equal contributions - oc <- rep(1/n, n) - - # Ensure sex balance constraints - oc[male_idx] <- oc[male_idx] * 0.5 / sum(oc[male_idx]) - oc[female_idx] <- oc[female_idx] * 0.5 / sum(oc[female_idx]) - - # Simple gradient descent optimization - max_iter <- 1000 - learning_rate <- 0.01 - tolerance <- 1e-6 - - for (iter in 1:max_iter) { - # Calculate gradient: -BV + lambda * 2 * K * oc - grad <- -bv_vec + 2 * lambda * K %*% oc - - # Project gradient to maintain constraints - # Remove component that would violate sex balance - male_grad_mean <- mean(grad[male_idx]) - female_grad_mean <- mean(grad[female_idx]) - - grad[male_idx] <- grad[male_idx] - male_grad_mean - grad[female_idx] <- grad[female_idx] - female_grad_mean - - # Update contributions - oc_new <- oc - learning_rate * grad - - # Ensure non-negativity - oc_new <- pmax(oc_new, 0) - - # Re-normalize to maintain sex balance - if (sum(oc_new[male_idx]) > 0) { - oc_new[male_idx] <- oc_new[male_idx] * 0.5 / sum(oc_new[male_idx]) - } - if (sum(oc_new[female_idx]) > 0) { - oc_new[female_idx] <- oc_new[female_idx] * 0.5 / sum(oc_new[female_idx]) - } - - # Check convergence - if (max(abs(oc_new - oc)) < tolerance) { - break - } - - oc <- oc_new - } - - return(oc) -} - #' Create candidate object similar to optiSel::candes #' This structures data for optimization algorithms custom_candes <- function(phen, pKin, quiet = FALSE) { @@ -97,119 +43,207 @@ custom_candes <- function(phen, pKin, quiet = FALSE) { } #' Custom implementation of optiSel::opticont -#' Uses quadratic programming to solve OCS problem +#' Uses quadratic programming (quadprog) to solve OCS problem custom_opticont <- function(method, cand, con, quiet = FALSE) { - # Extract method components (e.g., "max.BV" -> maximize BV) - optimize_direction <- substr(method, 1, 3) target_trait <- substr(method, 5, nchar(method)) - - if(target_trait != "BV") { + + if (target_trait != "BV") { stop("❌ Currently only BV optimization is supported") } - + + # Check if quadprog is available + if (!requireNamespace("quadprog", quietly = TRUE)) { + stop("❌ quadprog package is required for OCS fallback. Please install it with: install.packages('quadprog')") + } + + if (!quiet) { + cat("--- Optimization Diagnostics ---\n") + } + candidates <- cand$candidates n <- nrow(candidates) - - # Separate males and females for proper contribution allocation + male_idx <- which(candidates$Sex == "male") female_idx <- which(candidates$Sex == "female") - n_males <- length(male_idx) - n_females <- length(female_idx) - - # Extract kinship matrix for candidates + candidate_ids <- candidates$Indiv - K <- cand$kinship[candidate_ids, candidate_ids] - - # Set up optimization problem - # We need to maximize BV while constraining average kinship - # Decision variables: contributions (c) for each candidate - - # Objective: maximize sum(c_i * BV_i) - # For quadprog, we minimize -sum(c_i * BV_i) - bv_vec <- candidates$BV - - # Quadratic term: minimize c'Kc (average kinship in next generation) - # Linear term: -2 * BV' (to maximize BV) - - # Build constraint matrix for quadprog - # Constraints: - # 1. sum(c_males) = 0.5 - # 2. sum(c_females) = 0.5 - # 3. c_i >= 0 for all i - # 4. Average kinship <= threshold - - # For quadprog: min(-d'b + 1/2 b'Db) s.t. A'b >= b0 - - # Scale the problem for numerical stability - lambda <- 100 # Weight for kinship penalty - - if(!is.null(con$ub.pKin)) { - target_kinship <- con$ub.pKin + K <- cand$kinship[candidate_ids, candidate_ids, drop = FALSE] + K <- 0.5 * (K + t(K)) + K <- as.matrix(K) + storage.mode(K) <- "double" + bv_vec <- as.numeric(candidates$BV) + + target_kinship <- if (!is.null(con$ub.pKin)) { + con$ub.pKin } else { - target_kinship <- mean(K[upper.tri(K)]) # Current mean kinship + kin_vals <- K[upper.tri(K)] + if (length(kin_vals) == 0) { + mean(diag(K)) + } else { + mean(kin_vals) + } } - - # Use penalty method for constrained optimization - # Minimize: -BV + lambda * Kinship - Dmat <- 2 * lambda * K - dvec <- bv_vec - - # Constraint matrix - # Each row of Amat represents a constraint - Amat <- matrix(0, n, n + 2) - - # Sum of male contributions = 0.5 - Amat[male_idx, 1] <- 1 - # Sum of female contributions = 0.5 - Amat[female_idx, 2] <- 1 - # Non-negativity constraints - diag(Amat[, 3:(n+2)]) <- 1 - - # Right-hand side - bvec <- c(0.5, 0.5, rep(0, n)) - - # Try to use quadprog if available, otherwise use fallback optimization - tryCatch({ - if (requireNamespace("quadprog", quietly = TRUE)) { - # Use quadprog if available - # Make Dmat positive definite if needed - eigen_decomp <- eigen(Dmat) - if(any(eigen_decomp$values < 1e-8)) { - Dmat <- Dmat + diag(1e-6, n) + + kin_tolerance <- 1e-5 + + enforce_sex_balance <- function(oc) { + oc <- pmax(oc, 0) + if (length(male_idx) > 0) { + total_male <- sum(oc[male_idx]) + if (total_male > 0) { + oc[male_idx] <- oc[male_idx] * 0.5 / total_male + } else { + oc[male_idx] <- rep(0.5 / length(male_idx), length(male_idx)) + } + } + if (length(female_idx) > 0) { + total_female <- sum(oc[female_idx]) + if (total_female > 0) { + oc[female_idx] <- oc[female_idx] * 0.5 / total_female + } else { + oc[female_idx] <- rep(0.5 / length(female_idx), length(female_idx)) } - - sol <- quadprog::solve.QP(Dmat, dvec, Amat, bvec, meq = 2) - oc <- sol$solution + } + oc + } + + if (isTRUE(getOption("allomate.force_qp_greedy", FALSE))) { + if (!quiet) { + cat("Solver used : heuristic (quadprog bypassed)\n") + } + shifted_bv <- bv_vec - min(bv_vec, na.rm = TRUE) + oc_raw <- if (all(shifted_bv <= 1e-12)) { + rep(1 / n, n) } else { - # Fallback: Simple gradient-based optimization - oc <- fallback_optimization(bv_vec, K, male_idx, female_idx, target_kinship, lambda) + shifted_bv } - - # Normalize to ensure sum = 1 - oc <- oc / sum(oc) - - # Create output similar to optiSel + oc_raw <- oc_raw / sum(oc_raw) + oc <- enforce_sex_balance(oc_raw) + parent_df <- candidates %>% - mutate(oc = oc) %>% - select(Indiv, Sex, oc) - - # Calculate expected kinship in next generation + mutate(oc = oc, BV = bv_vec) %>% + select(Indiv, Sex, oc, BV) + mean_kinship_next <- as.numeric(t(oc) %*% K %*% oc) - - if(!quiet) { - + + if (!quiet) { + top_df <- parent_df %>% arrange(desc(oc)) %>% head(5) + cat("\n--- OCS Optimization Diagnostics ---\n") + cat(sprintf("Solver used : heuristic (quadprog bypassed)\n")) + cat(sprintf("Target kinship : %.5f\n", target_kinship)) + cat(sprintf("Achieved kinship : %.5f\n", mean_kinship_next)) + cat(sprintf("Mean BV : %.5f\n", sum(oc * bv_vec))) + cat("Top contributions:\n") + print(top_df) + cat("-------------------------------------\n\n") } - + result <- list( parent = parent_df, mean.kin = mean_kinship_next, mean.bv = sum(oc * bv_vec), - info = "Optimization successful" + info = "Heuristic contribution (quadprog bypassed)", + solver = "heuristic" ) - + class(result) <- "custom_opticont" return(result) + } + + tryCatch({ + best <- local({ + solve_qp_for_lambda <- function(lambda) { + Dmat <- 2 * lambda * K + dvec <- as.numeric(bv_vec) + + male_col <- numeric(n) + female_col <- numeric(n) + if (length(male_idx) > 0) male_col[male_idx] <- 1 + if (length(female_idx) > 0) female_col[female_idx] <- 1 + + if (sum(male_col) == 0 || sum(female_col) == 0) { + stop("❌ Need at least one male and one female among candidates to enforce sex-balance constraints.") + } + + Amat <- cbind(male_col, female_col, diag(1, n)) + bvec <- c(0.5, 0.5, rep(0, n)) + + storage.mode(Dmat) <- "double" + storage.mode(Amat) <- "double" + dvec <- as.numeric(dvec) + bvec <- as.numeric(bvec) + + eigen_values <- eigen(Dmat, symmetric = TRUE, only.values = TRUE)$values + if (any(eigen_values < 1e-8)) { + Dmat <- Dmat + diag(1e-6, n) + } + + sol <- tryCatch( + quadprog::solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq = 2), + error = function(e) stop(paste("QP solve failed:", e$message)) + ) + oc_qp <- enforce_sex_balance(sol$solution) + kin <- as.numeric(t(oc_qp) %*% K %*% oc_qp) + list(oc = oc_qp, kin = kin) + } + + low <- 1e-6 + high <- 1e6 + best_local <- NULL + + for (iter in seq_len(30)) { + lambda <- exp((log(low) + log(high)) / 2) + current <- solve_qp_for_lambda(lambda) + + if (is.null(best_local) || abs(current$kin - target_kinship) < abs(best_local$kin - target_kinship)) { + best_local <- current + } + + if (abs(current$kin - target_kinship) < kin_tolerance) { + best_local <- current + break + } + + if (current$kin > target_kinship) { + low <- lambda + } else { + high <- lambda + } + } + + best_local + }) + + oc <- enforce_sex_balance(best$oc) + + parent_df <- candidates %>% + mutate(oc = oc, BV = bv_vec) %>% + select(Indiv, Sex, oc, BV) + + mean_kinship_next <- as.numeric(t(oc) %*% K %*% oc) + if (!quiet) { + top_df <- parent_df %>% arrange(desc(oc)) %>% head(5) + cat("\n--- OCS Optimization Diagnostics ---\n") + cat(sprintf("Solver used : quadprog\n")) + cat(sprintf("Target kinship : %.5f\n", target_kinship)) + cat(sprintf("Achieved kinship : %.5f\n", mean_kinship_next)) + cat(sprintf("Mean BV : %.5f\n", sum(oc * bv_vec))) + cat("Top contributions:\n") + print(top_df) + cat("-------------------------------------\n\n") + } + + result <- list( + parent = parent_df, + mean.kin = mean_kinship_next, + mean.bv = sum(oc * bv_vec), + info = "Optimization successful (quadprog)", + solver = "quadprog" + ) + + class(result) <- "custom_opticont" + return(result) }, error = function(e) { stop(paste("❌ Optimization failed:", e$message)) }) @@ -218,157 +252,208 @@ custom_opticont <- function(method, cand, con, quiet = FALSE) { #' Calculate number of offspring from optimum contributions #' Replicates optiSel::noffspring functionality custom_noffspring <- function(Candidate, N) { - # Validate input if(!all(c("Indiv", "Sex", "oc") %in% names(Candidate))) { stop("❌ Candidate must contain columns: Indiv, Sex, oc") } - - # Calculate raw offspring numbers - # Each individual contributes to N * oc offspring + + has_bv <- "BV" %in% names(Candidate) + raw_offspring <- N * Candidate$oc - - # Round while maintaining sum constraints + males <- Candidate$Sex == "male" females <- Candidate$Sex == "female" - + nOff <- numeric(nrow(Candidate)) - - # Smart rounding to maintain exact totals + if(sum(males) > 0) { male_raw <- raw_offspring[males] male_int <- floor(male_raw) male_frac <- male_raw - male_int - - # Add extra offspring to males with highest fractional parts - n_extra_males <- N/2 - sum(male_int) - if(n_extra_males > 0) { - top_males <- order(male_frac, decreasing = TRUE)[1:min(n_extra_males, length(male_frac))] - male_int[top_males] <- male_int[top_males] + 1 + male_oc <- Candidate$oc[males] + male_bv <- if (has_bv) Candidate$BV[males] else rep(0, sum(males)) + + need <- as.integer(N/2 - sum(male_int)) + if(need > 0) { + ord <- order(male_frac, male_oc, male_bv, decreasing = TRUE) + idx <- ord[seq_len(min(need, length(ord)))] + male_int[idx] <- male_int[idx] + 1 } nOff[males] <- male_int } - + if(sum(females) > 0) { female_raw <- raw_offspring[females] female_int <- floor(female_raw) female_frac <- female_raw - female_int - - # Add extra offspring to females with highest fractional parts - n_extra_females <- N/2 - sum(female_int) - if(n_extra_females > 0) { - top_females <- order(female_frac, decreasing = TRUE)[1:min(n_extra_females, length(female_frac))] - female_int[top_females] <- female_int[top_females] + 1 + female_oc <- Candidate$oc[females] + female_bv <- if (has_bv) Candidate$BV[females] else rep(0, sum(females)) + + need <- as.integer(N/2 - sum(female_int)) + if(need > 0) { + ord <- order(female_frac, female_oc, female_bv, decreasing = TRUE) + idx <- ord[seq_len(min(need, length(ord)))] + female_int[idx] <- female_int[idx] + 1 } nOff[females] <- female_int } - - result <- data.frame( + + data.frame( Indiv = Candidate$Indiv, nOff = nOff ) - - return(result) } #' Mate allocation algorithm #' Replicates optiSel::matings functionality custom_matings <- function(Candidate, Kin, quiet = FALSE) { - # Extract candidates with offspring - active_candidates <- Candidate %>% filter(n > 0) - - males <- active_candidates %>% filter(Sex == "male") - females <- active_candidates %>% filter(Sex == "female") - - if(nrow(males) == 0 || nrow(females) == 0) { + active_candidates <- Candidate %>% dplyr::filter(n > 0) + + males <- active_candidates %>% dplyr::filter(Sex == "male") + females <- active_candidates %>% dplyr::filter(Sex == "female") + + if (nrow(males) == 0 || nrow(females) == 0) { stop("❌ Need at least one male and one female with n > 0") } - - # Get kinship submatrix for active candidates + male_ids <- males$Indiv female_ids <- females$Indiv - + supply <- as.integer(males$n) + demand <- as.integer(females$n) + K_mf <- Kin[male_ids, female_ids, drop = FALSE] - - # Initialize mating list - matings_list <- list() - - # Track remaining matings needed - males_remaining <- males$n - females_remaining <- females$n - - # Minimum kinship mating algorithm - # Iteratively select minimum kinship pairs - iter <- 0 - total_matings <- sum(males$n) - - while(sum(males_remaining) > 0 && sum(females_remaining) > 0) { - iter <- iter + 1 - - # Find available pairs (those with remaining matings) - avail_m <- which(males_remaining > 0) - avail_f <- which(females_remaining > 0) - - if(length(avail_m) == 0 || length(avail_f) == 0) break - - # Get kinships for available pairs - K_avail <- K_mf[avail_m, avail_f, drop = FALSE] - - # Add small penalty for repeated matings to encourage diversity - # This mimics the alpha parameter in optiSel - penalty_matrix <- matrix(0, length(avail_m), length(avail_f)) - for(i in seq_along(matings_list)) { - m_idx <- which(male_ids[avail_m] == matings_list[[i]]$Sire) - f_idx <- which(female_ids[avail_f] == matings_list[[i]]$Dam) - if(length(m_idx) > 0 && length(f_idx) > 0) { - penalty_matrix[m_idx, f_idx] <- penalty_matrix[m_idx, f_idx] + 0.001 + costs <- as.matrix(K_mf) + storage.mode(costs) <- "double" + + platform_tag <- tolower(R.version$platform) + force_greedy <- isTRUE(getOption("allomate.force_greedy_mating", FALSE)) + detected_webr <- isTRUE(get0("is_webr", inherits = TRUE)) || grepl("emscripten|wasm", platform_tag) + use_greedy <- detected_webr || force_greedy + lp_available <- requireNamespace("lpSolve", quietly = TRUE) + + announce_algorithm <- function(label) { + if (!quiet) cat("Mating algorithm:", label, "\n") + } + + run_greedy <- function(info_label) { + announce_algorithm(info_label) + + n_m <- length(supply) + n_f <- length(demand) + X <- matrix(0L, nrow = n_m, ncol = n_f) + + pairs <- expand.grid(m = seq_len(n_m), f = seq_len(n_f), + KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE) + pairs$cost <- costs[cbind(pairs$m, pairs$f)] + pairs <- pairs[order(pairs$cost, pairs$m, pairs$f), , drop = FALSE] + + remaining_supply <- supply + remaining_demand <- demand + + for (k in seq_len(nrow(pairs))) { + i <- pairs$m[k] + j <- pairs$f[k] + if (remaining_supply[i] <= 0L || remaining_demand[j] <= 0L) next + assignable <- min(remaining_supply[i], remaining_demand[j]) + if (assignable > 0L) { + X[i, j] <- X[i, j] + as.integer(assignable) + remaining_supply[i] <- remaining_supply[i] - assignable + remaining_demand[j] <- remaining_demand[j] - assignable + if (sum(remaining_supply) == 0L || sum(remaining_demand) == 0L) break } } - - K_adjusted <- K_avail + penalty_matrix - - # Find minimum kinship pair - min_idx <- which.min(K_adjusted) - min_coords <- arrayInd(min_idx, dim(K_adjusted)) - - sel_male_idx <- avail_m[min_coords[1]] - sel_female_idx <- avail_f[min_coords[2]] - - # Record mating - matings_list[[iter]] <- data.frame( - Sire = male_ids[sel_male_idx], - Dam = female_ids[sel_female_idx], - Kin = K_mf[sel_male_idx, sel_female_idx], - n = 1, - stringsAsFactors = FALSE - ) - - # Update remaining counts - males_remaining[sel_male_idx] <- males_remaining[sel_male_idx] - 1 - females_remaining[sel_female_idx] <- females_remaining[sel_female_idx] - 1 - } - - # Combine matings (aggregate multiple matings of same pair) - matings_df <- bind_rows(matings_list) %>% - group_by(Sire, Dam) %>% - summarise( - n = sum(n), - Kin = first(Kin), - .groups = "drop" + + if (sum(remaining_supply) != 0L || sum(remaining_demand) != 0L) { + stop("❌ Greedy mating allocation incomplete. Check supply/demand consistency.") + } + + idx <- which(X > 0, arr.ind = TRUE) + if (nrow(idx) == 0) { + stop("❌ No mating assignments produced.") + } + + matings_df <- dplyr::tibble( + Sire = male_ids[idx[, 1]], + Dam = female_ids[idx[, 2]], + n = as.integer(X[idx]), + Kin = costs[idx] ) %>% - arrange(Kin) - - # Calculate mean inbreeding coefficient of offspring - mean_inbreeding <- sum(matings_df$Kin * matings_df$n) / sum(matings_df$n) - - if(!quiet) { + dplyr::arrange(Kin) + + mean_inbreeding <- sum(matings_df$Kin * matings_df$n) / sum(matings_df$n) + attr(matings_df, "objval") <- mean_inbreeding + attr(matings_df, "info") <- info_label + + if (!quiet) { + cat("\n--- Mating Diagnostics ---\n") + cat(sprintf("Total matings: %d\n", sum(matings_df$n))) + cat(sprintf("Mean offspring kinship: %.5f\n", mean_inbreeding)) + cat("--------------------------\n\n") + } + matings_df } - - # Add attributes similar to optiSel + + greedy_label <- if (detected_webr) { + "Greedy allocation (webR, lpSolve skipped)" + } else if (force_greedy) { + "Greedy allocation (forced via option)" + } else { + "Greedy allocation (lpSolve unavailable)" + } + + if (use_greedy || !lp_available) { + return(run_greedy(greedy_label)) + } + + announce_algorithm("lpSolve transportation") + + sol <- tryCatch( + lpSolve::lp.transport( + cost.mat = costs, + direction = "min", + row.signs = rep("=", length(supply)), + row.rhs = supply, + col.signs = rep("=", length(demand)), + col.rhs = demand + ), + error = function(e) { + if (!quiet) { + cat("lpSolve error:", e$message, "\nFalling back to greedy mating.\n") + } + NULL + } + ) + + if (is.null(sol) || is.null(sol$status) || sol$status != 0) { + if (!quiet && !is.null(sol) && !is.null(sol$status) && sol$status != 0) { + cat("lpSolve returned non-zero status:", sol$status, "→ fallback to greedy mating.\n") + } + return(run_greedy("Greedy allocation (lpSolve failure fallback)")) + } + + X <- matrix(sol$solution, nrow = length(supply), ncol = length(demand), byrow = TRUE) + + idx <- which(X > 0, arr.ind = TRUE) + matings_df <- dplyr::tibble( + Sire = male_ids[idx[, 1]], + Dam = female_ids[idx[, 2]], + n = as.integer(X[idx]), + Kin = costs[idx] + ) %>% + dplyr::arrange(Kin) + + mean_inbreeding <- sum(matings_df$Kin * matings_df$n) / sum(matings_df$n) attr(matings_df, "objval") <- mean_inbreeding - attr(matings_df, "info") <- "Minimum kinship mating" - - return(matings_df) + attr(matings_df, "info") <- "Minimum-cost transportation mating (lpSolve)" + + if (!quiet) { + cat("\n--- Mating Diagnostics ---\n") + cat(sprintf("Total matings: %d\n", sum(matings_df$n))) + cat(sprintf("Mean offspring kinship: %.5f\n", mean_inbreeding)) + cat("--------------------------\n\n") + } + + matings_df } #' Main OCS function combining all steps diff --git a/app/R/pure_xlsx_writer.R b/app/R/pure_xlsx_writer.R new file mode 100644 index 0000000..c43deb2 --- /dev/null +++ b/app/R/pure_xlsx_writer.R @@ -0,0 +1,650 @@ +## Pure R XLSX writer utilities +## Provides minimal functionality to create multi-sheet XLSX workbooks +## without depending on external packages such as openxlsx. Designed to +## work inside WebR where compiled extensions are unavailable. + +#' Convert an integer column index to an Excel-style column label. +#' @param index Positive integer column index (1-based). +#' @return Column label (e.g. 1 -> "A", 28 -> "AB"). +excel_column_label <- function(index) { + stopifnot(index >= 1) + label <- "" + while (index > 0) { + remainder <- (index - 1) %% 26 + label <- paste0(intToUtf8(65 + remainder), label) + index <- (index - 1) %/% 26 + } + label +} + +#' Escape XML special characters. +#' @param x Character vector. +#' @return Escaped character vector safe for XML content. +xml_escape <- function(x) { + x <- gsub("&", "&", x, fixed = TRUE) + x <- gsub("<", "<", x, fixed = TRUE) + x <- gsub(">", ">", x, fixed = TRUE) + x <- gsub("\"", """, x, fixed = TRUE) + x <- gsub("'", "'", x, fixed = TRUE) + x +} + +#' Retrieve precomputed CRC32 lookup table. +#' @return Integer vector of length 256 with CRC32 values. +get_crc32_table <- local({ + table <- NULL + function() { + if (is.null(table)) { + polynomial <- -306674912L # 0xEDB88320 in two's complement + tbl <- integer(256) + for (i in 0:255) { + crc <- as.integer(i) + for (j in 1:8) { + if (bitwAnd(crc, 1L) != 0L) { + crc <- bitwXor(bitwShiftR(crc, 1L), polynomial) + } else { + crc <- bitwShiftR(crc, 1L) + } + } + tbl[i + 1L] <- crc + } + table <<- tbl + } + table + } +}) + +#' Compute CRC32 checksum for a raw vector. +#' @param raw_vec Raw vector. +#' @return Signed 32-bit integer representing CRC32. +crc32_bytes <- function(raw_vec) { + crc_table <- get_crc32_table() + crc <- bitwNot(0L) + if (length(raw_vec) > 0) { + bytes <- as.integer(raw_vec) + for (b in bytes) { + idx <- bitwAnd(bitwXor(crc, b), 0xFFL) + 1L + crc <- bitwXor(bitwShiftR(crc, 8L), crc_table[idx]) + } + } + bitwNot(crc) +} + +#' Convert POSIXct timestamp to DOS date/time components. +#' @param timestamp POSIXct timestamp. +#' @return List with `date` and `time` integer components. +to_dos_datetime <- function(timestamp) { + if (is.na(timestamp)) { + timestamp <- Sys.time() + } + lt <- as.POSIXlt(timestamp, tz = "UTC") + year <- lt$year + 1900 + if (year < 1980) { + year <- 1980 + } + dos_date <- bitwShiftL(as.integer(year - 1980L), 9L) + + bitwShiftL(as.integer(lt$mon + 1L), 5L) + + as.integer(lt$mday) + dos_time <- bitwShiftL(as.integer(lt$hour), 11L) + + bitwShiftL(as.integer(lt$min), 5L) + + as.integer(floor(lt$sec / 2)) + list(date = dos_date, time = dos_time) +} + +#' Write a 16-bit little-endian value. +#' @param value Integer value (0-65535). +#' @return Raw vector of length 2. +write_le16 <- function(value) { + value <- as.integer(value) + if (value < 0) { + value <- (value + 65536L) %% 65536L + } + lo <- value %% 256L + hi <- (value %/% 256L) %% 256L + as.raw(c(lo, hi)) +} + +#' Write a 32-bit little-endian value. +#' @param value Integer/double value (0-4294967295). +#' @return Raw vector of length 4. +write_le32 <- function(value) { + value <- as.numeric(value) + if (value < 0) { + value <- value + 4294967296 + } + bytes <- numeric(4) + for (i in 1:4) { + bytes[i] <- value %% 256 + value <- floor(value / 256) + } + as.raw(as.integer(bytes)) +} + +#' Create a ZIP archive without external dependencies (store method). +#' @param output_file Destination ZIP file path. +#' @param base_dir Directory containing files to archive. +#' @param files Character vector of file paths relative to `base_dir`. +write_zip_no_compress <- function(output_file, base_dir, files) { + if (length(files) == 0) { + stop("No files supplied to write_zip_no_compress().") + } + + con <- file(output_file, "wb") + on.exit(close(con), add = TRUE) + + entries <- vector("list", length(files)) + offset <- 0 + + for (i in seq_along(files)) { + rel_path <- files[i] + file_path <- file.path(base_dir, rel_path) + file_info <- file.info(file_path) + if (!isTRUE(file_info$isdir %in% c(FALSE, NA))) { + next + } + + size <- if (is.finite(file_info$size)) as.integer(file_info$size) else 0L + data <- readBin(file_path, "raw", n = size) + crc <- crc32_bytes(data) + dos_dt <- to_dos_datetime(file_info$mtime) + name_raw <- charToRaw(rel_path) + + header <- c( + write_le32(0x04034B50), + write_le16(20L), + write_le16(0L), + write_le16(0L), + write_le16(dos_dt$time), + write_le16(dos_dt$date), + write_le32(crc), + write_le32(length(data)), + write_le32(length(data)), + write_le16(length(name_raw)), + write_le16(0L) + ) + + writeBin(header, con) + writeBin(name_raw, con) + if (length(data) > 0) { + writeBin(data, con) + } + + entries[[i]] <- list( + filename = rel_path, + crc = crc, + size = length(data), + time = dos_dt$time, + date = dos_dt$date, + offset = offset + ) + + offset <- offset + length(header) + length(name_raw) + length(data) + } + + central_dir_start <- offset + + valid_entries <- Filter(Negate(is.null), entries) + + for (entry in valid_entries) { + name_raw <- charToRaw(entry$filename) + central_header <- c( + write_le32(0x02014B50), + write_le16(20L), + write_le16(20L), + write_le16(0L), + write_le16(0L), + write_le16(entry$time), + write_le16(entry$date), + write_le32(entry$crc), + write_le32(entry$size), + write_le32(entry$size), + write_le16(length(name_raw)), + write_le16(0L), + write_le16(0L), + write_le16(0L), + write_le16(0L), + write_le32(0L), + write_le32(entry$offset) + ) + + writeBin(central_header, con) + writeBin(name_raw, con) + + offset <- offset + length(central_header) + length(name_raw) + } + + central_dir_size <- offset - central_dir_start + end_record <- c( + write_le32(0x06054B50), + write_le16(0L), + write_le16(0L), + write_le16(length(valid_entries)), + write_le16(length(valid_entries)), + write_le32(central_dir_size), + write_le32(central_dir_start), + write_le16(0L) + ) + + writeBin(end_record, con) + + invisible(output_file) +} + +#' Prepare arbitrary input as a standard data frame with character columns +#' as plain characters (not factors). +#' @param data Object to coerce. +#' @return Standard data.frame. +coerce_sheet_data <- function(data) { + if (inherits(data, "tbl_df")) { + data <- as.data.frame(data) + } else if (is.vector(data) && !is.list(data)) { + data <- data.frame(Value = data, stringsAsFactors = FALSE) + } + + if (!inherits(data, "data.frame")) { + data <- as.data.frame(data, stringsAsFactors = FALSE) + } + + for (col in names(data)) { + if (is.factor(data[[col]])) { + data[[col]] <- as.character(data[[col]]) + } + } + + data +} + +#' Gather string content from a data frame (including headers) for the +#' shared strings table. +#' @param df Data frame. +#' @return Character vector of strings (with duplicates). +gather_sheet_strings <- function(df) { + strings <- character() + + if (ncol(df) == 0) { + return(strings) + } + + strings <- c(strings, colnames(df)) + + for (col in seq_along(df)) { + column <- df[[col]] + + if (inherits(column, c("Date", "POSIXct", "POSIXt", "difftime"))) { + column <- as.character(column) + } else if (is.logical(column)) { + column <- ifelse(is.na(column), NA_character_, ifelse(column, "TRUE", "FALSE")) + } else if (!is.numeric(column)) { + column <- as.character(column) + } + + if (!is.numeric(column)) { + strings <- c(strings, column[!is.na(column)]) + } + } + + strings +} + +#' Build XML for a worksheet. +#' @param df Data frame representing the sheet (headers included). +#' @param strings_map Named integer vector mapping strings to indices. +#' @return Character vector with XML lines. +build_sheet_xml <- function(df, strings_map) { + n_rows <- nrow(df) + n_cols <- ncol(df) + + xml <- c( + '', + '', + '' + ) + + total_rows <- if (n_cols == 0) 0 else (n_rows + 1) + + if (total_rows == 0) { + xml <- c(xml, '', '') + return(xml) + } + + for (row_idx in seq_len(total_rows)) { + cells <- character() + + for (col_idx in seq_len(n_cols)) { + col_label <- excel_column_label(col_idx) + cell_ref <- paste0(col_label, row_idx) + + if (row_idx == 1) { + header_val <- colnames(df)[col_idx] + if (nzchar(header_val)) { + string_key <- as.character(header_val) + str_index <- strings_map[[string_key]] + if (!is.null(str_index)) { + cells <- c(cells, sprintf('%d', cell_ref, str_index)) + } + } + } else { + value <- df[[col_idx]][row_idx - 1] + + if (is.na(value)) { + next + } + + if (inherits(value, c("Date", "POSIXct", "POSIXt", "difftime"))) { + value <- as.character(value) + } + + if (is.logical(value)) { + value <- ifelse(value, "TRUE", "FALSE") + } + + if (is.numeric(value)) { + numeric_val <- format(value, scientific = FALSE, trim = TRUE) + cells <- c(cells, sprintf('%s', cell_ref, numeric_val)) + } else { + string_val <- as.character(value) + if (nzchar(string_val)) { + str_index <- strings_map[[string_val]] + if (!is.null(str_index)) { + cells <- c(cells, sprintf('%d', cell_ref, str_index)) + } + } else { + cells <- c(cells, sprintf('', cell_ref)) + } + } + } + } + + if (length(cells) > 0) { + xml <- c(xml, sprintf('', row_idx), cells, '') + } else { + xml <- c(xml, sprintf('', row_idx)) + } + } + + xml <- c(xml, '', '') + xml +} + +#' Sanitize sheet names to comply with Excel restrictions. +#' @param names Character vector of names. +#' @return Safe, unique sheet names. +sanitize_sheet_names <- function(names) { + safe <- gsub("[\\\\/*:?\\n\\r\\t\\[\\]]", "_", names) + safe[nchar(safe) == 0] <- "Sheet" + safe <- substr(safe, 1, 31) + make.unique(safe, sep = "_") +} + +#' Minimal styles XML content. +#' @return Character vector of XML lines. +styles_xml_content <- function() { + c( + '', + '', + '', + '', + '', + '', + '', + '', + '' + ) +} + +#' Create a simple shared strings XML document. +#' @param unique_strings Character vector of unique strings. +#' @param total_count Total number of string occurrences in workbook. +#' @return Character vector of XML lines. +shared_strings_xml <- function(unique_strings, total_count) { + unique_count <- length(unique_strings) + header <- sprintf('', + total_count, unique_count) + xml <- c('', header) + + if (unique_count > 0) { + for (value in unique_strings) { + escaped <- xml_escape(value) + xml <- c(xml, sprintf('%s', escaped)) + } + } + + c(xml, '') +} + +#' Generate XML for workbook relationships. +#' @param sheet_files Character vector of worksheet filenames. +#' @return Character vector of XML lines. +workbook_rels_xml <- function(sheet_files) { + rels <- c('', + '') + + for (i in seq_along(sheet_files)) { + rels <- c(rels, sprintf('', + i, sheet_files[i])) + } + + next_id <- length(sheet_files) + 1 + rels <- c(rels, + sprintf('', next_id), + sprintf('', next_id + 1), + '') + rels +} + +#' Generate workbook XML. +#' @param sheet_names Character vector of sheet names. +#' @return Character vector of XML lines. +workbook_xml <- function(sheet_names) { + xml <- c('', + '', + '') + + for (i in seq_along(sheet_names)) { + xml <- c(xml, sprintf('', xml_escape(sheet_names[i]), i, i)) + } + + c(xml, '', '') +} + +#' Generate [Content_Types].xml content. +#' @param sheet_files Character vector of worksheet filenames. +#' @return Character vector of XML lines. +content_types_xml <- function(sheet_files) { + xml <- c('', + '', + '', + '', + '', + '', + '', + '', + '') + + for (file in sheet_files) { + xml <- c(xml, sprintf('', file)) + } + + c(xml, '') +} + +#' Generate docProps/core.xml content. +#' @param creator Workbook creator name. +#' @param timestamp Timestamp string in W3CDTF. +#' @return Character vector of XML lines. +core_props_xml <- function(creator, timestamp) { + c( + '', + '', + sprintf('%s', xml_escape(creator)), + sprintf('%s', xml_escape(creator)), + sprintf('%s', timestamp), + sprintf('%s', timestamp), + '' + ) +} + +#' Generate docProps/app.xml content. +#' @param sheet_names Character vector of sheet names. +#' @return Character vector of XML lines. +app_props_xml <- function(sheet_names) { + sheet_count <- length(sheet_names) + titles <- paste(sprintf('%s', xml_escape(sheet_names)), collapse = '') + + c( + '', + '', + 'AlloMate', + '0', + 'false', + 'Worksheets', + sprintf('%d', sheet_count), + '', + sprintf('%s', sheet_count, titles), + '', + 'false', + 'false', + 'false', + '16.0000', + '' + ) +} + +#' Generate top-level relationships XML. +#' @return Character vector of XML lines. +root_relationships_xml <- function() { + c( + '', + '', + '', + '', + '', + '' + ) +} + +#' Write an XLSX file using only base R functionality. +#' @param file Destination file path. +#' @param sheets Named list of sheet data. Each element is coerced to a data frame. +#' @param creator Optional creator/author string. +write_xlsx_pure <- function(file, sheets, creator = "AlloMate") { + if (is.null(names(sheets)) || any(names(sheets) == "")) { + stop("All sheets must be named for write_xlsx_pure().") + } + + sheet_names <- sanitize_sheet_names(names(sheets)) + sheet_data <- lapply(sheets, coerce_sheet_data) + + string_pool <- character() + for (df in sheet_data) { + string_pool <- c(string_pool, gather_sheet_strings(df)) + } + + unique_strings <- unique(string_pool) + string_count <- length(string_pool) + + if (length(unique_strings) > 0) { + strings_map <- setNames(seq_along(unique_strings) - 1L, unique_strings) + } else { + strings_map <- setNames(integer(), character()) + } + + timestamp <- format(Sys.time(), "%Y-%m-%dT%H:%M:%SZ", tz = "UTC") + + tmp_dir <- tempfile("allomate_xlsx_") + dir.create(tmp_dir) + on.exit(unlink(tmp_dir, recursive = TRUE, force = TRUE), add = TRUE) + + dir.create(file.path(tmp_dir, "_rels")) + dir.create(file.path(tmp_dir, "docProps")) + dir.create(file.path(tmp_dir, "xl")) + dir.create(file.path(tmp_dir, "xl", "_rels")) + dir.create(file.path(tmp_dir, "xl", "worksheets")) + + sheet_files <- character(length(sheet_data)) + + for (i in seq_along(sheet_data)) { + sheet_file <- sprintf("sheet%d.xml", i) + sheet_path <- file.path(tmp_dir, "xl", "worksheets", sheet_file) + sheet_files[i] <- sheet_file + + xml_lines <- build_sheet_xml(sheet_data[[i]], strings_map) + writeLines(xml_lines, sheet_path, useBytes = TRUE) + } + + writeLines(styles_xml_content(), file.path(tmp_dir, "xl", "styles.xml"), useBytes = TRUE) + writeLines(workbook_xml(sheet_names), file.path(tmp_dir, "xl", "workbook.xml"), useBytes = TRUE) + writeLines(workbook_rels_xml(sheet_files), file.path(tmp_dir, "xl", "_rels", "workbook.xml.rels"), useBytes = TRUE) + writeLines(shared_strings_xml(unique_strings, string_count), file.path(tmp_dir, "xl", "sharedStrings.xml"), useBytes = TRUE) + writeLines(root_relationships_xml(), file.path(tmp_dir, "_rels", ".rels"), useBytes = TRUE) + writeLines(content_types_xml(sheet_files), file.path(tmp_dir, "[Content_Types].xml"), useBytes = TRUE) + writeLines(core_props_xml(creator, timestamp), file.path(tmp_dir, "docProps", "core.xml"), useBytes = TRUE) + writeLines(app_props_xml(sheet_names), file.path(tmp_dir, "docProps", "app.xml"), useBytes = TRUE) + + if (file.exists(file)) { + unlink(file) + } + + files_to_zip <- list.files(path = tmp_dir, recursive = TRUE, include.dirs = FALSE) + + if (length(files_to_zip) == 0) { + stop("No files generated for XLSX archive.") + } + + zip_temp <- tempfile(fileext = ".zip") + zip_success <- FALSE + zip_error <- NULL + + try_system_zip <- nzchar(Sys.which("zip")) || identical(.Platform$OS.type, "windows") + + if (try_system_zip) { + old_wd <- getwd() + on.exit(setwd(old_wd), add = TRUE) + setwd(tmp_dir) + + zip_success <- tryCatch({ + utils::zip(zipfile = zip_temp, files = files_to_zip, flags = "-q") + file.exists(zip_temp) + }, warning = function(w) { + zip_error <<- w$message + FALSE + }, error = function(e) { + zip_error <<- e$message + FALSE + }) + + setwd(old_wd) + } + + if (zip_success && file.exists(zip_temp)) { + if (!file.copy(zip_temp, file, overwrite = TRUE)) { + stop("Unable to write XLSX file to destination path.") + } + unlink(zip_temp) + return(invisible(file)) + } + + if (file.exists(zip_temp)) { + unlink(zip_temp) + } + + fallback_success <- tryCatch({ + write_zip_no_compress(file, tmp_dir, files_to_zip) + TRUE + }, error = function(e) { + zip_error <<- e$message + FALSE + }) + + if (!fallback_success) { + msg <- "Failed to create ZIP archive for XLSX output." + if (!is.null(zip_error)) { + msg <- paste(msg, "Details:", zip_error) + } + stop(msg) + } + + invisible(file) +} + + diff --git a/app/R/ui_helpers.R b/app/R/ui_helpers.R index 5f30376..0766ee5 100644 --- a/app/R/ui_helpers.R +++ b/app/R/ui_helpers.R @@ -32,23 +32,23 @@ create_ocs_trait_inputs <- function(n) { #' Generate package status text #' @return Formatted status text for display generate_package_status <- function() { - status_text <- "" + optisel_flag <- exists("optisel_available") && isTRUE(optisel_available) + fallback_flag <- exists("custom_ocs_available") && isTRUE(custom_ocs_available) + quadprog_flag <- exists("quadprog_available") && isTRUE(quadprog_available) - if (exists("is_webr") && is_webr) { - status_text <- paste(status_text, "🌐 WebR environment detected\n", sep = "") - } - - if (exists("optisel_available") && optisel_available) { - status_text <- paste(status_text, "✅ optiSel package is available - OCS functionality enabled", sep = "") - } else if (exists("custom_ocs_available") && custom_ocs_available) { - status_text <- paste(status_text, "✅ Custom OCS fallback is available - OCS functionality enabled", sep = "") + if (optisel_flag) { + status_text <- "✅ Optimum Contribution Selection: Ready\n📦 Current solver: optiSel" + } else if (fallback_flag && quadprog_flag) { + status_text <- "✅ Optimum Contribution Selection: Ready\n📦 Current solver: quadprog fallback" if (exists("is_webr") && is_webr) { - status_text <- paste(status_text, "\n📦 Using custom implementation (optiSel not available in WebR)", sep = "") + status_text <- paste(status_text, "\nℹ️ optiSel not available in WebR", sep = "") } else { - status_text <- paste(status_text, "\n📦 Using custom implementation (optiSel not installed)", sep = "") + status_text <- paste(status_text, "\nℹ️ optiSel not installed", sep = "") } + } else if (fallback_flag && !quadprog_flag) { + status_text <- "❌ Optimum Contribution Selection: Not Ready\n⚠️ quadprog package required for fallback\nℹ️ Install with: install.packages('quadprog')" } else { - status_text <- paste(status_text, "❌ OCS functionality not available - Neither optiSel nor fallback could be loaded", sep = "") + status_text <- "❌ Optimum Contribution Selection: Not Ready\n⚠️ Neither optiSel nor fallback available" } status_text diff --git a/app/R/utils.R b/app/R/utils.R index 0df919b..f230b1b 100644 --- a/app/R/utils.R +++ b/app/R/utils.R @@ -23,7 +23,22 @@ fallback_kinship <- function(ped) { #' @param file Uploaded file object #' @return List with candidates data frame and male/female ID vectors read_candidates <- function(file) { - df <- readr::read_table(file$datapath) + df <- readr::read_table(file$datapath, show_col_types = FALSE) + names(df) <- tolower(names(df)) + + if (!"id" %in% names(df) && "candidate" %in% names(df)) { + names(df)[names(df) == "candidate"] <- "id" + } + + required <- c("id", "sex") + missing_cols <- setdiff(required, names(df)) + if (length(missing_cols) > 0) { + stop(sprintf("CANDIDATES: missing required column(s): %s", paste(missing_cols, collapse = ", "))) + } + + df$id <- as.character(df$id) + df$sex <- toupper(as.character(df$sex)) + list( candidates = df, males = filter(df, sex == "M") %>% pull(id), @@ -34,8 +49,22 @@ read_candidates <- function(file) { #' Clean and validate pedigree data #' @param ped Raw pedigree data frame #' @return Cleaned pedigree object for kinship calculation -clean_pedigree <- function(ped) { - final_ped <- ped %>% +clean_pedigree <- function(ped, return_stats = FALSE) { + ped_chr <- ped %>% + mutate(across(c(id, sire, dam), as.character)) + + is_missing_parent <- function(x) { + is.na(x) | x == "" | x == "0" + } + + total_records <- nrow(ped_chr) + unknown_parent_count <- sum(is_missing_parent(ped_chr$sire) | is_missing_parent(ped_chr$dam), na.rm = TRUE) + circular_rows <- (ped_chr$id == ped_chr$sire) | (ped_chr$id == ped_chr$dam) + circular_reference_count <- sum(circular_rows, na.rm = TRUE) + # Count only the extra duplicates (not the first occurrence we'll keep) + duplicates_removed <- sum(duplicated(ped_chr$id)) + + final_ped <- ped_chr %>% mutate(across(c(id, sire, dam), as.factor)) %>% mutate(sex = case_when(id %in% sire ~ 0, id %in% dam ~ 1, TRUE ~ 2)) %>% { @@ -47,9 +76,8 @@ clean_pedigree <- function(ped) { parents_fixed } %>% { - # Remove duplicates (same as original) - doubled <- table(.$id)[table(.$id) > 1] %>% names() - .[!.$id %in% doubled, ] + # Remove duplicate rows, keeping first occurrence + .[!duplicated(.$id), ] } %>% { # Remove circular dependencies (same as original) @@ -70,7 +98,17 @@ clean_pedigree <- function(ped) { } else { fallback_pedigree(id, dadid, momid, sex, missid = "0") }) - + + if (return_stats) { + stats <- list( + records_loaded = total_records, + unknown_parent_count = unknown_parent_count, + circular_reference_count = circular_reference_count, + duplicates_removed = duplicates_removed + ) + return(list(pedigree = final_ped, stats = stats)) + } + return(final_ped) } diff --git a/app/global.R b/app/global.R index fe0912c..6d1a027 100644 --- a/app/global.R +++ b/app/global.R @@ -17,17 +17,15 @@ is_shiny_server <- grepl("^/home/web_user/", getwd()) || grepl("^/tmp/", getwd() # Load non-tidyverse packages first required_packages <- c("shiny", "shinyjs", "DT", "openxlsx") -# Try to load quadprog (optional - we have fallback) -quadprog_available <- FALSE -tryCatch({ - if (!require(quadprog, quietly = TRUE)) { - install.packages("quadprog", repos = "https://cran.rstudio.com/", quiet = TRUE) - } - library(quadprog, quietly = TRUE) - quadprog_available <- TRUE -}, error = function(e) { - warning(paste("Package quadprog not available (will use fallback optimization):", e$message)) -}) +openxlsx_available <- FALSE + +# Try to detect quadprog (optional - we have fallback) +quadprog_available <- requireNamespace("quadprog", quietly = TRUE) +if (quadprog_available) { + message("✅ quadprog available") +} else { + warning("Package quadprog not available (fallback optimization will be unavailable).") +} # Try to load kinship2 (optional - we have fallback) kinship2_available <- FALSE @@ -41,13 +39,30 @@ tryCatch({ warning(paste("Package kinship2 not available (will use fallback):", e$message)) }) +# Try to load lpSolve for transportation-based mating +tryCatch({ + if (!require(lpSolve, quietly = TRUE)) { + install.packages("lpSolve", repos = "https://cran.rstudio.com/", quiet = TRUE) + } + library(lpSolve, quietly = TRUE) + message("✅ lpSolve loaded successfully") +}, error = function(e) { + warning(paste("Package lpSolve not available (transportation mating will be unavailable):", e$message)) +}) + # Load non-tidyverse packages first for (pkg in required_packages) { tryCatch({ library(pkg, character.only = TRUE) message(paste("✅ Loaded package:", pkg)) + if (identical(pkg, "openxlsx")) { + openxlsx_available <<- TRUE + } }, error = function(e) { warning(paste("Package", pkg, "not available:", e$message)) + if (identical(pkg, "openxlsx")) { + openxlsx_available <<- FALSE + } }) } @@ -60,12 +75,30 @@ tryCatch({ }) # Detect WebR environment -is_webr <- exists("webr") && !is.null(webr) +is_webr <- (exists("webr") && !is.null(webr)) || + grepl("emscripten|wasm", tolower(R.version$platform)) + +# Configure defaults for webR runtime quirks +if (is_webr) { + options(allomate.force_greedy_mating = TRUE) + options(allomate.force_qp_greedy = FALSE) + message("🌐 webR environment detected; enabling greedy mating fallback.") +} else { + options(allomate.force_greedy_mating = FALSE) + options(allomate.force_qp_greedy = FALSE) +} # Try to install and load optiSel - with custom fallback optisel_available <- FALSE custom_ocs_available <- FALSE +functions_path <- if (app_dir) "R/load_functions.R" else "app/R/load_functions.R" +if (!file.exists(functions_path)) { + stop(paste("Functions file not found at:", functions_path)) +} + +source(functions_path) + tryCatch({ if (!require(optiSel, quietly = TRUE)) { install.packages("optiSel", repos = "https://cran.rstudio.com/", quiet = TRUE) @@ -74,34 +107,13 @@ tryCatch({ optisel_available <- TRUE message("✅ optiSel loaded successfully") }, error = function(e) { - message("⚠️ optiSel not available - loading custom OCS fallback") - - # Load all organized functions (which includes the fallback) - tryCatch({ - # Determine correct path based on current directory - if (app_dir) { - functions_path <- "R/load_functions.R" - } else { - functions_path <- "app/R/load_functions.R" - } - - # Check if the functions file exists before sourcing - if (!file.exists(functions_path)) { - stop(paste("Functions file not found at:", functions_path)) - } - - source(functions_path) - - # Check if the fallback functions were loaded and flag was set - if (exists("custom_ocs_available") && custom_ocs_available) { - message("✅ Custom OCS fallback loaded successfully") - message("📦 OCS functionality enabled via custom fallback") - } else { - message("⚠️ Custom OCS fallback not available after loading functions") - } - }, error = function(func_error) { - message("❌ Could not load organized functions:", func_error$message) - }) + message("⚠️ optiSel not available - using custom OCS fallback") + if (exists("custom_ocs_available") && custom_ocs_available) { + message("✅ Custom OCS fallback loaded successfully") + message("📦 OCS functionality enabled via custom fallback") + } else { + message("⚠️ Custom OCS fallback not available after loading functions") + } }) # Set global flags for app behavior diff --git a/app/rsconnect/shinyapps.io/tnst4x-aj0ackerman/allomate.dcf b/app/rsconnect/shinyapps.io/tnst4x-aj0ackerman/allomate.dcf index ee9df06..b50d280 100644 --- a/app/rsconnect/shinyapps.io/tnst4x-aj0ackerman/allomate.dcf +++ b/app/rsconnect/shinyapps.io/tnst4x-aj0ackerman/allomate.dcf @@ -5,6 +5,6 @@ account: tnst4x-aj0ackerman server: shinyapps.io hostUrl: https://api.shinyapps.io/v1 appId: 15539609 -bundleId: 10943464 +bundleId: 11136719 url: https://tnst4x-aj0ackerman.shinyapps.io/allomate/ version: 1 diff --git a/app/server.R b/app/server.R index 9912783..2abd0e9 100644 --- a/app/server.R +++ b/app/server.R @@ -265,10 +265,186 @@ server <- function(input, output, session) { }) trait_counter <- reactiveVal(1) + candidate_status <- reactiveVal(list(ok = FALSE, error = NULL)) ocs_results_reactive <- reactiveVal() # For OCS results ebv_results_reactive <- reactiveVal() # For EBV results error_message <- reactiveVal("") # For tracking errors + pedigree_validation_stats <- reactiveVal(NULL) + # Export: Prebuild workbook so Chrome downloads reliably + export_cache <- reactiveVal(NULL) + + generate_export_xlsx <- function(dest_file) { + use_openxlsx <- exists("openxlsx_available") && isTRUE(openxlsx_available) + safe_char <- function(x) if (is.null(x) || length(x) == 0) NA_character_ else as.character(x) + + readme_text <- c( + "📊 AlloMate Complete Results Report", + "", + "This Excel file contains all results from your AlloMate analysis:", + "", + "📋 Worksheets included:", + "1. README - This overview and explanation", + "2. Filtered Results - Crosses meeting criteria (positive EBVs, kinship below threshold)", + "3. EBV Matrix - Complete matrix view with masked values", + "4. OCS Candidates - Selected candidates from Optimum Contribution Selection", + "5. Mating Plan - Recommended mating pairs from OCS", + "6. Parameters - Analysis parameters used", + "", + "🔍 Data Details:", + "- Filtered Results: Only crosses with positive EBVs and kinship below threshold", + "- EBV Matrix: All possible crosses with values masked for failed criteria", + "- OCS Results: Available only after running Optimum Contribution Selection", + "", + "📅 Generated on:", as.character(Sys.Date()) + ) + + ebv_results <- ebv_results_reactive() + mat_for_excel <- NULL + filtered_results_df <- NULL + ebv_matrix_df <- NULL + + if (!is.null(ebv_results)) { + m_ids <- unique(ebv_results$full_results$Male) + f_ids <- unique(ebv_results$full_results$Female) + mat_for_excel <- matrix(NA_real_, nrow = length(m_ids), ncol = length(f_ids), + dimnames = list(m_ids, f_ids)) + + if (nrow(ebv_results$filt_results_matrix) > 0) { + for (i in seq_len(nrow(ebv_results$filt_results_matrix))) { + m <- ebv_results$filt_results_matrix$Male[i] + f <- ebv_results$filt_results_matrix$Female[i] + val <- ebv_results$filt_results_matrix$EBV[i] + if (!is.na(m) && !is.na(f) && + m %in% rownames(mat_for_excel) && f %in% colnames(mat_for_excel)) { + mat_for_excel[m, f] <- val + } + } + } + + filtered_results_df <- as.data.frame(ebv_results$filt_results_table) + + row_labels <- rownames(mat_for_excel) + if (is.null(row_labels)) row_labels <- seq_len(nrow(mat_for_excel)) + ebv_matrix_df <- data.frame( + Male = row_labels, + mat_for_excel, + check.names = FALSE, + stringsAsFactors = FALSE + ) + } + + formatted_results <- NULL + if (!is.null(ocs_results_reactive())) { + formatted_results <- tryCatch( + format_ocs_results(ocs_results_reactive()), + error = function(e) { + message("format_ocs_results error: ", e$message) + NULL + } + ) + } + + params_data <- data.frame( + Parameter = c("Kinship Threshold", "Desired Inbreeding Rate", "Number of Offspring", "Analysis Date"), + Value = c( + safe_char(input$thresh), + safe_char(input$inbreeding_rate), + safe_char(input$num_offspring), + as.character(Sys.Date()) + ), + stringsAsFactors = FALSE + ) + + ok <- FALSE + tryCatch({ + if (use_openxlsx) { + wb <- openxlsx::createWorkbook() + + openxlsx::addWorksheet(wb, "README") + openxlsx::writeData(wb, "README", readme_text) + + if (!is.null(ebv_results)) { + openxlsx::addWorksheet(wb, "Filtered Results") + openxlsx::writeData(wb, "Filtered Results", filtered_results_df, rowNames = TRUE) + + openxlsx::addWorksheet(wb, "EBV Matrix") + openxlsx::writeData(wb, "EBV Matrix", ebv_matrix_df, rowNames = FALSE) + } + + if (!is.null(formatted_results)) { + openxlsx::addWorksheet(wb, "OCS Candidates") + openxlsx::writeData(wb, "OCS Candidates", formatted_results$candidate_table, rowNames = FALSE) + + openxlsx::addWorksheet(wb, "Mating Plan") + openxlsx::writeData(wb, "Mating Plan", formatted_results$mating_table, rowNames = FALSE) + } + + openxlsx::addWorksheet(wb, "Parameters") + openxlsx::writeData(wb, "Parameters", params_data, rowNames = FALSE) + + openxlsx::saveWorkbook(wb, dest_file, overwrite = TRUE) + } else { + sheets <- list( + README = data.frame(Text = readme_text, stringsAsFactors = FALSE), + Parameters = params_data + ) + if (!is.null(filtered_results_df)) sheets$`Filtered Results` <- filtered_results_df + if (!is.null(ebv_matrix_df)) sheets$`EBV Matrix` <- ebv_matrix_df + if (!is.null(formatted_results)) { + sheets$`OCS Candidates` <- as.data.frame(formatted_results$candidate_table) + sheets$`Mating Plan` <- as.data.frame(formatted_results$mating_table) + } + write_xlsx_pure(dest_file, sheets) + } + ok <- TRUE + }, error = function(e) { + message("Export error: ", e$message) + writeLines(paste("Export failed:", e$message), dest_file, useBytes = TRUE) + }) + ok + } + + observeEvent({ + list( + ebv_results_reactive(), + ocs_results_reactive(), + input$thresh, + input$inbreeding_rate, + input$num_offspring + ) + }, { + tmp <- tempfile(pattern = "allomate_export_", fileext = ".xlsx") + if (generate_export_xlsx(tmp) && file.exists(tmp) && file.info(tmp)$size > 0) { + old <- export_cache() + export_cache(tmp) + if (!is.null(old) && file.exists(old)) unlink(old, force = TRUE) + } else if (file.exists(tmp)) { + unlink(tmp, force = TRUE) + } + }, ignoreNULL = FALSE) + + session$onSessionEnded(function() { + cache <- export_cache() + if (!is.null(cache) && file.exists(cache)) unlink(cache, force = TRUE) + }) + + output$download_all_results <- downloadHandler( + filename = function() { + paste0("AlloMate_Complete_Results-", Sys.Date(), ".xlsx") + }, + content = function(file) { + cache <- export_cache() + if (!is.null(cache) && file.exists(cache)) { + file.copy(cache, file, overwrite = TRUE) + } else { + generate_export_xlsx(file) + } + }, + contentType = "application/vnd.openxmlformats-officedocument.spreadsheetml.sheet" + ) + outputOptions(output, "download_all_results", suspendWhenHidden = FALSE) + # Dynamic startup guide output$dynamic_guide <- renderUI({ # Check for errors first @@ -289,23 +465,37 @@ server <- function(input, output, session) { has_kinship_threshold <- !is.null(input$thresh) && input$thresh < 1 has_traits <- trait_counter() > 0 && any(sapply(1:trait_counter(), function(i) !is.null(input[[paste0("trait_file_", i)]]))) has_ocs_results <- !is.null(ocs_results_reactive()) + + # Determine error type from error message, if not working properly, update the grep logic + pedigree_error <- current_error != "" && grepl("error|pedigree|kinship", current_error, ignore.case = TRUE) + trait_error <- current_error != "" && grepl("error|ebv|trait|weight", current_error, ignore.case = TRUE) + ocs_error <- current_error != "" && grepl("error|ocs", current_error, ignore.case = TRUE) + + cs <- candidate_status() + candidate_ready <- isTRUE(cs$ok) + candidate_error_flag <- !is.null(cs$error) + + # Helper to decide icon per step + get_step_icon <- function(has_item, error_flag) { + if (error_flag) { + "❌" + } else if (has_item) { + "✅" + } else { + "⬜" + } + } # Build step-by-step guide steps <- c() # Step 1: Candidates - if (has_candidates) { - steps <- c(steps, "

✅ Step 1: Upload your candidate list to begin the analysis

") - } else { - steps <- c(steps, "

Step 1: Upload your candidate list to begin the analysis

") - } + step1_icon <- get_step_icon(candidate_ready, candidate_error_flag) + steps <- c(steps, sprintf("

%s Step 1: Upload your candidate list to begin the analysis

", step1_icon)) # Step 2: Pedigree - if (has_pedigree) { - steps <- c(steps, "

✅ Step 2: Upload your pedigree file for kinship calculations

") - } else { - steps <- c(steps, "

Step 2: Upload your pedigree file for kinship calculations

") - } + step2_icon <- get_step_icon(has_pedigree, pedigree_error) + steps <- c(steps, sprintf("

%s Step 2: Upload your pedigree file for kinship calculations

", step2_icon)) # Step 3: Kinship threshold if (has_pedigree) { @@ -315,18 +505,12 @@ server <- function(input, output, session) { } # Step 4: Traits - if (has_traits) { - steps <- c(steps, "

✅ Step 4: Add trait files and weights for breeding value analysis

") - } else { - steps <- c(steps, "

Step 4: Add trait files and weights for breeding value analysis

") - } + step4_icon <- get_step_icon(has_traits, trait_error) + steps <- c(steps, sprintf("

%s Step 4: Add trait files and weights for breeding value analysis

", step4_icon)) # Step 5: OCS - if (has_ocs_results) { - steps <- c(steps, "

✅ Step 5: Configure OCS parameters and run analysis

") - } else { - steps <- c(steps, "

Step 5: Configure OCS parameters and run analysis

") - } + step5_icon <- get_step_icon(has_ocs_results, ocs_error) + steps <- c(steps, sprintf("

%s Step 5: Configure OCS parameters and run analysis

", step5_icon)) # Add completion message if all steps are done if (has_ocs_results) { @@ -340,6 +524,188 @@ server <- function(input, output, session) { HTML(paste(steps, collapse = "")) }) + # File status display + output$file_status_display <- renderUI({ + # Check what files/data are ready + has_candidates <- !is.null(input$candidate_file) + has_pedigree <- !is.null(input$pedigree_file) + has_ebv <- !is.null(ebv_results_reactive()) + has_ocs <- !is.null(ocs_results_reactive()) + + # Check for errors + current_error <- error_message() + has_error <- current_error != "" + + cs <- candidate_status() + candidate_ready <- isTRUE(cs$ok) + candidate_error_flag <- !is.null(cs$error) + + # Determine error type from error message + pedigree_error <- has_error && grepl("pedigree|kinship", current_error, ignore.case = TRUE) + ebv_error <- has_error && grepl("ebv|trait|weight", current_error, ignore.case = TRUE) + ocs_error <- has_error && grepl("ocs", current_error, ignore.case = TRUE) + + # Helper function to determine status icon and text + get_status <- function(has_data, uploaded, has_specific_error, ready_text = "Ready", pending_text = "Not uploaded", error_text = "Error") { + if (has_specific_error) { + list(icon = "❌", text = paste0("", error_text, "")) + } else if (has_data) { + list(icon = "✅", text = paste0("", ready_text, "")) + } else if (uploaded) { + list(icon = "⬜", text = paste0("Processing...")) + } else { + list(icon = "⬜", text = paste0("", pending_text, "")) + } + } + + # Candidate status with explicit tracking + candidate_status_ui <- if (candidate_error_flag) { + list(icon = "❌", text = "Error") + } else if (candidate_ready) { + list(icon = "✅", text = "Ready") + } else if (has_candidates) { + list(icon = "⬜", text = "Processing...") + } else { + list(icon = "⬜", text = "Not uploaded") + } + + # Pedigree status + pedigree_status <- get_status(has_pedigree, has_pedigree, pedigree_error, + "Ready", "Not uploaded", "Error") + + # EBV status - can error if uploaded but processing failed + ebv_status <- get_status(has_ebv, + has_candidates && has_pedigree, + ebv_error, + "Ready", "Pending", "Error") + + # OCS status + ocs_status <- get_status(has_ocs, FALSE, ocs_error, + "Ready", "Not run", "Error") + + # Create status lines with emojis + status_lines <- c( + paste0( + candidate_status_ui$icon, + " Candidate List: ", + candidate_status_ui$text + ), + paste0( + pedigree_status$icon, + " Pedigree Data: ", + pedigree_status$text + ), + paste0( + ebv_status$icon, + " EBV Matrix: ", + ebv_status$text + ), + paste0( + ocs_status$icon, + " OCS Results: ", + ocs_status$text + ) + ) + + # Check if all files are ready for download + all_ready <- candidate_ready && has_pedigree && has_ebv + + download_status <- if (has_error) { + "
+

⚠️ Error detected. Check the startup guide above for details.

+
" + } else if (all_ready) { + "
+

✅ Download ready! All core data is available.

+
" + } else { + "
+

Upload required files to enable download.

+
" + } + + HTML(paste0( + "
", + paste(status_lines, collapse = "
"), + "
", + download_status + )) + }) + + output$pedigree_status_display <- renderUI({ + stats <- pedigree_validation_stats() + if (is.null(stats)) { + return(NULL) + } + + get_count <- function(val) { + if (is.null(val) || is.na(val)) { + 0L + } else { + as.integer(val) + } + } + + format_count <- function(val) { + format(get_count(val), big.mark = ",", scientific = FALSE) + } + + records <- format_count(stats$records_loaded) + unknown_count <- get_count(stats$unknown_parent_count) + circular_count <- get_count(stats$circular_reference_count) + missing_count <- get_count(stats$missing_candidates) + duplicates <- get_count(stats$duplicates_removed) + + green_box <- paste0( + "
", + "✅ ", records, " records loaded", + "
" + ) + + # Build yellow box warnings only for non-zero counts + yellow_warnings <- c() + if (unknown_count > 0) { + yellow_warnings <- c(yellow_warnings, + paste0("

⚠️ ", + format_count(stats$unknown_parent_count), + " individuals with unknown parent(s) (treated as founders)

")) + } + if (circular_count > 0) { + yellow_warnings <- c(yellow_warnings, + paste0("

⚠️ ", + format_count(stats$circular_reference_count), + " circular references detected and broken at earliest generation

")) + } + if (missing_count > 0) { + yellow_warnings <- c(yellow_warnings, + paste0("

⚠️ ", + format_count(stats$missing_candidates), + " selection candidates missing from pedigree

")) + } + + yellow_box <- if (length(yellow_warnings) > 0) { + paste0( + "
", + paste(yellow_warnings, collapse = ""), + "
" + ) + } else { + "" + } + + red_box <- if (duplicates > 0) { + paste0( + "
", + "❌ ", format(duplicates, big.mark = ",", scientific = FALSE), " duplicates removed", + "
" + ) + } else { + "" + } + + HTML(paste0(green_box, yellow_box, red_box)) + }) + observeEvent(input$add_trait, { trait_counter(trait_counter() + 1) }) @@ -362,7 +728,23 @@ server <- function(input, output, session) { candidates_data <- reactive({ req(input$candidate_file) - read_candidates(input$candidate_file) + + tryCatch({ + res <- read_candidates(input$candidate_file) + candidate_status(list(ok = TRUE, error = NULL)) + + current_err <- error_message() + if (current_err != "" && grepl("^Error processing candidates", current_err)) { + error_message("") + } + + res + }, error = function(e) { + candidate_status(list(ok = FALSE, error = e$message)) + error_message(paste0("Error processing candidates: ", e$message)) + validate(need(FALSE, e$message)) + NULL + }) }) pedigree_data <- reactiveVal(NULL) @@ -374,7 +756,17 @@ server <- function(input, output, session) { tryCatch({ raw_ped <- readr::read_table(input$pedigree_file$datapath) - final_ped <- clean_pedigree(raw_ped) + cleaned_ped <- clean_pedigree(raw_ped, return_stats = TRUE) + final_ped <- cleaned_ped$pedigree + + # Check for candidates missing from pedigree + candidate_ids <- candidates_data()$candidates$id + pedigree_ids <- as.character(raw_ped$id) + missing_candidates <- sum(!candidate_ids %in% pedigree_ids) + + # Add missing candidates count to stats + cleaned_ped$stats$missing_candidates <- missing_candidates + kinship_res <- compute_kinship_matrix(final_ped, males, females) output$quadrants_table <- DT::renderDT({ @@ -385,11 +777,13 @@ server <- function(input, output, session) { }) pedigree_data(list(results = kinship_res$results, quads = kinship_res$quads)) + pedigree_validation_stats(cleaned_ped$stats) error_message("") # Clear any previous errors output$message1 <- renderText("✅ Kinship matrix generated successfully.") }, error = function(e) { error_message(paste0("Error processing pedigree: ", e$message)) + pedigree_validation_stats(NULL) output$message1 <- renderText( paste0("❌ Error processing pedigree: Make sure your pedigree is clean and valid.\n", "Original error: ", e$message) @@ -487,90 +881,19 @@ server <- function(input, output, session) { "✅ EBV matrix generated successfully." ) - output$download_all_results <- downloadHandler( - filename = function() { - paste0("AlloMate_Complete_Results-", Sys.Date(), ".xlsx") - }, - content = function(file) { - wb <- openxlsx::createWorkbook() - - # Add comprehensive README worksheet - openxlsx::addWorksheet(wb, "README") - readme_text <- c( - "📊 AlloMate Complete Results Report", - "", - "This Excel file contains all results from your AlloMate analysis:", - "", - "📋 Worksheets included:", - "1. README - This overview and explanation", - "2. Filtered Results - Crosses meeting criteria (positive EBVs, kinship below threshold)", - "3. EBV Matrix - Complete matrix view with masked values", - "4. OCS Candidates - Selected candidates from Optimum Contribution Selection", - "5. Mating Plan - Recommended mating pairs from OCS", - "6. Parameters - Analysis parameters used", - "", - "🔍 Data Details:", - "- Filtered Results: Only crosses with positive EBVs and kinship below threshold", - "- EBV Matrix: All possible crosses with values masked for failed criteria", - "- OCS Results: Available only after running Optimum Contribution Selection", - "", - "📅 Generated on:", as.character(Sys.Date()) - ) - openxlsx::writeData(wb, "README", readme_text) - - # Add kinship/EBV results if available - ebv_results <- ebv_results_reactive() - if (!is.null(ebv_results)) { - openxlsx::addWorksheet(wb, "Filtered Results") - openxlsx::writeData(wb, "Filtered Results", ebv_results$filt_results_table, rowNames = TRUE) - - openxlsx::addWorksheet(wb, "EBV Matrix") - m_ids <- unique(ebv_results$full_results$Male) - f_ids <- unique(ebv_results$full_results$Female) - mat_for_excel <- matrix(NA_real_, nrow = length(m_ids), ncol = length(f_ids), - dimnames = list(m_ids, f_ids)) - - for (i in seq_len(nrow(ebv_results$filt_results_matrix))) { - m <- ebv_results$filt_results_matrix$Male[i] - f <- ebv_results$filt_results_matrix$Female[i] - val <- ebv_results$filt_results_matrix$EBV[i] - mat_for_excel[m, f] <- val - } - openxlsx::writeData(wb, "EBV Matrix", mat_for_excel, rowNames = TRUE) - } - - # Add OCS results if available - if (exists("ocs_results_reactive") && !is.null(ocs_results_reactive())) { - results <- ocs_results_reactive() - formatted_results <- format_ocs_results(results) - - openxlsx::addWorksheet(wb, "OCS Candidates") - openxlsx::writeData(wb, "OCS Candidates", formatted_results$candidate_table, rowNames = FALSE) - - openxlsx::addWorksheet(wb, "Mating Plan") - openxlsx::writeData(wb, "Mating Plan", formatted_results$mating_table, rowNames = FALSE) - } - - # Add parameters worksheet - openxlsx::addWorksheet(wb, "Parameters") - params_data <- data.frame( - Parameter = c("Kinship Threshold", "Desired Inbreeding Rate", "Number of Offspring", "Analysis Date"), - Value = c( - as.character(input$thresh), - as.character(input$inbreeding_rate), - as.character(input$num_offspring), - as.character(Sys.Date()) - ) - ) - openxlsx::writeData(wb, "Parameters", params_data, rowNames = FALSE) - - openxlsx::saveWorkbook(wb, file, overwrite = TRUE) - } - ) + # download handler defined at top-level above }) #### OCS Server Logic #### + observeEvent(input$force_greedy_mating, { + options(allomate.force_greedy_mating = isTRUE(input$force_greedy_mating)) + }, ignoreNULL = FALSE) + + observeEvent(input$force_qp_greedy, { + options(allomate.force_qp_greedy = isTRUE(input$force_qp_greedy)) + }, ignoreNULL = FALSE) + observeEvent(input$run_ocs_btn, { req(input$pedigree_file, input$candidate_file) @@ -586,6 +909,9 @@ server <- function(input, output, session) { return(NULL) } + options(allomate.force_greedy_mating = isTRUE(input$force_greedy_mating)) + options(allomate.force_qp_greedy = isTRUE(input$force_qp_greedy)) + tryCatch({ ped_data <- read.table(input$pedigree_file$datapath, header = TRUE) candidates <- read.table(input$candidate_file$datapath, header = TRUE) @@ -610,6 +936,9 @@ server <- function(input, output, session) { joint_ebvs <- calculate_index(ebv_result$joint_ebvs, ebv_result$rel_weights) candidates <- left_join(candidates, joint_ebvs, by = c("id" = "ID")) + # Set seed for reproducibility when testing optiSel vs fallback + set.seed(42) + results <- run_ocs( candidates_df = candidates, kinship_matrix = kinship_matrix, @@ -653,6 +982,19 @@ server <- function(input, output, session) { ) }) + output$ocs_solver_note <- renderUI({ + req(ocs_results_reactive()) + formatted_results <- format_ocs_results(ocs_results_reactive()) + info <- formatted_results$summary_stats$mating_info + if (is.null(info) || is.na(info) || info == "") { + return(NULL) + } + div( + style = "margin: 10px 0; padding: 10px; background-color: #fff8e1; border-left: 4px solid #ffb300; font-size: 13px;", + tags$strong("Solver note: "), info + ) + }) + #### R Code Display and Download #### # HTML escape function for code display diff --git a/app/ui.R b/app/ui.R index 0124103..9749266 100644 --- a/app/ui.R +++ b/app/ui.R @@ -2,6 +2,13 @@ library(shiny) library(shinyjs) library(DT) +# Workaround for Chromium Issue 468227 (Chrome download attribute bug) +downloadButton <- function(...) { + tag <- shiny::downloadButton(...) + tag$attribs$download <- NULL + tag +} + ui <- function(request) { fluidPage( useShinyjs(), # enables shinyjs functions @@ -230,7 +237,7 @@ ui <- function(request) { var $target = $(target); if ($target.length) { $('.help-content').animate({ - scrollTop: $target.offset().top - 100 + scrollTop: $target.offset().top - 200 }, 800); } }); @@ -250,8 +257,8 @@ ui <- function(request) { # Logos banner (left) tags$img( - src = "logos.png", - style = "width: 70%; height: auto;" + src = "logos2.png", + style = "width: 67%; height: auto;" ) ), @@ -264,7 +271,6 @@ ui <- function(request) { style = "background-color: #ffffff; border: 1px solid #dee2e6; padding: 10px; margin-bottom: 15px; border-radius: 5px;", h4("🚀 Getting Started"), htmlOutput("dynamic_guide"), - verbatimTextOutput("package_status_text"), conditionalPanel( condition = "output.webr_detected", div( @@ -291,6 +297,7 @@ ui <- function(request) { h5("Calculate kinship matrix"), fileInput("pedigree_file", "Upload pedigree file", accept = ".txt"), + uiOutput("pedigree_status_display"), h5("Set kinship threshold"), numericInput("thresh", "Max kinship allowed between mates:", @@ -318,11 +325,30 @@ ui <- function(request) { h4("🎯 Optimum Contribution Selection", style = "color: #856404; margin-bottom: 15px; border-bottom: 1px solid #ffeaa7; padding-bottom: 8px;"), p("Configure breeding objectives and constraints:", style = "color: #6c757d; font-size: 12px; margin-bottom: 15px;"), + verbatimTextOutput("package_status_text"), + h5("Breeding Objectives"), numericInput("inbreeding_rate", "Desired Inbreeding Rate", value = 0.05, min = 0.01, max = 0.2, step = 0.01), numericInput("num_offspring", "Number of Offspring", value = 100, min = 10, step = 1), + checkboxInput( + "force_greedy_mating", + "Use greedy mating (browser-safe)", + value = isTRUE(getOption("allomate.force_greedy_mating", FALSE)) + ), + checkboxInput( + "force_qp_greedy", + "Bypass quadprog (heuristic contributions)", + value = isTRUE(getOption("allomate.force_qp_greedy", FALSE)) + ), + helpText( + style = "color: #6c757d; font-size: 11px;", + strong("Greedy mating:"), " Browser-safe, near-optimal results.", + br(), + strong("Bypass quadprog:"), " ", + span(style = "color: #dc3545;", "⚠️ Testing only—ignores inbreeding constraints!") + ), actionButton("run_ocs_btn", "Run OCS", style = "margin-top: 15px; width: 100%; background-color: #856404; color: white; border: none; padding: 10px; border-radius: 5px;") @@ -334,6 +360,14 @@ ui <- function(request) { p("Download all results in a single Excel file with multiple tabs:", style = "color: #6c757d; font-size: 12px; margin-bottom: 15px;"), downloadButton("download_all_results", "📥 Export All Results", style = "width: 100%; background-color: #28a745; color: white; border: none; padding: 10px; border-radius: 5px; margin-bottom: 10px;"), + + # File status display box + div( + style = "background-color: #f8f9fa; border: 1px solid #dee2e6; padding: 12px; margin-bottom: 10px; border-radius: 5px;", + h5("📋 File Status", style = "color: #495057; margin-top: 0; margin-bottom: 10px; font-size: 14px;"), + htmlOutput("file_status_display") + ), + actionButton("view_r_code_btn", "📝 View R Code", style = "width: 100%; background-color: #17a2b8; color: white; border: none; padding: 10px; border-radius: 5px;") ) @@ -351,6 +385,7 @@ ui <- function(request) { ), tabPanel("Optimum Contribution Selection", DTOutput("ocs_candidate_table"), + uiOutput("ocs_solver_note"), br(), DTOutput("ocs_mating_table") ), diff --git a/app/www/logos2.png b/app/www/logos2.png new file mode 100644 index 0000000..6c177bd Binary files /dev/null and b/app/www/logos2.png differ diff --git a/docs/AlloMate/app.json b/docs/AlloMate/app.json index bd5a015..4538178 100644 --- a/docs/AlloMate/app.json +++ b/docs/AlloMate/app.json @@ -1 +1 @@ -[{"name":"server.R","content":"# All packages are loaded in global.R\n# optiSel availability is checked in global.R\n\n# All functions are now loaded from the organized functions folder\n# See functions/load_all_functions.R for details\n\n#### Server Function ####\n\nserver <- function(input, output, session) {\n \n # Package status output\n output$package_status_text <- renderText({\n generate_package_status()\n })\n \n # WebR detection output\n output$webr_detected <- renderText({\n if (is_webr_environment()) {\n \"WebR detected\"\n } else {\n \"\"\n }\n })\n outputOptions(output, \"webr_detected\", suspendWhenHidden = FALSE)\n \n # Help button functionality\n observeEvent(input$help_btn, {\n updateTabsetPanel(session, \"main_tabs\", selected = \"Help\")\n })\n \n # View R Code button functionality\n observeEvent(input$view_r_code_btn, {\n updateTabsetPanel(session, \"main_tabs\", selected = \"R Code\")\n })\n \n # Back to top functionality\n observeEvent(input$back_to_top, {\n runjs(\"document.querySelector('.help-content').scrollTop = 0;\")\n })\n \n # Enhanced markdown to HTML conversion function\n markdown_to_html <- function(markdown_text) {\n # Split into lines for processing\n lines <- strsplit(markdown_text, \"\\n\")[[1]]\n html_lines <- character(length(lines))\n in_code_block <- FALSE\n in_list <- FALSE\n list_type <- \"\"\n \n for (i in seq_along(lines)) {\n line <- lines[i]\n \n # Handle code blocks\n if (grepl(\"^```\", line)) {\n if (!in_code_block) {\n html_lines[i] <- \"
\"\n          in_code_block <- TRUE\n        } else {\n          html_lines[i] <- \"<\/code><\/pre>\"\n          in_code_block <- FALSE\n        }\n        next\n      }\n      \n      if (in_code_block) {\n        html_lines[i] <- paste0(\"\", line, \"<\/span>\")\n        next\n      }\n      \n      # Handle headers\n      if (grepl(\"^# \", line)) {\n        level <- nchar(gsub(\"^(#+).*\", \"\\\\1\", line))\n        content <- gsub(\"^#+ \", \"\", line)\n        # Add emoji support for main headers\n        if (level == 1) {\n          content <- paste0(\"🧬 \", content)\n        } else if (level == 2) {\n          content <- paste0(\"🚀 \", content)\n        } else if (level == 3) {\n          content <- paste0(\"📁 \", content)\n        } else if (level == 4) {\n          content <- paste0(\"🛠️ \", content)\n        } else if (level == 5) {\n          content <- paste0(\"📊 \", content)\n        } else if (level == 6) {\n          content <- paste0(\"🔧 \", content)\n        }\n        # Create anchor ID for TOC linking\n        anchor <- tolower(gsub(\"[^a-zA-Z0-9\\\\s]\", \"\", content))\n        anchor <- gsub(\"\\\\s+\", \"-\", anchor)\n        \n        html_lines[i] <- paste0(\"\", content, \"<\/h\", level, \">\")\n        next\n      }\n      \n      # Handle bold and italic\n      line <- gsub(\"\\\\*\\\\*(.+?)\\\\*\\\\*\", \"\\\\1<\/strong>\", line)\n      line <- gsub(\"\\\\*(.+?)\\\\*\", \"\\\\1<\/em>\", line)\n      \n      # Handle inline code\n      line <- gsub(\"`(.+?)`\", \"\\\\1<\/code>\", line)\n      \n      # Handle links\n      line <- gsub(\"\\\\[(.+?)\\\\]\\\\((.+?)\\\\)\", \"\\\\1<\/a>\", line)\n      \n      # Handle lists\n      if (grepl(\"^[*-] \", line)) {\n        content <- gsub(\"^[*-] \", \"\", line)\n        if (!in_list) {\n          html_lines[i] <- paste0(\"