From bc97cec9a6c63cebcd74138b62a63e1be11c356e Mon Sep 17 00:00:00 2001 From: ankurjuneja Date: Fri, 4 Sep 2026 11:13:04 -0700 Subject: [PATCH 1/3] Stream kinship output per family to bound peak memory --- .../pipeline/kinship/populateKinship.r | 201 ++++++++++++++---- 1 file changed, 159 insertions(+), 42 deletions(-) diff --git a/ehr/resources/pipeline/kinship/populateKinship.r b/ehr/resources/pipeline/kinship/populateKinship.r index 8ccf8e741..21e479099 100644 --- a/ehr/resources/pipeline/kinship/populateKinship.r +++ b/ehr/resources/pipeline/kinship/populateKinship.r @@ -5,68 +5,185 @@ ## # This R script will calculate and store kinship coefficients (aka. relatedness) for all animals in the colony. This is a large, sparse matrix. -# The matrix is converted into a very long 3-column dataframe (animal1, animal2, coefficient). This dataframe is output to a TSV file, -# which is normally imported into ehr.kinship by java code in GeneticCalculationsImportTask +# The matrix is written out as a very long 3-column TSV file (animal1, animal2, coefficient), which is normally imported into +# ehr.kinship by java code in GeneticCalculationsImportTask. +# +# This runs as a child of the web server, so its peak is charged against the same memory budget as the JVM. The pedigree is therefore +# processed one family at a time and rows are streamed to disk in column blocks, which bounds the peak by the largest single family's +# dense kinship matrix rather than by the whole colony. library(kinship2) library(getopt) -library(Matrix) -library(dplyr) + +# Minimum coefficient to emit; 0 disables the filter. Enabling this silently lowers colony-wide kinship averages, because those +# queries divide by an independent population count rather than by the rows present. The useful range is 2^-5 to 2^-6. +MIN_COEFFICIENT <- 0 + +# Target size of the column block copied out of the dense matrix at a time. Bounds the streaming transients. +BLOCK_BYTES <- 64e6 + +OUTPUT_FILE <- 'kinship.txt' +TEMP_FILE <- 'kinship.txt.part' spec <- matrix(c( 'inputFile', '-f', 1, 'character' ), ncol=4, byrow=TRUE) opts <- getopt(spec, commandArgs(trailingOnly = TRUE)) -allPed <- read.table(opts$inputFile, quote="\"") -colnames(allPed)<-c('Id', 'Dam', 'Sire', 'Gender', 'Species') +# Emits both directions of every related pair in one family. Walks the dense matrix in column blocks so the index and value +# transients stay bounded no matter how large the family is. +writeFamily <- function(kin, con) +{ + ids <- colnames(kin) + if (is.null(ids)) + stop('Kinship matrix has no column names') + + n <- length(ids) + + # Equivalent to the previous "^(?!addin).*$" filter, but over the ids rather than over every emitted row + keep <- !startsWith(ids, 'addin') + + blockCols <- max(1L, as.integer(BLOCK_BYTES %/% (8 * n))) + written <- 0 -allPed$Id[allPed$Id == ""] <- NA -allPed$Dam[allPed$Dam == ""] <- NA -allPed$Sire[allPed$Sire == ""] <- NA -allPed$Gender[allPed$Gender == "" | is.na(allPed$Gender)] <- 3 # 3 = unknown + for (start in seq.int(1L, n, by = blockCols)) + { + end <- min(start + blockCols - 1L, n) -allPed$Species <- as.character(allPed$Species) -allPed$Species[is.na(allPed$Species)] <- c('Unknown') -allPed$Species <- as.factor(allPed$Species) + block <- kin[, start:end, drop = FALSE] + mask <- block != 0 + # which(arr.ind) and block[mask] both walk in column-major order, so these line up row for row + nz <- which(mask, arr.ind = TRUE, useNames = FALSE) + vals <- block[mask] + rm(block, mask) + + if (nrow(nz) == 0) + { + rm(nz, vals) + next + } + + i <- nz[, 1L] + j <- nz[, 2L] + (start - 1L) + rm(nz) + + # Upper triangle only: the matrix is symmetric, so the mirror row is emitted below rather than stored. This also drops + # the diagonal, and the importer discards self-pairs anyway. + sel <- i < j & keep[i] & keep[j] + if (MIN_COEFFICIENT > 0) + sel <- sel & vals >= MIN_COEFFICIENT + + i <- i[sel] + j <- j[sel] + x <- vals[sel] + rm(vals, sel) + + if (length(i) > 0) + { + # Both directions are required: the importer never symmetrises and the ONPRC aggregates join on Id and Id2 separately + rows <- data.frame(Id = c(ids[i], ids[j]), Id2 = c(ids[j], ids[i]), coefficient = c(x, x), stringsAsFactors = FALSE) + write.table(rows, file = con, row.names = FALSE, col.names = FALSE, quote = FALSE, sep = '\t') + written <- written + nrow(rows) + rm(rows) + } + + rm(i, j, x) + } -if (any(allPed$Species == 'Unknown')) { - print(paste0('There are ', sum(allPed$Species == 'Unknown'), ' Ids with species = Unknown')) + written } -newRecords <- NULL -for (species in unique(allPed$Species)){ - allRecordsForSpecies <- allPed[allPed$Species %in% species,] - print(paste0('Processing species: ', species, ', with ', nrow(allRecordsForSpecies), ' IDs')) - if (nrow(allRecordsForSpecies) == 1) { - print('single record, skipping') - next +# Kept in its own frame so the dense matrix is released before the next family allocates its own +processFamily <- function(fixed, idx, con) +{ + sub <- fixed[idx, ] + ped <- with(sub, pedigree(id = id, dadid = dadid, momid = momid, sex = sex, missid = 0)) + kin <- kinship(ped) + rm(sub, ped) + + writeFamily(kin, con) +} + +main <- function() +{ + allPed <- read.table(opts$inputFile, quote="\"") + colnames(allPed)<-c('Id', 'Dam', 'Sire', 'Gender', 'Species') + + allPed$Id[allPed$Id == ""] <- NA + allPed$Dam[allPed$Dam == ""] <- NA + allPed$Sire[allPed$Sire == ""] <- NA + allPed$Gender[allPed$Gender == "" | is.na(allPed$Gender)] <- 3 # 3 = unknown + + allPed$Species <- as.character(allPed$Species) + allPed$Species[is.na(allPed$Species)] <- c('Unknown') + allPed$Species <- as.factor(allPed$Species) + + if (any(allPed$Species == 'Unknown')) { + print(paste0('There are ', sum(allPed$Species == 'Unknown'), ' Ids with species = Unknown')) } - # Add missing parents for accurate kinship calculations - fixedRecords <- with(allRecordsForSpecies, fixParents(id = Id, dadid = Sire, momid = Dam, sex = Gender)) + if (file.exists(TEMP_FILE)) + unlink(TEMP_FILE) + + con <- file(TEMP_FILE, open = 'wt') + on.exit(if (!is.null(con)) close(con), add = TRUE) - # Kinship is expecting records to be sorted IAW it's own pedigree function - recordsForSpecies <- with(fixedRecords, pedigree(id = id, dadid = dadid, momid = momid, sex = sex, missid = 0)) + # The header is retained because the importer detects it by testing whether the third field equals "coefficient" + writeLines(paste('Id', 'Id2', 'coefficient', sep = '\t'), con) - temp.kin <- kinship(recordsForSpecies) + totalRows <- 0 + for (species in unique(allPed$Species)){ + allRecordsForSpecies <- allPed[allPed$Species %in% species,] + print(paste0('Processing species: ', species, ', with ', nrow(allRecordsForSpecies), ' IDs')) + if (nrow(allRecordsForSpecies) == 1) { + print('single record, skipping') + next + } - # Add rownames to make matrix symmetric, which is required downstream - rownames(temp.kin) <- colnames(temp.kin) + # Add missing parents for accurate kinship calculations + fixed <- with(allRecordsForSpecies, fixParents(id = Id, dadid = Sire, momid = Dam, sex = Gender)) + rm(allRecordsForSpecies) + + # Split into unrelated families. Kinship across families is zero and was never emitted, so the output is unchanged, + # but the dense matrix is now capped at the largest family instead of the whole species. + fam <- makefamid(fixed$id, fixed$dadid, fixed$momid) + famIds <- unique(fam[fam != 0]) # 0 means no relatives, so there is nothing but a self-pair to emit + print(paste0(' ', length(famIds), ' families, largest ', max(c(0, tabulate(fam[fam != 0]))))) + + speciesRows <- 0 + for (f in famIds) + { + idx <- which(fam == f) + if (length(idx) < 2) + next + + speciesRows <- speciesRows + processFamily(fixed, idx, con) + + # Only worth forcing for families whose matrix is large enough to be worth returning to the OS; + # a species can have hundreds of tiny families and R collects those on its own + if (length(idx) > 1000) + gc(verbose = FALSE) + } + + print(paste0(' wrote ', speciesRows, ' rows')) + totalRows <- totalRows + speciesRows + rm(fixed, fam) + } - # Convert kinship matrix to a triplet list of two ids and their coefficient - summaryDf <- as.data.frame(summary(as(temp.kin, "dgCMatrix"))) - idList <- rownames(temp.kin) - temp.tri <- data.frame(Id=idList[summaryDf$i], Id2=idList[summaryDf$j], coefficient=summaryDf$x) + close(con) + con <- NULL - # Now filter out parents added for kinship calculation - temp.tri <- dplyr::filter(temp.tri, grepl("^(?!addin).*$", Id, perl = TRUE)) - temp.tri <- dplyr::filter(temp.tri, grepl("^(?!addin).*$", Id2, perl = TRUE)) - temp.tri <- merge(temp.tri, allRecordsForSpecies[c('Id', 'Species')], by = 'Id', all.x = TRUE) + print(paste0('Total kinship records: ', totalRows)) + if (totalRows == 0) + stop('No kinship records were generated, so the output was not written. The import deletes all existing rows before reading it.') - newRecords <- dplyr::bind_rows(newRecords,temp.tri) + # Rename only once the file is known to be complete. The importer deletes every existing row before reading it and only checks + # that it has three lines, so a partial file would silently replace good data with partial data. The old file is removed first + # because file.rename does not reliably replace an existing destination on Windows; a crash in that window leaves no file at all, + # which the importer already reports as an error, rather than a truncated one, which it does not. + if (file.exists(OUTPUT_FILE)) + unlink(OUTPUT_FILE) + if (!file.rename(TEMP_FILE, OUTPUT_FILE)) + stop(paste0('Unable to rename ', TEMP_FILE, ' to ', OUTPUT_FILE)) } -# write TSV to disk -print(paste0('Total kinship records: ', nrow(newRecords))) -newRecords <- dplyr::arrange(newRecords, Id, Id2) -write.table(newRecords, file = "kinship.txt", append = FALSE, row.names = FALSE, quote = FALSE, sep = '\t') \ No newline at end of file +main() From f447facf5ac9238a9e4c65dde75e84c710e2e721 Mon Sep 17 00:00:00 2001 From: ankurjuneja Date: Fri, 4 Sep 2026 12:59:35 -0700 Subject: [PATCH 2/3] claude code review changes --- .../pipeline/kinship/populateKinship.r | 62 +++++++++++++++---- 1 file changed, 50 insertions(+), 12 deletions(-) diff --git a/ehr/resources/pipeline/kinship/populateKinship.r b/ehr/resources/pipeline/kinship/populateKinship.r index 21e479099..ff12b0b85 100644 --- a/ehr/resources/pipeline/kinship/populateKinship.r +++ b/ehr/resources/pipeline/kinship/populateKinship.r @@ -8,9 +8,14 @@ # The matrix is written out as a very long 3-column TSV file (animal1, animal2, coefficient), which is normally imported into # ehr.kinship by java code in GeneticCalculationsImportTask. # -# This runs as a child of the web server, so its peak is charged against the same memory budget as the JVM. The pedigree is therefore -# processed one family at a time and rows are streamed to disk in column blocks, which bounds the peak by the largest single family's -# dense kinship matrix rather than by the whole colony. +# This runs as a child of the web server, so its peak is charged against the same memory budget as the JVM. The saving comes from +# streaming: rows are appended to disk as each matrix is produced, rather than accumulated into one data frame for the whole colony +# and sorted at the end. +# +# The pedigree is also split into unrelated families, but that is worth little by itself. In a managed breeding colony almost every +# animal descends from a shared founder set, so the largest family is most of the species; at the three colonies measured it was 76%, +# 98% and 100% of the largest species. The split is kept because it is what allows each matrix to be written and released before the +# next is allocated, not because it meaningfully lowers the largest allocation. library(kinship2) library(getopt) @@ -18,7 +23,9 @@ library(getopt) # queries divide by an independent population count rather than by the rows present. The useful range is 2^-5 to 2^-6. MIN_COEFFICIENT <- 0 -# Target size of the column block copied out of the dense matrix at a time. Bounds the streaming transients. +# Size of the column slab copied out of the dense matrix at a time. This is not the peak added by the streaming loop: the index +# vector, the value vector and the frame handed to write.table are all derived from the slab, so a slab whose cells are mostly +# nonzero costs roughly four times this figure while it is being written. Halve it if a site needs the transients smaller. BLOCK_BYTES <- 64e6 OUTPUT_FILE <- 'kinship.txt' @@ -29,6 +36,14 @@ spec <- matrix(c( ), ncol=4, byrow=TRUE) opts <- getopt(spec, commandArgs(trailingOnly = TRUE)) +# One direction of a set of pairs. Kept separate so the frame and the character image write.table builds from it -- together the +# largest transient in the loop below -- are only as long as the pair count, rather than twice that. +writePairs <- function(a, b, x, con) +{ + write.table(data.frame(Id = a, Id2 = b, coefficient = x, stringsAsFactors = FALSE), + file = con, row.names = FALSE, col.names = FALSE, quote = FALSE, sep = '\t') +} + # Emits both directions of every related pair in one family. Walks the dense matrix in column blocks so the index and value # transients stay bounded no matter how large the family is. writeFamily <- function(kin, con) @@ -51,6 +66,13 @@ writeFamily <- function(kin, con) block <- kin[, start:end, drop = FALSE] mask <- block != 0 + + # which() skips an NA in the mask but block[mask] returns an element for it, so a single NA would shift every value + # against its ids from that point on and write wrong coefficients with no error. kinship() should never produce one; + # this turns the case into a visible failure rather than silent corruption. anyNA() scans without allocating. + if (anyNA(mask)) + stop('Kinship matrix contains NA') + # which(arr.ind) and block[mask] both walk in column-major order, so these line up row for row nz <- which(mask, arr.ind = TRUE, useNames = FALSE) vals <- block[mask] @@ -80,10 +102,9 @@ writeFamily <- function(kin, con) if (length(i) > 0) { # Both directions are required: the importer never symmetrises and the ONPRC aggregates join on Id and Id2 separately - rows <- data.frame(Id = c(ids[i], ids[j]), Id2 = c(ids[j], ids[i]), coefficient = c(x, x), stringsAsFactors = FALSE) - write.table(rows, file = con, row.names = FALSE, col.names = FALSE, quote = FALSE, sep = '\t') - written <- written + nrow(rows) - rm(rows) + writePairs(ids[i], ids[j], x, con) + writePairs(ids[j], ids[i], x, con) + written <- written + 2 * length(i) } rm(i, j, x) @@ -131,6 +152,7 @@ main <- function() writeLines(paste('Id', 'Id2', 'coefficient', sep = '\t'), con) totalRows <- 0 + totalFamilies <- 0 for (species in unique(allPed$Species)){ allRecordsForSpecies <- allPed[allPed$Species %in% species,] print(paste0('Processing species: ', species, ', with ', nrow(allRecordsForSpecies), ' IDs')) @@ -143,8 +165,7 @@ main <- function() fixed <- with(allRecordsForSpecies, fixParents(id = Id, dadid = Sire, momid = Dam, sex = Gender)) rm(allRecordsForSpecies) - # Split into unrelated families. Kinship across families is zero and was never emitted, so the output is unchanged, - # but the dense matrix is now capped at the largest family instead of the whole species. + # Split into unrelated families. Kinship across families is zero and was never emitted, so the output is unchanged. fam <- makefamid(fixed$id, fixed$dadid, fixed$momid) famIds <- unique(fam[fam != 0]) # 0 means no relatives, so there is nothing but a self-pair to emit print(paste0(' ', length(famIds), ' families, largest ', max(c(0, tabulate(fam[fam != 0]))))) @@ -156,6 +177,7 @@ main <- function() if (length(idx) < 2) next + totalFamilies <- totalFamilies + 1 speciesRows <- speciesRows + processFamily(fixed, idx, con) # Only worth forcing for families whose matrix is large enough to be worth returning to the OS; @@ -169,12 +191,28 @@ main <- function() rm(fixed, fam) } + if (totalRows == 0) + { + # Families were found but produced nothing, so something is wrong upstream of the write. Stopping leaves the previous + # kinship.txt in place, which matters because the import deletes every existing row before reading the file. + if (totalFamilies > 0) + stop(paste0('Found ', totalFamilies, ' related families but produced no kinship records, so the output was not written.')) + + # No parent-child links anywhere in the pedigree. That is a legitimate state -- demo containers, newly provisioned sites + # and colonies loaded without parentage all reach it -- and the previous version of this script succeeded there because the + # matrix diagonal padded the file past the importer's three-line minimum. Emit that same diagonal so those sites keep + # working: the importer discards self-pairs, so this correctly clears ehr.kinship instead of failing the nightly job. + # With no relationships every animal is non-inbred, so its self-kinship is exactly 0.5. + print('No parent-child relationships found in the pedigree; writing self-pairs only, which the import discards') + selfIds <- unique(allPed$Id[!is.na(allPed$Id)]) + writePairs(selfIds, selfIds, rep(0.5, length(selfIds)), con) + totalRows <- length(selfIds) + } + close(con) con <- NULL print(paste0('Total kinship records: ', totalRows)) - if (totalRows == 0) - stop('No kinship records were generated, so the output was not written. The import deletes all existing rows before reading it.') # Rename only once the file is known to be complete. The importer deletes every existing row before reading it and only checks # that it has three lines, so a partial file would silently replace good data with partial data. The old file is removed first From 09b3ef4b236e5874ccd2907e68d0da8ce67a6938 Mon Sep 17 00:00:00 2001 From: Marty Pradere Date: Wed, 16 Sep 2026 04:50:28 -0600 Subject: [PATCH 3/3] Bound kinship peak memory and remove the kinship2 dependency The nightly job runs as a child of the web server, so its peak is charged against the JVM's budget. Computing one family at a time and running the recursion in-script rather than through kinship2 cuts the measured peak 2.5x to 3.8x across five production pedigrees, lets the largest complete at all, and stops coefficients varying with the kinship2 version a site has installed. --- .../pipeline/kinship/populateKinship.r | 237 ++++++++++++++---- ehr/webapp/WEB-INF/ehr/ehrContext.xml | 2 +- 2 files changed, 195 insertions(+), 44 deletions(-) diff --git a/ehr/resources/pipeline/kinship/populateKinship.r b/ehr/resources/pipeline/kinship/populateKinship.r index ff12b0b85..418350fbf 100644 --- a/ehr/resources/pipeline/kinship/populateKinship.r +++ b/ehr/resources/pipeline/kinship/populateKinship.r @@ -8,26 +8,29 @@ # The matrix is written out as a very long 3-column TSV file (animal1, animal2, coefficient), which is normally imported into # ehr.kinship by java code in GeneticCalculationsImportTask. # -# This runs as a child of the web server, so its peak is charged against the same memory budget as the JVM. The saving comes from -# streaming: rows are appended to disk as each matrix is produced, rather than accumulated into one data frame for the whole colony -# and sorted at the end. +# This runs as a child of the web server, so its peak is charged against the same memory budget as the JVM. Two things bound it. +# Each family -- animals joined by a chain of parent-child links -- is computed as its own dense matrix and streamed to disk before +# the next is allocated. And the kinship recursion runs here rather than through kinship2::kinship(), which allocates (n+1)^2 +# doubles, copies them to n^2, and updates a whole generation at once; at a 24,000-animal family that is three times the matrix. # -# The pedigree is also split into unrelated families, but that is worth little by itself. In a managed breeding colony almost every -# animal descends from a shared founder set, so the largest family is most of the species; at the three colonies measured it was 76%, -# 98% and 100% of the largest species. The split is kept because it is what allows each matrix to be written and released before the -# next is allocated, not because it meaningfully lowers the largest allocation. -library(kinship2) +# In a managed breeding colony almost every animal descends from a shared founder set, so the largest family is most of the species +# (76% to 100% at the colonies measured). The split is what lets one matrix be released before the next, not a way to make the +# largest one small. library(getopt) # Minimum coefficient to emit; 0 disables the filter. Enabling this silently lowers colony-wide kinship averages, because those -# queries divide by an independent population count rather than by the rows present. The useful range is 2^-5 to 2^-6. +# queries divide by an independent population count rather than by the rows present: at 2^-6 it removed 62% of the rows and a +# third of the total kinship at one colony measured. MIN_COEFFICIENT <- 0 -# Size of the column slab copied out of the dense matrix at a time. This is not the peak added by the streaming loop: the index -# vector, the value vector and the frame handed to write.table are all derived from the slab, so a slab whose cells are mostly -# nonzero costs roughly four times this figure while it is being written. Halve it if a site needs the transients smaller. +# Columns of the dense matrix copied out at a time. The index and value vectors and the frame handed to write.table are all derived +# from the slab, so a mostly-nonzero slab costs roughly four times this while it is written. Smaller is not cheaper: 4e6 and 512e6 +# both raised the peak at a 25,000-animal colony. BLOCK_BYTES <- 64e6 +# Rows of one generation updated per pass in kinshipDense, bounding that pass's temporaries to CHUNK_ROWS x n doubles +CHUNK_ROWS <- 512L + OUTPUT_FILE <- 'kinship.txt' TEMP_FILE <- 'kinship.txt.part' @@ -36,6 +39,159 @@ spec <- matrix(c( ), ncol=4, byrow=TRUE) opts <- getopt(spec, commandArgs(trailingOnly = TRUE)) +# Row index of each animal's parent, 0 where the parent is unknown +parentIndex <- function(id, parent) +{ + match(parent, id, nomatch = 0L) +} + +# The conditions that fail the run: a blank id, an id with two rows, and an id recorded as both a dam and a sire or as its own parent +validatePedigree <- function(id, dam, sire) +{ + if (anyNA(id) || any(grepl('^ *$', id))) + stop('Pedigree contains a blank id') + + dup <- unique(id[duplicated(id)]) + if (length(dup) > 0) + stop(paste('Duplicate subject id:', paste(head(dup, 6), collapse = ' '))) + + both <- intersect(dam[!is.na(dam)], sire[!is.na(sire)]) + if (length(both) > 0) + stop(paste('Id recorded as both a dam and a sire:', paste(head(both, 6), collapse = ' '))) + + # Caught here rather than as a loop in generationDepth, which cannot say which animal is at fault + self <- id[(!is.na(dam) & id == dam) | (!is.na(sire) & id == sire)] + if (length(self) > 0) + stop(paste('Id recorded as its own parent:', paste(head(self, 6), collapse = ' '))) +} + +# Warned rather than fatal: autosomal kinship never reads gender, so a contradiction here cannot move a coefficient. An unknown +# gender is a gap rather than a contradiction, so only an explicitly male dam or female sire is reported. +reportParentSex <- function(species, id, dam, sire, gender) +{ + femaleSires <- unique(id[id %in% sire & !is.na(gender) & gender == 2]) + maleDams <- unique(id[id %in% dam & !is.na(gender) & gender == 1]) + + if (length(femaleSires) > 0) + warning(paste0(species, ': ', length(femaleSires), ' id(s) recorded female but used as a sire: ', + paste(head(femaleSires, 6), collapse = ' ')), immediate. = TRUE, call. = FALSE) + + if (length(maleDams) > 0) + warning(paste0(species, ': ', length(maleDams), ' id(s) recorded male but used as a dam: ', + paste(head(maleDams, 6), collapse = ' ')), immediate. = TRUE, call. = FALSE) +} + +# Appends a founder row for every parent that is referenced but has no row of its own, sires first then dams. These animals reach +# the output like any other; a parent that is simply unknown is left unknown, which the recursion handles directly. +closePedigree <- function(id, dam, sire) +{ + absentSires <- unique(sire[!is.na(sire) & !(sire %in% id)]) + id <- c(id, absentSires) + absentDams <- unique(dam[!is.na(dam) & !(dam %in% id)]) + id <- c(id, absentDams) + + added <- length(absentSires) + length(absentDams) + list(id = id, dam = c(dam, rep(NA_character_, added)), sire = c(sire, rep(NA_character_, added))) +} + +# Family label per animal: the smallest row index reachable through parent-child links. Kinship between different families is +# exactly zero, so each can be computed and written on its own. +familyIds <- function(dadIdx, momIdx) +{ + n <- length(dadIdx) + fam <- seq_len(n) + hasDad <- dadIdx > 0L + hasMom <- momIdx > 0L + + pullUp <- function(lab, has, idx) + { + if (!any(has)) + return(lab) + up <- tapply(lab[has], idx[has], min) + k <- as.integer(names(up)) + lab[k] <- pmin(lab[k], as.integer(up)) + lab + } + + repeat + { + new <- fam + new[hasDad] <- pmin(new[hasDad], new[dadIdx[hasDad]]) + new[hasMom] <- pmin(new[hasMom], new[momIdx[hasMom]]) + new <- pullUp(new, hasDad, dadIdx) + new <- pullUp(new, hasMom, momIdx) + if (identical(new, fam)) + break + fam <- new + } + + fam +} + +# Generation of each animal: founders are 0, everyone else is one more than their deepest parent. Stops if the links loop. +generationDepth <- function(dadIdx, momIdx) +{ + n <- length(dadIdx) + depth <- integer(n) + parents <- which(dadIdx == 0L & momIdx == 0L) + + for (i in seq_len(n)) + { + child <- match(momIdx, parents, nomatch = 0L) + match(dadIdx, parents, nomatch = 0L) + if (all(child == 0L)) + break + if (i == n) + stop('Impossible pedigree: an animal is its own ancestor') + parents <- which(child > 0L) + depth[parents] <- i + } + + if (any(depth == 0L & (dadIdx > 0L | momIdx > 0L))) + stop('Impossible pedigree: an animal is its own ancestor') + + depth +} + +# Dense kinship matrix for one family. Founders start at 0.5 on the diagonal; each generation's rows and then columns become the +# mean of its parents', and the diagonal becomes (1 + kinship of the parents) / 2. An unknown parent contributes a zero row, and +# (x + 0) / 2 is the same double as x / 2, so that case is written directly. Rows within a generation are independent, so they +# are updated in chunks to bound the temporaries. +kinshipDense <- function(dadIdx, momIdx, depth, chunk = CHUNK_ROWS) +{ + n <- length(dadIdx) + if (n == 1L) + return(matrix(0.5, 1L, 1L)) + + kmat <- diag(0.5, n) + + for (gen in seq_len(max(depth))) + { + indx <- which(depth == gen) + chunks <- split(indx, ceiling(seq_along(indx) / chunk)) + + for (ch in chunks) + { + mm <- momIdx[ch]; dd <- dadIdx[ch] + both <- mm != 0L & dd != 0L; mo <- mm != 0L & dd == 0L; do <- mm == 0L & dd != 0L + if (any(both)) kmat[ch[both], ] <- (kmat[mm[both], , drop = FALSE] + kmat[dd[both], , drop = FALSE]) / 2 + if (any(mo)) kmat[ch[mo], ] <- kmat[mm[mo], , drop = FALSE] / 2 + if (any(do)) kmat[ch[do], ] <- kmat[dd[do], , drop = FALSE] / 2 + } + for (ch in chunks) + { + mm <- momIdx[ch]; dd <- dadIdx[ch] + both <- mm != 0L & dd != 0L; mo <- mm != 0L & dd == 0L; do <- mm == 0L & dd != 0L + if (any(both)) kmat[, ch[both]] <- (kmat[, mm[both], drop = FALSE] + kmat[, dd[both], drop = FALSE]) / 2 + if (any(mo)) kmat[, ch[mo]] <- kmat[, mm[mo], drop = FALSE] / 2 + if (any(do)) kmat[, ch[do]] <- kmat[, dd[do], drop = FALSE] / 2 + } + for (j in indx) + kmat[j, j] <- if (momIdx[j] != 0L && dadIdx[j] != 0L) (1 + kmat[momIdx[j], dadIdx[j]]) / 2 else 0.5 + } + + kmat +} + # One direction of a set of pairs. Kept separate so the frame and the character image write.table builds from it -- together the # largest transient in the loop below -- are only as long as the pair count, rather than twice that. writePairs <- function(a, b, x, con) @@ -46,17 +202,9 @@ writePairs <- function(a, b, x, con) # Emits both directions of every related pair in one family. Walks the dense matrix in column blocks so the index and value # transients stay bounded no matter how large the family is. -writeFamily <- function(kin, con) +writeFamily <- function(kin, ids, con) { - ids <- colnames(kin) - if (is.null(ids)) - stop('Kinship matrix has no column names') - n <- length(ids) - - # Equivalent to the previous "^(?!addin).*$" filter, but over the ids rather than over every emitted row - keep <- !startsWith(ids, 'addin') - blockCols <- max(1L, as.integer(BLOCK_BYTES %/% (8 * n))) written <- 0 @@ -68,8 +216,7 @@ writeFamily <- function(kin, con) mask <- block != 0 # which() skips an NA in the mask but block[mask] returns an element for it, so a single NA would shift every value - # against its ids from that point on and write wrong coefficients with no error. kinship() should never produce one; - # this turns the case into a visible failure rather than silent corruption. anyNA() scans without allocating. + # against its ids from that point on and write wrong coefficients with no error. anyNA() scans without allocating. if (anyNA(mask)) stop('Kinship matrix contains NA') @@ -90,7 +237,7 @@ writeFamily <- function(kin, con) # Upper triangle only: the matrix is symmetric, so the mirror row is emitted below rather than stored. This also drops # the diagonal, and the importer discards self-pairs anyway. - sel <- i < j & keep[i] & keep[j] + sel <- i < j if (MIN_COEFFICIENT > 0) sel <- sel & vals >= MIN_COEFFICIENT @@ -114,25 +261,28 @@ writeFamily <- function(kin, con) } # Kept in its own frame so the dense matrix is released before the next family allocates its own -processFamily <- function(fixed, idx, con) +processFamily <- function(id, dam, sire, idx, con) { - sub <- fixed[idx, ] - ped <- with(sub, pedigree(id = id, dadid = dadid, momid = momid, sex = sex, missid = 0)) - kin <- kinship(ped) - rm(sub, ped) + ids <- id[idx] + dadIdx <- parentIndex(ids, sire[idx]) + momIdx <- parentIndex(ids, dam[idx]) - writeFamily(kin, con) + kin <- kinshipDense(dadIdx, momIdx, generationDepth(dadIdx, momIdx)) + writeFamily(kin, ids, con) } main <- function() { + # Gender is reported on when it contradicts a parent role; the kinship recursion itself never reads it allPed <- read.table(opts$inputFile, quote="\"") colnames(allPed)<-c('Id', 'Dam', 'Sire', 'Gender', 'Species') + allPed$Id <- as.character(allPed$Id) + allPed$Dam <- as.character(allPed$Dam) + allPed$Sire <- as.character(allPed$Sire) allPed$Id[allPed$Id == ""] <- NA allPed$Dam[allPed$Dam == ""] <- NA allPed$Sire[allPed$Sire == ""] <- NA - allPed$Gender[allPed$Gender == "" | is.na(allPed$Gender)] <- 3 # 3 = unknown allPed$Species <- as.character(allPed$Species) allPed$Species[is.na(allPed$Species)] <- c('Unknown') @@ -146,7 +296,8 @@ main <- function() unlink(TEMP_FILE) con <- file(TEMP_FILE, open = 'wt') - on.exit(if (!is.null(con)) close(con), add = TRUE) + # A successful run has already renamed the temp file away, so this only removes the partial file left by a failure + on.exit({ if (!is.null(con)) close(con); unlink(TEMP_FILE) }, add = TRUE) # The header is retained because the importer detects it by testing whether the third field equals "coefficient" writeLines(paste('Id', 'Id2', 'coefficient', sep = '\t'), con) @@ -161,24 +312,24 @@ main <- function() next } - # Add missing parents for accurate kinship calculations - fixed <- with(allRecordsForSpecies, fixParents(id = Id, dadid = Sire, momid = Dam, sex = Gender)) + validatePedigree(allRecordsForSpecies$Id, allRecordsForSpecies$Dam, allRecordsForSpecies$Sire) + with(allRecordsForSpecies, reportParentSex(species, Id, Dam, Sire, Gender)) + ped <- with(allRecordsForSpecies, closePedigree(Id, Dam, Sire)) rm(allRecordsForSpecies) - # Split into unrelated families. Kinship across families is zero and was never emitted, so the output is unchanged. - fam <- makefamid(fixed$id, fixed$dadid, fixed$momid) - famIds <- unique(fam[fam != 0]) # 0 means no relatives, so there is nothing but a self-pair to emit - print(paste0(' ', length(famIds), ' families, largest ', max(c(0, tabulate(fam[fam != 0]))))) + dadIdx <- parentIndex(ped$id, ped$sire) + momIdx <- parentIndex(ped$id, ped$dam) + fam <- familyIds(dadIdx, momIdx) + sizes <- tabulate(fam) + famIds <- which(sizes >= 2) # a family of one has nothing but a self-pair to emit + print(paste0(' ', length(famIds), ' families, largest ', max(c(0, sizes[famIds])))) speciesRows <- 0 for (f in famIds) { idx <- which(fam == f) - if (length(idx) < 2) - next - totalFamilies <- totalFamilies + 1 - speciesRows <- speciesRows + processFamily(fixed, idx, con) + speciesRows <- speciesRows + processFamily(ped$id, ped$dam, ped$sire, idx, con) # Only worth forcing for families whose matrix is large enough to be worth returning to the OS; # a species can have hundreds of tiny families and R collects those on its own @@ -188,7 +339,7 @@ main <- function() print(paste0(' wrote ', speciesRows, ' rows')) totalRows <- totalRows + speciesRows - rm(fixed, fam) + rm(ped, dadIdx, momIdx, fam) } if (totalRows == 0) diff --git a/ehr/webapp/WEB-INF/ehr/ehrContext.xml b/ehr/webapp/WEB-INF/ehr/ehrContext.xml index 29992da49..fd058163c 100644 --- a/ehr/webapp/WEB-INF/ehr/ehrContext.xml +++ b/ehr/webapp/WEB-INF/ehr/ehrContext.xml @@ -24,7 +24,7 @@ - +