From 45484bc23cea89db68a68ffc7f8d2865bb2585a4 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Tue, 4 Aug 2026 18:31:14 -0400 Subject: [PATCH 01/18] refactor(dataProcess): Stabilize RAM usage for parallel processing --- DESCRIPTION | 4 +- NAMESPACE | 11 + R/MSstatsSummarizeWithMultipleCores.R | 924 +++++++++++++++++++++++ R/dataProcess.R | 93 --- man/MSstatsSummarizeWithMultipleCores.Rd | 53 +- man/dot-MSstatsSummarizeSingleTMPV2.Rd | 48 ++ man/dot-buildProteinSlotV3.Rd | 22 + man/dot-buildSummarizeWorkerV6.Rd | 32 + man/dot-reconstructProteinDTV3.Rd | 21 + 9 files changed, 1107 insertions(+), 101 deletions(-) create mode 100644 R/MSstatsSummarizeWithMultipleCores.R create mode 100644 man/dot-MSstatsSummarizeSingleTMPV2.Rd create mode 100644 man/dot-buildProteinSlotV3.Rd create mode 100644 man/dot-buildSummarizeWorkerV6.Rd create mode 100644 man/dot-reconstructProteinDTV3.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 2bb30c59..0971acfc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -38,7 +38,9 @@ Imports: methods, statmod, parallel, - rlang + rlang, + matter, + BiocParallel Suggests: BiocStyle, knitr, diff --git a/NAMESPACE b/NAMESPACE index e674f5c6..ba6c1354 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,7 @@ export(MSstatsSelectFeatures) export(MSstatsSummarizationOutput) export(MSstatsSummarizeSingleLinear) export(MSstatsSummarizeSingleTMP) +export(MSstatsSummarizeWithMultipleCores) export(MSstatsSummarizeWithSingleCore) export(MZMinetoMSstatsFormat) export(MaxQtoMSstatsFormat) @@ -54,6 +55,12 @@ import(data.table) import(ggplot2) import(limma) import(lme4) +importFrom(BiocParallel,bpisup) +importFrom(BiocParallel,bplapply) +importFrom(BiocParallel,bpnworkers) +importFrom(BiocParallel,bpprogressbar) +importFrom(BiocParallel,bpstart) +importFrom(BiocParallel,bpstop) importFrom(MASS,rlm) importFrom(MSstatsConvert,DIANNtoMSstatsFormat) importFrom(MSstatsConvert,DIAUmpiretoMSstatsFormat) @@ -74,9 +81,12 @@ importFrom(MSstatsConvert,SkylinetoMSstatsFormat) importFrom(MSstatsConvert,SpectronauttoMSstatsFormat) importFrom(Rcpp,sourceCpp) importFrom(data.table,as.data.table) +importFrom(data.table,data.table) +importFrom(data.table,fifelse) importFrom(data.table,melt) importFrom(data.table,rbindlist) importFrom(data.table,setDT) +importFrom(data.table,setDTthreads) importFrom(data.table,uniqueN) importFrom(ggrepel,geom_text_repel) importFrom(gplots,heatmap.2) @@ -97,6 +107,7 @@ importFrom(htmltools,tagList) importFrom(limma,squeezeVar) importFrom(lme4,lmer) importFrom(marray,maPalette) +importFrom(matter,SnowfastParam) importFrom(methods,is) importFrom(parallel,clusterExport) importFrom(parallel,makeCluster) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R new file mode 100644 index 00000000..a36b0cb9 --- /dev/null +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -0,0 +1,924 @@ +## ── Memory-monitoring helpers ───────────────────────────────────────────────── + +# RSS of the current process in MB. +# On Linux reads /proc/self/status (VmRSS); elsewhere falls back to gc() counts. +.memMB <- function() { + if (file.exists("/proc/self/status")) { + ln <- readLines("/proc/self/status", warn = FALSE) + m <- grep("^VmRSS:", ln, value = TRUE) + if (length(m)) + return(as.numeric(gsub("[^0-9]", "", m[1L])) / 1024) + } + g <- gc(reset = FALSE) + (g["Ncells", "used"] * 8L + g["Vcells", "used"] * 8L) / 1024^2 +} + +# Cross-platform peak-RSS reader. Reflects the true lifetime peak of the +# calling process, regardless of when you call it — no polling required. +.peakRSS_MB <- function() { + if (file.exists("/proc/self/status")) { + # Linux: VmHWM = kernel-maintained peak resident set size + ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) + if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024) + } + # macOS (also works as a Linux fallback): POSIX getrusage() ru_maxrss + # is likewise a lifetime high-water mark, just different units per OS. + if (!exists(".rusage_maxrss_mb_impl", mode = "function")) { + Rcpp::cppFunction( + depends = "Rcpp", + includes = "#include ", + code = " + double rusage_maxrss_mb_impl() { + struct rusage ru; getrusage(RUSAGE_SELF, &ru); + #ifdef __APPLE__ + return (double) ru.ru_maxrss / (1024.0*1024.0); // bytes -> MB + #else + return (double) ru.ru_maxrss / 1024.0; // KB -> MB + #endif + }") + assign(".rusage_maxrss_mb_impl", rusage_maxrss_mb_impl, envir = .GlobalEnv) + } + .rusage_maxrss_mb_impl() +} + +# Print a formatted memory report to stderr via message(). +# checkpoints: named numeric vector of RSS snapshots (MB). +# worker_mems: numeric vector of per-worker peak RSS values (may be NA). +# elapsed: total wall-clock seconds (NULL to omit). +.printMemReport <- function(fn_name, checkpoints, worker_mems = NULL, + elapsed = NULL) { + w <- 65L + hr <- strrep("─", w) + fmt_mb <- function(x) if (is.na(x)) " n/a" else sprintf("%7.1f", x) + fmt_delta <- function(now, prev) { + if (is.na(now) || is.na(prev)) return("") + sprintf(" (%+.1f MB)", now - prev) + } + lines <- c( + hr, + sprintf(" MSstats Memory Report — %s", fn_name), + hr, + sprintf(" %-36s %7s %s", "Checkpoint", "RSS MB", "Delta") + ) + prev <- NA_real_ + for (nm in names(checkpoints)) { + val <- checkpoints[[nm]] + lines <- c(lines, + sprintf(" %-36s %s%s", nm, fmt_mb(val), fmt_delta(val, prev))) + prev <- val + } + if (!is.null(worker_mems)) { + wm <- worker_mems[!is.na(worker_mems)] + if (length(wm)) { + lines <- c(lines, "", + sprintf(" Worker RSS min / mean / max : %.1f / %.1f / %.1f MB", + min(wm), mean(wm), max(wm))) + } + } + if (!is.null(elapsed)) + lines <- c(lines, sprintf(" Total elapsed: %.1f s", as.numeric(elapsed))) + lines <- c(lines, hr) + message(paste(lines, collapse = "\n")) +} + +## ── V3 internal helpers ─────────────────────────────────────────────────────── +## +## Packed double-vector layout for one protein slot (all column-major matrices): +## +## pos 1 : slot_k (protein index, cast to double) +## pos 2 .. FL*R+1 : newABUNDANCE (FL × R) +## pos FL*R+2 .. 2*FL*R+1 : ABUNDANCE (FL × R; NA for unlabeled/TMP) +## pos 2*FL*R+2 .. 3*FL*R+1 : censored (FL × R; 0.0/1.0) +## pos 3*FL*R+2 .. 4*FL*R+1 : cen (FL × R; 1-censored event flag) +## pos 4*FL*R+2 .. 5*FL*R+1 : ANOMALYSCORES (FL × R; NA if unused) +## pos 5*FL*R+2 .. 5*FL*R+FL+1 : n_obs (FL; per feature-label) +## pos 5*FL*R+FL+2.. 5*FL*R+FL+R+1: n_obs_run (R; per run) +## pos 5*FL*R+FL+R+2..5*FL*R+FL+2R+1: prop_features (R; per run) +## +## Total length: 1 + 5*FL*R + FL + 2*R +## ───────────────────────────────────────────────────────────────────────────── + +#' Build the packed double vector and lightweight metadata for one protein slot +#' +#' @param dt data.table rows belonging to one protein (or protein × label) slot +#' @param slot_k integer position of this slot in the global protein list +#' @param all_runs character vector of all run names in global order +#' @return list with elements \code{packed} (double vector) and \code{meta} (list) +#' @keywords internal +.buildProteinSlotV3 <- function(dt, slot_k, all_runs) { + + R <- length(all_runs) + + # ── Unique (FEATURE, LABEL) → "effective features", stable ordering ─────── + has_peptide <- "PEPTIDE" %in% colnames(dt) + flp <- unique(dt[, .( + FEATURE = as.character(FEATURE), + LABEL = as.character(LABEL), + PEPTIDE = if (has_peptide) as.character(PEPTIDE) else as.character(FEATURE) + )]) + data.table::setorder(flp, FEATURE, LABEL) + FL <- nrow(flp) + + # ── Index maps: feature-label → row index, run name → column index ───────── + fi <- match( + paste(as.character(dt$FEATURE), as.character(dt$LABEL), sep = "\t"), + paste(flp$FEATURE, flp$LABEL, sep = "\t")) + ri_map <- seq_len(R); names(ri_map) <- all_runs + ri <- ri_map[as.character(dt$RUN)] + valid <- !is.na(fi) & !is.na(ri) + + # ── Helper: allocate FL × R matrix and fill from dt column ─────────────── + fill_mat <- function(col_vec, fill = NA_real_) { + m <- matrix(fill, nrow = FL, ncol = R) + m[cbind(fi[valid], ri[valid])] <- as.double(col_vec[valid]) + m + } + + # ── Build the five numeric matrices ─────────────────────────────────────── + mat_newABU <- fill_mat(dt$newABUNDANCE) + + has_ABUNDANCE <- "ABUNDANCE" %in% colnames(dt) + mat_ABU <- if (has_ABUNDANCE) fill_mat(dt$ABUNDANCE) else + matrix(NA_real_, FL, R) + + has_censored <- "censored" %in% colnames(dt) + mat_cens <- if (has_censored) fill_mat(as.double(dt$censored)) else + matrix(0.0, FL, R) # treat as non-censored when column absent + + has_cen <- "cen" %in% colnames(dt) + mat_cen <- if (has_cen) fill_mat(dt$cen) else + matrix(NA_real_, FL, R) + + has_anom <- "ANOMALYSCORES" %in% colnames(dt) && + !all(is.na(dt$ANOMALYSCORES)) + mat_anom <- if (has_anom) fill_mat(dt$ANOMALYSCORES) else + matrix(NA_real_, FL, R) + + # ── Per-feature scalar: n_obs (constant within FEATURE × LABEL) ─────────── + n_obs_by_fl <- dt[, .(n_obs = as.double(n_obs[1L])), + by = .(FEATURE = as.character(FEATURE), + LABEL = as.character(LABEL))] + flp_nobs <- n_obs_by_fl[flp, on = c("FEATURE", "LABEL")] + data.table::setorder(flp_nobs, FEATURE, LABEL) + n_obs_vec <- flp_nobs$n_obs + + # ── Per-run scalars: n_obs_run, prop_features (constant within RUN) ─────── + run_vals <- dt[, .(n_obs_run = as.double(n_obs_run[1L]), + prop_features = as.double(prop_features[1L])), + by = .(RUN = as.character(RUN))] + run_full <- run_vals[data.table::data.table(RUN = all_runs), on = "RUN"] + n_obs_run_vec <- run_full$n_obs_run + prop_features_vec <- run_full$prop_features + + # ── Pack ───────────────────────────────────────────────────────────────── + packed <- c( + as.double(slot_k), + as.vector(mat_newABU), + as.vector(mat_ABU), + as.vector(mat_cens), + as.vector(mat_cen), + as.vector(mat_anom), + n_obs_vec, + n_obs_run_vec, + prop_features_vec + ) + + # ── Metadata (string labels, flags; kept in main-process RAM) ──────────── + meta <- list( + PROTEIN = as.character(dt$PROTEIN[1L]), + feat_label_pep = as.data.frame(flp), # FL × 3: FEATURE, LABEL, PEPTIDE + runs = all_runs, + FL = FL, + R = R, + is_labeled_ref = "is_labeled_ref" %in% colnames(dt) && + isTRUE(any(dt$is_labeled_ref, na.rm = TRUE)), + has_ABUNDANCE = has_ABUNDANCE, + has_censored = has_censored, + has_cen = has_cen, + has_anom = has_anom, + add_ref_covariate = "ref_covariate" %in% colnames(dt) + ) + + list(packed = packed, meta = meta) +} + + +#' Reconstruct a per-protein data.table from a V3 packed double vector +#' +#' @param packed double vector produced by \code{.buildProteinSlotV3} +#' @param meta metadata list from \code{.buildProteinSlotV3} +#' @return data.table compatible with \code{MSstatsSummarizeSingleTMP} / +#' \code{MSstatsSummarizeSingleLinear} +#' @keywords internal +.reconstructProteinDTV3 <- function(packed, meta) { + + FL <- meta$FL + R <- meta$R + n_mat <- FL * R + + # ── Unpack sections (1-indexed; position 1 is slot_k header) ───────────── + i <- 2L + extract_mat <- function() { + m <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R) + i <<- i + n_mat + m + } + mat_newABU <- extract_mat() + mat_ABU <- extract_mat() + mat_cens <- extract_mat() + mat_cen <- extract_mat() + mat_anom <- extract_mat() + + n_obs_vec <- packed[i:(i + FL - 1L)]; i <- i + FL + n_obs_run_vec <- packed[i:(i + R - 1L)]; i <- i + R + prop_feat_vec <- packed[i:(i + R - 1L)] + + # ── Column-major melt: rep/each match the matrix column-major ordering ──── + flp <- meta$feat_label_pep # data.frame: FEATURE, LABEL, PEPTIDE (FL rows) + runs <- meta$runs + n_rows <- FL * R + + dt <- data.table::data.table( + PROTEIN = rep(meta$PROTEIN, n_rows), + FEATURE = rep(flp$FEATURE, times = R), + LABEL = rep(flp$LABEL, times = R), + PEPTIDE = rep(flp$PEPTIDE, times = R), + RUN = rep(runs, each = FL), + newABUNDANCE = as.vector(mat_newABU), + n_obs = as.integer(rep(n_obs_vec, times = R)), + n_obs_run = as.integer(rep(n_obs_run_vec, each = FL)), + prop_features = rep(prop_feat_vec, each = FL) + ) + + # Optional: ABUNDANCE (labeled linear model uses raw log-intensities) + if (meta$has_ABUNDANCE) { + dt[, ABUNDANCE := as.vector(mat_ABU)] + } + + # censored: stored as 0.0/1.0 doubles; NA → treat as non-censored + dt[, censored := { + v <- as.vector(mat_cens) + if (meta$has_censored) as.logical(v > 0.5) else rep(FALSE, n_rows) + }] + + # cen: survival event indicator (1 = observed, 0 = left-censored) + if (meta$has_cen) { + dt[, cen := as.vector(mat_cen)] + } + + # ANOMALYSCORES: NA when unused (linear model skips anomaly weighting) + dt[, ANOMALYSCORES := as.vector(mat_anom)] + + # SRM-specific columns derived from LABEL and RUN + if (meta$is_labeled_ref) { + dt[, is_labeled_ref := (LABEL == "H")] + if (meta$add_ref_covariate) { + dt[, ref_covariate := factor( + data.table::fifelse(LABEL == "L", as.character(RUN), "0"))] + } + } + + dt +} + + + + + +## ── V4: matrix-native TMP ───────────────────────────────────────────────────── +## +## .MSstatsSummarizeSingleTMPV2 takes the V3 packed double vector + meta list +## directly and produces TMP results without ever building the full long-format +## data.table or calling dcast. +## +## Compared with MSstatsSummarizeSingleTMP the hot path is: +## V3: packed → long-DT (melt) → filter → AFT fit → dcast → TMP +## V4: packed → FL×R matrix → mask → (AFT fit if needed) → t(mat) → TMP +## +## The dcast savings are most significant for high-feature proteins (>100 +## features) and large run counts (>50 runs), where the intermediate FL*R +## data.table allocation dominates. +## ───────────────────────────────────────────────────────────────────────────── + + +#' Summarize a single protein with TMP directly from a V3 packed double vector +#' +#' Bypasses the long-format \code{data.table} reconstruction and the +#' \code{dcast} inside \code{.fitTukey}. The FL×R packed matrices are +#' operated on directly: +#' \enumerate{ +#' \item Row/column masking replaces the \code{n_obs}/\code{n_obs_run} filter. +#' \item AFT survival fitting uses a lazily-built minimal \code{data.table} +#' (constructed only when \code{impute=TRUE} and censored values are +#' present). +#' \item TMP is applied via a vectorized scatter into a +#' \code{(LABEL×RUN) × FEATURE} wide matrix followed by +#' \code{median_polish_summary} — no \code{dcast} round-trip. +#' } +#' +#' @param packed double vector produced by \code{.buildProteinSlotV3} +#' @param meta metadata list from \code{.buildProteinSlotV3} +#' @param impute logical; impute censored values with AFT survival model +#' @param censored_symbol \code{"0"}, \code{"NA"}, or \code{NULL} +#' @param remove50missing logical; skip proteins where all runs are >50\% +#' missing +#' @param aft_iterations integer; max AFT iterations +#' @return \code{list(result_dt, survival_dt)} matching the format of +#' \code{MSstatsSummarizeSingleTMP} +#' @keywords internal +.MSstatsSummarizeSingleTMPV2 <- function(packed, meta, + impute, censored_symbol, remove50missing, aft_iterations = 90L) +{ + FL <- meta$FL + R <- meta$R + n_mat <- FL * R + PROTEIN <- meta$PROTEIN + flp <- meta$feat_label_pep # data.frame: FEATURE, LABEL, PEPTIDE + runs <- meta$runs + + # ── 1. Unpack matrices ───────────────────────────────────────────────────── + # Layout (see top-of-file comment block): + # pos 1 : slot_k (skip) + # pos 2..FL*R+1 : newABUNDANCE (FL × R, column-major) + # pos FL*R+2.. : ABUNDANCE (skip — not needed for TMP) + # pos 2*FL*R+2.. : censored (FL × R) + # pos 3*FL*R+2.. : cen (FL × R) + # pos 4*FL*R+2.. : ANOMALYSCORES(skip — not needed for TMP) + # then: n_obs(FL), n_obs_run(R), prop_features(R) + i <- 2L + mat_newABU <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R); i <- i + n_mat + i <- i + n_mat # skip ABUNDANCE + mat_cens <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R); i <- i + n_mat + mat_cen <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R); i <- i + n_mat + i <- i + n_mat # skip ANOMALYSCORES + n_obs_vec <- packed[i:(i + FL - 1L)]; i <- i + FL + n_obs_run_vec <- packed[i:(i + R - 1L)]; i <- i + R + prop_feat_vec <- packed[i:(i + R - 1L)] + + # ── 2. Row/column masking (mirrors the n_obs / n_obs_run filter) ─────────── + keep_feat <- !is.na(n_obs_vec) & n_obs_vec > 1 + keep_runs <- !is.na(n_obs_run_vec) & n_obs_run_vec > 0 + + mat_newABU <- mat_newABU[keep_feat, keep_runs, drop = FALSE] + mat_cens <- mat_cens [keep_feat, keep_runs, drop = FALSE] + mat_cen <- mat_cen [keep_feat, keep_runs, drop = FALSE] + + flp_filt <- flp[keep_feat, , drop = FALSE] + runs_filt <- runs[keep_runs] + prop_feat_filt <- prop_feat_vec[keep_runs] + FL_filt <- nrow(mat_newABU) + R_filt <- ncol(mat_newABU) + + # ── 3. Empty-matrix early exit ───────────────────────────────────────────── + if (FL_filt == 0L || R_filt == 0L) { + msg <- paste("Can't summarize for protein", PROTEIN, + "because all measurements are missing or censored.") + try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) + try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) + return(list(NULL, NULL)) + } + + # ── 4. remove50missing check (independent of imputation) ────────────────── + if (remove50missing && + all(prop_feat_filt <= 0.5 | is.na(prop_feat_filt))) { + msg <- paste("Can't summarize for protein", PROTEIN, + "because all runs have more than 50% missing values and", + "are removed with the option, remove50missing=TRUE.") + try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) + try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) + return(list(NULL, NULL)) + } + + # ── 5. AFT survival imputation (lazy: only if any censored present) ──────── + n_rows_filt <- FL_filt * R_filt + orig_abu_vec <- as.vector(mat_newABU) + + surv_cols <- c("newABUNDANCE", "cen", "RUN", "FEATURE", "ref_covariate") + any_censored <- meta$has_censored && any(mat_cens > 0.5, na.rm = TRUE) + + if (impute && any_censored && meta$has_cen) { + # Build minimal long-format DT — only columns needed by .fitSurvival + surv_dt <- data.table::data.table( + newABUNDANCE = orig_abu_vec, + cen = as.vector(mat_cen), + censored = as.logical(as.vector(mat_cens) > 0.5), + FEATURE = factor(rep(flp_filt$FEATURE, times = R_filt)), + RUN = factor(rep(runs_filt, each = FL_filt)), + LABEL = rep(as.character(flp_filt$LABEL), times = R_filt) + ) + if (meta$add_ref_covariate) { + surv_dt[, ref_covariate := factor( + data.table::fifelse(LABEL == "L", as.character(RUN), "0"))] + } + + fit_cols <- intersect(surv_cols, colnames(surv_dt)) + fit_data <- if (meta$is_labeled_ref) { + surv_dt[LABEL != "H", fit_cols, with = FALSE] + } else { + surv_dt[, fit_cols, with = FALSE] + } + + converged <- TRUE + survival_fit <- withCallingHandlers({ + MSstats:::.fitSurvival(fit_data, aft_iterations) + }, warning = function(w) { + if (grepl("converge", conditionMessage(w), ignore.case = TRUE)) { + message("Convergence warning caught: ", conditionMessage(w)) + converged <<- FALSE + } + }) + + predicted_all <- if (converged) { + predict(survival_fit, newdata = surv_dt) + } else { + rep(NA_real_, n_rows_filt) + } + + use_imputed <- if (meta$is_labeled_ref) { + surv_dt$censored & surv_dt$LABEL != "H" + } else { + surv_dt$censored + } + + imputed_abu <- ifelse(use_imputed, predicted_all, orig_abu_vec) + surv_dt[, predicted := ifelse(use_imputed, predicted_all, NA_real_)] + surv_dt[, newABUNDANCE := imputed_abu] + mat_newABU <- matrix(imputed_abu, nrow = FL_filt, ncol = R_filt) + survival <- surv_dt[, intersect(c(surv_cols, "LABEL", "predicted"), + colnames(surv_dt)), with = FALSE] + + } else { + # No imputation: build survival DT from original values + surv_dt <- data.table::data.table( + newABUNDANCE = orig_abu_vec, + cen = if (meta$has_cen) as.vector(mat_cen) + else rep(NA_real_, n_rows_filt), + FEATURE = rep(as.character(flp_filt$FEATURE), times = R_filt), + RUN = rep(runs_filt, each = FL_filt), + LABEL = rep(as.character(flp_filt$LABEL), times = R_filt), + predicted = NA # logical NA, matching MSstatsSummarizeSingleTMP's + # bare `survival[, predicted := NA]` in its + # equivalent non-imputed branch (dataProcess.R:586) + ) + survival <- surv_dt[, intersect(c(surv_cols, "LABEL", "predicted"), + colnames(surv_dt)), with = FALSE] + } + + # ── 6. Post-imputation .isSummarizable check ─────────────────────────────── + if (all(is.na(mat_newABU) | mat_newABU == 0)) { + msg <- paste("Can't summarize for protein", PROTEIN, + "because all measurements are missing or censored.") + try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) + try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) + return(list(NULL, NULL)) + } + + # ── 7. TMP via vectorized scatter into (LABEL×RUN) × FEATURE wide matrix ─── + # + # Mirrors .fitTukey: dcast(LABEL + RUN ~ FEATURE, value.var="newABUNDANCE") + # Row order in wide_mat: (lbl1,run1),(lbl1,run2),...,(lbl2,run1),... + # where labels are sorted — matches dcast alphabetical ordering for + # character LABEL values. + if (FL_filt == 1L) { + # ── Single (FEATURE, LABEL) row: skip TMP, use values directly ──────── + if (meta$is_labeled_ref) { + h_sel <- as.character(flp_filt$LABEL) == "H" + l_sel <- as.character(flp_filt$LABEL) == "L" + if (any(h_sel) && any(l_sel)) { + h_vals <- as.vector(mat_newABU[h_sel, , drop = FALSE]) + l_vals <- as.vector(mat_newABU[l_sel, , drop = FALSE]) + h_median <- stats::median(h_vals, na.rm = TRUE) + result <- data.table::data.table( + Protein = PROTEIN, + LABEL = "L", + RUN = runs_filt, + LogIntensities = l_vals - h_vals + h_median + ) + } else { + result <- data.table::data.table( + Protein = PROTEIN, + LABEL = as.character(flp_filt$LABEL[1L]), + RUN = runs_filt, + LogIntensities = as.vector(mat_newABU) + ) + } + } else { + result <- data.table::data.table( + Protein = PROTEIN, + LABEL = rep(as.character(flp_filt$LABEL), times = R_filt), + RUN = rep(runs_filt, each = FL_filt), + LogIntensities = as.vector(mat_newABU) + ) + } + } else { + # ── Multi-feature: scatter into wide matrix, apply TMP ───────────────── + unique_labels <- sort(unique(as.character(flp_filt$LABEL))) + unique_feats <- sort(unique(as.character(flp_filt$FEATURE))) + n_lbl <- length(unique_labels) + n_feat <- length(unique_feats) + + lbl_idx_map <- match(as.character(flp_filt$LABEL), unique_labels) + feat_idx_map <- match(as.character(flp_filt$FEATURE), unique_feats) + + # Vectorized scatter: mat[fl, run] → wide[(lbl_idx-1)*R + run, feat_idx] + fl_rep <- rep(seq_len(FL_filt), times = R_filt) + run_rep <- rep(seq_len(R_filt), each = FL_filt) + wide_row <- (lbl_idx_map[fl_rep] - 1L) * R_filt + run_rep + wide_col <- feat_idx_map[fl_rep] + + wide_mat <- matrix(NA_real_, nrow = n_lbl * R_filt, ncol = n_feat) + wide_mat[cbind(wide_row, wide_col)] <- as.vector(mat_newABU) + + tmp_fitted <- MSstats:::median_polish_summary(wide_mat) + + # Row → (LABEL, RUN) mapping: label index cycles every R_filt rows + result_labels <- rep(unique_labels, each = R_filt) + result_runs <- rep(runs_filt, times = n_lbl) + + if (meta$is_labeled_ref) { + h_li <- match("H", unique_labels) + l_li <- match("L", unique_labels) + if (!is.na(h_li) && !is.na(l_li)) { + h_rows <- (h_li - 1L) * R_filt + seq_len(R_filt) + l_rows <- (l_li - 1L) * R_filt + seq_len(R_filt) + h_med <- stats::median(tmp_fitted[h_rows], na.rm = TRUE) + result <- data.table::data.table( + Protein = PROTEIN, + LABEL = "L", + RUN = runs_filt, + LogIntensities = tmp_fitted[l_rows] - tmp_fitted[h_rows] + h_med + ) + } else { + result <- data.table::data.table( + Protein = PROTEIN, + LABEL = result_labels, + RUN = result_runs, + LogIntensities = tmp_fitted + ) + } + } else { + result <- data.table::data.table( + Protein = PROTEIN, + LABEL = result_labels, + RUN = result_runs, + LogIntensities = tmp_fitted + ) + } + } + + list(result, survival) +} + +#' Build the per-record worker closure for +#' \code{MSstatsSummarizeWithMultipleCoresV6} +#' +#' Defined at package top level — not nested via \code{local()} inside +#' \code{MSstatsSummarizeWithMultipleCoresV6} — so the returned closure's +#' enclosing environment chain is this factory's own (small) call frame plus +#' the package namespace. \code{BiocParallel}/\code{matter} serialize a +#' closure's entire enclosing environment chain to ship it to each socket +#' worker, not just the variables the closure body actually references. A +#' closure built via \code{local()} inside +#' \code{MSstatsSummarizeWithMultipleCoresV6} would have that function's own +#' evaluation frame in its chain — which holds \code{input}, +#' \code{protein_records}, and other run-scale objects — so every worker +#' would receive and retain a serialized copy of them even though +#' \code{.worker} never touches them. Building the closure here instead keeps +#' its captured state to only the scalar run parameters. +#' +#' @keywords internal +.buildSummarizeWorkerV6 <- function( + use_TMP, impute, censored_symbol, remove50missing, + aft_iterations, equal_variance +) { + reconstruct_ <- .reconstructProteinDTV3 + use_TMP_ <- use_TMP + impute_ <- impute + censored_symbol_ <- censored_symbol + remove50missing_ <- remove50missing + aft_iterations_ <- aft_iterations + equal_variance_ <- equal_variance + + function(record) { + packed <- record$packed + meta <- record$meta + result <- if (use_TMP_) { + .MSstatsSummarizeSingleTMPV2( + packed, meta, + impute_, censored_symbol_, + remove50missing_, aft_iterations_) + } else { + dt <- reconstruct_(packed, meta) + MSstatsSummarizeSingleLinear( + dt, impute_, censored_symbol_, + remove50missing_, aft_iterations_, + equal_variances = equal_variance_) + } + # Normalize column types/levels at the V6 worker boundary rather than + # inside .MSstatsSummarizeSingleTMPV2/MSstatsSummarizeSingleLinear: + # those are shared with V3-V5 and dataProcess.R respectively, and the + # TMP path in particular carries RUN as plain character (from + # meta$runs) and cen as double (from the packed matrix). Re-leveling + # against meta$runs (rather than a bare factor(RUN)) preserves the + # numeric run order that levels(input$RUN) already established + # upstream — a bare factor() call would instead re-sort alphabetically + # ("1","10","11","2",... for >=10 runs). as.character() first strips + # any existing level order before re-leveling, since .MSstatsSummarize + # SingleTMPV2's survival table (result[[2L]]) already wraps RUN in a + # bare factor() for the imputed branch — that call has the exact same + # alphabetical-resort bug even though it looks pre-typed. droplevels() + # then matches MSstatsSummarizeSingleTMP/SingleCore, whose factor(RUN) + # only ever sees — and so only ever keeps — the runs present for that + # protein, since meta$runs carries every run across the whole input. + # + # FEATURE gets the same treatment: .MSstatsSummarizeSingleTMPV2's + # non-imputed survival branch leaves FEATURE as plain character (only + # the imputed branch wraps it in a bare factor()), so result[[2L]]$ + # FEATURE isn't reliably a factor the way SingleCore's is (single_ + # protein[, FEATURE := factor(FEATURE)] before survival is sliced off + # it). + # + # Deliberately NOT leveling against meta$feat_label_pep$FEATURE's + # existing row order here: that order comes from .buildProteinSlotV3's + # data.table::setorder(flp, FEATURE, LABEL), and data.table sorts + # character columns in the C-locale (byte order: all uppercase before + # any lowercase) for platform-independence — whereas SingleCore's + # bare factor(FEATURE) goes through base R's factor()/sort(), which + # use the session's collation locale (e.g. en_US.UTF-8, where case is + # interleaved: "a" < "A" < "b" < "B"). Any FEATURE strings with mixed + # case (e.g. modification tags) sort differently between the two, so + # re-deriving the level order with a base R sort() reproduces + # SingleCore's order instead of inheriting data.table's. + feature_levels <- sort(unique(meta$feat_label_pep$FEATURE)) + for (idx in 1:2) { + if (!is.null(result[[idx]])) { + if ("RUN" %in% colnames(result[[idx]])) + result[[idx]][, RUN := droplevels(factor(as.character(RUN), levels = meta$runs))] + if ("FEATURE" %in% colnames(result[[idx]])) + result[[idx]][, FEATURE := droplevels(factor(as.character(FEATURE), levels = feature_levels))] + } + } + if (!is.null(result[[2L]]) && "cen" %in% colnames(result[[2L]])) + result[[2L]][, cen := as.integer(cen)] + result + } +} + +#' Per-worker peak-RSS query task for \code{MSstatsSummarizeWithMultipleCoresV6} +#' +#' Dispatched once per worker (via \code{seq_len(bpnworkers(BPPARAM))}) after +#' the main summarization \code{bplapply} call, while the persistent workers +#' are still alive, so each worker reports its own true lifetime-peak RSS +#' rather than a snapshot taken mid-run. +#' +#' @keywords internal +#' @noRd +.reportWorkerPeakV6 <- function(i) { + list(worker = i, pid = Sys.getpid(), peak_mb = .peakRSS_MB()) +} + +# Per-worker warm-up task for MSstatsSummarizeWithMultipleCoresV6: loads +# MSstats once per persistent worker process and pins data.table to a single +# thread. Defined at top level rather than inline inside +# MSstatsSummarizeWithMultipleCoresV6: even though its body never references +# `input`/`protein_records`, a closure defined inline there would still carry +# that function's evaluation frame in its enclosing environment chain, and +# BiocParallel would serialize that whole frame — including those run-scale +# objects — to every worker just to ship this no-op task. +.warmupV6Worker <- function(i) { + library(MSstats, quietly = TRUE, warn.conflicts = FALSE) + data.table::setDTthreads(1) + NULL +} + + +#' Feature-level data summarization via socket-dispatched protein records (V6) +#' +#' Fixes a hidden RAM cost in \code{MSstatsSummarizeWithMultipleCoresV5}: V5's +#' worker closure captured \code{meta_list} — the metadata for \emph{every} +#' protein — by reference, so each worker received and retained the full +#' metadata set for the whole run regardless of how many proteins were +#' actually assigned to it. Bounding the per-task \emph{packed-vector} payload +#' (e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce that, +#' since the metadata was baked into the closure once, not re-sliced per task. +#' +#' V6 pairs each protein's packed double vector with only its own metadata +#' into a single list element (\code{list(packed = ..., meta = ...)}) and +#' dispatches that combined list through \code{bplapply}. \code{BiocParallel} +#' then only serializes and sends the records actually assigned to a given +#' task, so a worker never holds metadata for proteins outside its own +#' task(s) — unlike V5, RAM scales with the batch actually being processed. +#' +#' Progress is reported by turning on \code{SnowfastParam}'s built-in +#' \code{progressbar}, not by having workers talk back to the parent. +#' \code{BiocParallel} already ticks that progress bar from inside the +#' manager process itself, once per task result it collects — the same +#' receive that would happen regardless, over the same single +#' \code{bplapply} call. Workers never see the flag and send nothing extra +#' because of it, so this adds no serialization and no IPC beyond what an +#' unmonitored run already does. Reporting granularity therefore tracks +#' however many tasks the run is already split into (see +#' \code{max_proteins_per_worker} below): with the default \code{tasks = 0} +#' that's one step per worker. +#' +#' @param input feature-level data processed by dataProcess subfunctions +#' @param method summarization method: "linear" or "TMP" +#' @param equal_variance only for summaryMethod = "linear". Default is TRUE. +#' Logical variable for whether the model should account for heterogeneous variation +#' among intensities from different features. Default is TRUE, which assume equal +#' variance among intensities from features. FALSE means that we cannot assume +#' equal variance among intensities from features, then we will account for +#' heterogeneous variation from different features. +#' @param censored_symbol Indicates how censored missing values are encoded in +#' the 'Intensity' column. 'NA' (default) treats all NA intensities as +#' left-censored (i.e., below the limit of detection). '0' treats zero +#' intensities as left-censored; in this case NA intensities are assumed to be +#' missing at random and are not censored. Skyline output should use '0'. NULL +#' assumes that all missing values are missing at random — no values are treated +#' as censored, and imputation is disabled (impute is ignored). +#' @param remove50missing only for summaryMethod = "TMP". TRUE removes the proteins +#' where every run has at least 50\% missing values for each peptide. FALSE is default. +#' @param impute only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'. +#' TRUE (default) imputes censored missing values using an Accelerated Failure +#' Time model. FALSE excludes censored observations from summarization entirely, +#' treating them as missing at random; no imputed values are introduced. +#' Has no effect when censored_symbol = NULL, since no values are considered censored. +#' @param numberOfCores Number of cores for parallel processing. When > 1, +#' a logfile named `MSstats_dataProcess_log_progress.log` is created to +#' track progress. Only works for Linux & Mac OS. Default is 1. +#' @param max_proteins_per_worker integer; caps how many protein records +#' (packed vector + its own metadata) are bundled into a single +#' \code{bplapply} task sent to one worker. Translated into +#' \code{SnowfastParam(tasks = ...)}: with \code{N} proteins this becomes +#' \code{tasks = ceiling(N / max_proteins_per_worker)}. Only applied when +#' \code{BPPARAM} is \code{NULL}; ignored if the caller supplies +#' \code{BPPARAM} directly. Default \code{0} reproduces \code{tasks = 0}: +#' \code{X} is divided as evenly as possible across \code{numberOfCores} +#' workers. +#' +#' @return A named list with one element per protein slot, identical in +#' structure to \code{MSstatsSummarizeWithMultipleCores}. +#' +#' @importFrom matter SnowfastParam +#' @importFrom BiocParallel bplapply bpstart bpstop bpisup bpnworkers bpprogressbar +#' @importFrom data.table data.table fifelse setDTthreads +#' @importFrom stats median predict +#' +#' @export +MSstatsSummarizeWithMultipleCores <- function( + input, + method, + impute, + censored_symbol, + remove50missing, + equal_variance, + numberOfCores = 1L, + aft_iterations = 90L, + verbose = FALSE, + BPPARAM = NULL, + track_memory = FALSE, + max_proteins_per_worker = 0L +) { + # ── 0. Single-core fallback ──────────────────────────────────────────────── + if (numberOfCores <= 1L && is.null(BPPARAM)) { + return(MSstatsSummarizeWithSingleCore( + input, method, impute, censored_symbol, + remove50missing, equal_variance, aft_iterations)) + } + + t_start <- proc.time()[["elapsed"]] + mem_log <- list() + + # ── 1. Split input by protein slot ──────────────────────────────────────── + is_labeled_reference <- "is_labeled_ref" %in% colnames(input) && + any(input$is_labeled_ref, na.rm = TRUE) + split_keys <- if (is_labeled_reference) list(input$PROTEIN) else + list(input$PROTEIN, input$LABEL) + protein_indices <- split(seq_len(nrow(input)), split_keys) + protein_ids <- names(protein_indices) + num_proteins <- length(protein_indices) + + # Sort on input$RUN's own type (integer/numeric, typically) before + # converting to character. Sorting the character form directly — as a + # naive `sort(as.character(...))` would — collates lexicographically + # ("1","10","11","2",...) whenever RUN isn't already a pre-existing + # factor, e.g. when `input` comes straight from fread() and RUN reads in + # as integer. This mirrors what factor() itself does internally for a + # non-factor input: unique() + order() on the original values, then + # as.character() only at the end. + all_runs <- if (is.factor(input$RUN)) levels(input$RUN) else + as.character(sort(unique(input$RUN))) + + getOption("MSstatsLog")("INFO", + paste0("V6: packing ", num_proteins, " proteins × ", + length(all_runs), " runs into per-protein records")) + + # ── 2. Pack each protein into a (packed vector, own metadata) record ────── + # + # Unlike V5's parallel packed_list/meta_list arrays, each protein's + # metadata travels bundled with its own packed vector. bplapply/BPPARAM + # slice this combined list into tasks, so a worker's closure never needs + # — and never receives — metadata for proteins outside its assigned + # task(s). + protein_records <- vector("list", num_proteins) + for (k in seq_len(num_proteins)) { + slot <- .buildProteinSlotV3( + input[protein_indices[[k]], ], k, all_runs) + protein_records[[k]] <- list(packed = slot$packed, meta = slot$meta) + } + + payload_mb <- sum(vapply(protein_records, + function(r) length(r$packed), integer(1))) * 8 / 1024^2 + getOption("MSstatsLog")("INFO", + paste0("V6: dispatching via sockets (", + format(round(payload_mb, 1)), + " MB total packed payload; metadata sharded per task)")) + + # ── 3. Worker closure ───────────────────────────────────────────────────── + # + # No captured meta_list: each task's records already carry their own + # metadata, so the closure only needs the scalar run parameters. + # + # Built via .buildSummarizeWorkerV6() (defined at package top level) + # rather than local() here, so the closure's enclosing environment chain + # never includes this function's own frame — which holds `input`, + # `protein_records`, etc. — and BiocParallel doesn't serialize those + # run-scale objects to every worker. + use_TMP <- identical(method, "TMP") + + .worker <- .buildSummarizeWorkerV6( + use_TMP, impute, censored_symbol, remove50missing, + aft_iterations, equal_variance) + + # ── 4. Dispatch ─────────────────────────────────────────────────────────── + # + # tasks controls how many protein records are bundled into one bplapply + # task (i.e. one message sent to one worker). tasks == 0 (default) leaves + # BiocParallel's own behavior in place: X divided as evenly as possible + # over numberOfCores workers. When max_proteins_per_worker > 0, tasks is + # sized so no task exceeds that many records — and because metadata now + # travels with each record instead of being fully captured by the + # closure, this actually bounds peak worker RAM. + if (is.null(BPPARAM)) { + tasks <- if (max_proteins_per_worker > 0L) { + as.integer(ceiling(num_proteins / max_proteins_per_worker)) + } else { + 0L + } + getOption("MSstatsLog")("INFO", + paste0("V6: dispatching as ", + if (tasks > 0L) tasks else numberOfCores, + " task(s) (max_proteins_per_worker = ", + max_proteins_per_worker, ")")) + BPPARAM <- matter::SnowfastParam( + workers = numberOfCores, + tasks = tasks, + progressbar = TRUE, + force.GC = TRUE, + stop.on.error = FALSE) + } + + # ── Cluster setup ───────────────────────────────────────────────────────── + # Load MSstats once per persistent worker process instead of once per + # protein record, and pin data.table to a single thread per worker — + # otherwise each of the numberOfCores workers would independently + # auto-detect its own DT thread pool, oversubscribing the node's cores. + started_here <- !BiocParallel::bpisup(BPPARAM) + if (started_here) { + BiocParallel::bpstart(BPPARAM) + on.exit(BiocParallel::bpstop(BPPARAM), add = TRUE) + } + # Progress bar is on BPPARAM itself, so it would otherwise also print for + # these two trivial one-task-per-worker calls. Toggle it off around them + # and restore whatever it was (set above) so only the main summarization + # bplapply — the one call whose progress is actually informative — shows + # a bar. + show_progress <- BiocParallel::bpprogressbar(BPPARAM) + + BiocParallel::bpprogressbar(BPPARAM) <- FALSE + BiocParallel::bplapply( + seq_len(BiocParallel::bpnworkers(BPPARAM)), + .warmupV6Worker, BPPARAM = BPPARAM) + BiocParallel::bpprogressbar(BPPARAM) <- show_progress + + results <- BiocParallel::bplapply(protein_records, .worker, BPPARAM = BPPARAM) + names(results) <- protein_ids + + worker_peaks <- NULL + if (track_memory) { + BiocParallel::bpprogressbar(BPPARAM) <- FALSE + worker_peaks <- BiocParallel::bplapply( + seq_len(BiocParallel::bpnworkers(BPPARAM)), + .reportWorkerPeakV6, BPPARAM = BPPARAM) # must run BEFORE bpstop() while workers are alive + BiocParallel::bpprogressbar(BPPARAM) <- show_progress + mem_log[["parent peak (main)"]] <- .peakRSS_MB() + worker_peak_mb <- vapply(worker_peaks, function(x) x$peak_mb, numeric(1L)) + .printMemReport( + "MSstatsSummarizeWithMultipleCoresV6", + mem_log, worker_peak_mb, + elapsed = proc.time()[["elapsed"]] - t_start) + } + + getOption("MSstatsLog")("INFO", "V6: summarization complete.") + results +} diff --git a/R/dataProcess.R b/R/dataProcess.R index 073d72b0..ef23f518 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -181,99 +181,6 @@ dataProcess = function( output } -#' Feature-level data summarization with multiple cores -#' -#' @param input feature-level data processed by dataProcess subfunctions -#' @param method summarization method: "linear" or "TMP" -#' @param equal_variance only for summaryMethod = "linear". Default is TRUE. -#' Logical variable for whether the model should account for heterogeneous variation -#' among intensities from different features. Default is TRUE, which assume equal -#' variance among intensities from features. FALSE means that we cannot assume -#' equal variance among intensities from features, then we will account for -#' heterogeneous variation from different features. -#' @param censored_symbol Indicates how censored missing values are encoded in -#' the 'Intensity' column. 'NA' (default) treats all NA intensities as -#' left-censored (i.e., below the limit of detection). '0' treats zero -#' intensities as left-censored; in this case NA intensities are assumed to be -#' missing at random and are not censored. Skyline output should use '0'. NULL -#' assumes that all missing values are missing at random — no values are treated -#' as censored, and imputation is disabled (impute is ignored). -#' @param remove50missing only for summaryMethod = "TMP". TRUE removes the proteins -#' where every run has at least 50\% missing values for each peptide. FALSE is default. -#' @param impute only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'. -#' TRUE (default) imputes censored missing values using an Accelerated Failure -#' Time model. FALSE excludes censored observations from summarization entirely, -#' treating them as missing at random; no imputed values are introduced. -#' Has no effect when censored_symbol = NULL, since no values are considered censored. -#' @param numberOfCores Number of cores for parallel processing. When > 1, -#' a logfile named `MSstats_dataProcess_log_progress.log` is created to -#' track progress. Only works for Linux & Mac OS. Default is 1. -#' @param aft_iterations Number of iterations for AFT model fitting. Default is 90. -#' -#' @importFrom parallel makeCluster parLapply stopCluster clusterExport -#' -#' @return list of length one with run-level data. -#' -MSstatsSummarizeWithMultipleCores = function(input, method, impute, censored_symbol, - remove50missing, equal_variance, numberOfCores = 1, - aft_iterations = 90) { - if (numberOfCores > 1) { - is_labeled_reference = "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE) - if (is_labeled_reference) { - protein_indices = split(seq_len(nrow(input)), list(input$PROTEIN)) - } else { - protein_indices = split(seq_len(nrow(input)), list(input$PROTEIN, input$LABEL)) - } - num_proteins = length(protein_indices) - function_environment = environment() - cl = parallel::makeCluster(numberOfCores) - getOption("MSstatsLog")("INFO", - "Starting the cluster setup for summarization") - parallel::clusterExport(cl, c("MSstatsSummarizeSingleTMP", - "MSstatsSummarizeSingleLinear", - "input", "impute", "censored_symbol", - "remove50missing", "protein_indices", - "equal_variance", "aft_iterations"), - envir = function_environment) - cat(paste0("Number of proteins to process: ", num_proteins), - sep = "\n", file = "MSstats_dataProcess_log_progress.log") - if (method == "TMP") { - summarized_results = parallel::parLapply(cl, seq_len(num_proteins), function(i) { - if (i %% 100 == 0) { - cat("Finished processing an additional 100 proteins", - sep = "\n", file = "MSstats_dataProcess_log_progress.log", append = TRUE) - } - single_protein = input[protein_indices[[i]],] - MSstatsSummarizeSingleTMP( - single_protein, impute, censored_symbol, remove50missing, - aft_iterations) - }) - } else { - summarized_results = parallel::parLapply(cl, seq_len(num_proteins), function(i) { - if (i %% 100 == 0) { - cat("Finished processing an additional 100 proteins", - sep = "\n", file = "MSstats_dataProcess_log_progress.log", append = TRUE) - } - single_protein = input[protein_indices[[i]],] - MSstatsSummarizeSingleLinear( - single_protein, - impute, - censored_symbol, - remove50missing, - aft_iterations) - }) - } - parallel::stopCluster(cl) - return(summarized_results) - } else { - return(MSstatsSummarizeWithSingleCore(input, method, impute, - censored_symbol, - remove50missing, - equal_variance, - aft_iterations)) - } -} - #' Feature-level data summarization with 1 core #' #' @inheritParams MSstatsSummarizeWithMultipleCores diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index dfd7749d..dcf9ee4c 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -1,8 +1,8 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/dataProcess.R +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R \name{MSstatsSummarizeWithMultipleCores} \alias{MSstatsSummarizeWithMultipleCores} -\title{Feature-level data summarization with multiple cores} +\title{Feature-level data summarization via socket-dispatched protein records (V6)} \usage{ MSstatsSummarizeWithMultipleCores( input, @@ -11,8 +11,12 @@ MSstatsSummarizeWithMultipleCores( censored_symbol, remove50missing, equal_variance, - numberOfCores = 1, - aft_iterations = 90 + numberOfCores = 1L, + aft_iterations = 90L, + verbose = FALSE, + BPPARAM = NULL, + track_memory = FALSE, + max_proteins_per_worker = 0L ) } \arguments{ @@ -48,11 +52,46 @@ heterogeneous variation from different features.} a logfile named `MSstats_dataProcess_log_progress.log` is created to track progress. Only works for Linux & Mac OS. Default is 1.} -\item{aft_iterations}{Number of iterations for AFT model fitting. Default is 90.} +\item{max_proteins_per_worker}{integer; caps how many protein records +(packed vector + its own metadata) are bundled into a single +\code{bplapply} task sent to one worker. Translated into +\code{SnowfastParam(tasks = ...)}: with \code{N} proteins this becomes +\code{tasks = ceiling(N / max_proteins_per_worker)}. Only applied when +\code{BPPARAM} is \code{NULL}; ignored if the caller supplies +\code{BPPARAM} directly. Default \code{0} reproduces \code{tasks = 0}: +\code{X} is divided as evenly as possible across \code{numberOfCores} +workers.} } \value{ -list of length one with run-level data. +A named list with one element per protein slot, identical in + structure to \code{MSstatsSummarizeWithMultipleCores}. } \description{ -Feature-level data summarization with multiple cores +Fixes a hidden RAM cost in \code{MSstatsSummarizeWithMultipleCoresV5}: V5's +worker closure captured \code{meta_list} — the metadata for \emph{every} +protein — by reference, so each worker received and retained the full +metadata set for the whole run regardless of how many proteins were +actually assigned to it. Bounding the per-task \emph{packed-vector} payload +(e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce that, +since the metadata was baked into the closure once, not re-sliced per task. +} +\details{ +V6 pairs each protein's packed double vector with only its own metadata +into a single list element (\code{list(packed = ..., meta = ...)}) and +dispatches that combined list through \code{bplapply}. \code{BiocParallel} +then only serializes and sends the records actually assigned to a given +task, so a worker never holds metadata for proteins outside its own +task(s) — unlike V5, RAM scales with the batch actually being processed. + +Progress is reported by turning on \code{SnowfastParam}'s built-in +\code{progressbar}, not by having workers talk back to the parent. +\code{BiocParallel} already ticks that progress bar from inside the +manager process itself, once per task result it collects — the same +receive that would happen regardless, over the same single +\code{bplapply} call. Workers never see the flag and send nothing extra +because of it, so this adds no serialization and no IPC beyond what an +unmonitored run already does. Reporting granularity therefore tracks +however many tasks the run is already split into (see +\code{max_proteins_per_worker} below): with the default \code{tasks = 0} +that's one step per worker. } diff --git a/man/dot-MSstatsSummarizeSingleTMPV2.Rd b/man/dot-MSstatsSummarizeSingleTMPV2.Rd new file mode 100644 index 00000000..fb1f9060 --- /dev/null +++ b/man/dot-MSstatsSummarizeSingleTMPV2.Rd @@ -0,0 +1,48 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.MSstatsSummarizeSingleTMPV2} +\alias{.MSstatsSummarizeSingleTMPV2} +\title{Summarize a single protein with TMP directly from a V3 packed double vector} +\usage{ +.MSstatsSummarizeSingleTMPV2( + packed, + meta, + impute, + censored_symbol, + remove50missing, + aft_iterations = 90L +) +} +\arguments{ +\item{packed}{double vector produced by \code{.buildProteinSlotV3}} + +\item{meta}{metadata list from \code{.buildProteinSlotV3}} + +\item{impute}{logical; impute censored values with AFT survival model} + +\item{censored_symbol}{\code{"0"}, \code{"NA"}, or \code{NULL}} + +\item{remove50missing}{logical; skip proteins where all runs are >50\% +missing} + +\item{aft_iterations}{integer; max AFT iterations} +} +\value{ +\code{list(result_dt, survival_dt)} matching the format of + \code{MSstatsSummarizeSingleTMP} +} +\description{ +Bypasses the long-format \code{data.table} reconstruction and the +\code{dcast} inside \code{.fitTukey}. The FL×R packed matrices are +operated on directly: +\enumerate{ + \item Row/column masking replaces the \code{n_obs}/\code{n_obs_run} filter. + \item AFT survival fitting uses a lazily-built minimal \code{data.table} + (constructed only when \code{impute=TRUE} and censored values are + present). + \item TMP is applied via a vectorized scatter into a + \code{(LABEL×RUN) × FEATURE} wide matrix followed by + \code{median_polish_summary} — no \code{dcast} round-trip. +} +} +\keyword{internal} diff --git a/man/dot-buildProteinSlotV3.Rd b/man/dot-buildProteinSlotV3.Rd new file mode 100644 index 00000000..2b576d70 --- /dev/null +++ b/man/dot-buildProteinSlotV3.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.buildProteinSlotV3} +\alias{.buildProteinSlotV3} +\title{Build the packed double vector and lightweight metadata for one protein slot} +\usage{ +.buildProteinSlotV3(dt, slot_k, all_runs) +} +\arguments{ +\item{dt}{data.table rows belonging to one protein (or protein × label) slot} + +\item{slot_k}{integer position of this slot in the global protein list} + +\item{all_runs}{character vector of all run names in global order} +} +\value{ +list with elements \code{packed} (double vector) and \code{meta} (list) +} +\description{ +Build the packed double vector and lightweight metadata for one protein slot +} +\keyword{internal} diff --git a/man/dot-buildSummarizeWorkerV6.Rd b/man/dot-buildSummarizeWorkerV6.Rd new file mode 100644 index 00000000..acd91a5e --- /dev/null +++ b/man/dot-buildSummarizeWorkerV6.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.buildSummarizeWorkerV6} +\alias{.buildSummarizeWorkerV6} +\title{Build the per-record worker closure for +\code{MSstatsSummarizeWithMultipleCoresV6}} +\usage{ +.buildSummarizeWorkerV6( + use_TMP, + impute, + censored_symbol, + remove50missing, + aft_iterations, + equal_variance +) +} +\description{ +Defined at package top level — not nested via \code{local()} inside +\code{MSstatsSummarizeWithMultipleCoresV6} — so the returned closure's +enclosing environment chain is this factory's own (small) call frame plus +the package namespace. \code{BiocParallel}/\code{matter} serialize a +closure's entire enclosing environment chain to ship it to each socket +worker, not just the variables the closure body actually references. A +closure built via \code{local()} inside +\code{MSstatsSummarizeWithMultipleCoresV6} would have that function's own +evaluation frame in its chain — which holds \code{input}, +\code{protein_records}, and other run-scale objects — so every worker +would receive and retain a serialized copy of them even though +\code{.worker} never touches them. Building the closure here instead keeps +its captured state to only the scalar run parameters. +} +\keyword{internal} diff --git a/man/dot-reconstructProteinDTV3.Rd b/man/dot-reconstructProteinDTV3.Rd new file mode 100644 index 00000000..099a6f28 --- /dev/null +++ b/man/dot-reconstructProteinDTV3.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.reconstructProteinDTV3} +\alias{.reconstructProteinDTV3} +\title{Reconstruct a per-protein data.table from a V3 packed double vector} +\usage{ +.reconstructProteinDTV3(packed, meta) +} +\arguments{ +\item{packed}{double vector produced by \code{.buildProteinSlotV3}} + +\item{meta}{metadata list from \code{.buildProteinSlotV3}} +} +\value{ +data.table compatible with \code{MSstatsSummarizeSingleTMP} / + \code{MSstatsSummarizeSingleLinear} +} +\description{ +Reconstruct a per-protein data.table from a V3 packed double vector +} +\keyword{internal} From f885e0b3d6e367b3ec72dcaa1f7bec34d58e8536 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 10:12:10 -0400 Subject: [PATCH 02/18] refactor: Rename and reorganize code to be coherent in English --- R/MSstatsSummarizeWithMultipleCores.R | 760 +++++++++++++------------- 1 file changed, 384 insertions(+), 376 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index a36b0cb9..e4ffd6ce 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -2,7 +2,7 @@ # RSS of the current process in MB. # On Linux reads /proc/self/status (VmRSS); elsewhere falls back to gc() counts. -.memMB <- function() { +.current_rss_mb <- function() { if (file.exists("/proc/self/status")) { ln <- readLines("/proc/self/status", warn = FALSE) m <- grep("^VmRSS:", ln, value = TRUE) @@ -15,7 +15,7 @@ # Cross-platform peak-RSS reader. Reflects the true lifetime peak of the # calling process, regardless of when you call it — no polling required. -.peakRSS_MB <- function() { +.peak_rss_mb <- function() { if (file.exists("/proc/self/status")) { # Linux: VmHWM = kernel-maintained peak resident set size ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) @@ -42,50 +42,53 @@ } # Print a formatted memory report to stderr via message(). -# checkpoints: named numeric vector of RSS snapshots (MB). -# worker_mems: numeric vector of per-worker peak RSS values (may be NA). -# elapsed: total wall-clock seconds (NULL to omit). -.printMemReport <- function(fn_name, checkpoints, worker_mems = NULL, - elapsed = NULL) { - w <- 65L - hr <- strrep("─", w) - fmt_mb <- function(x) if (is.na(x)) " n/a" else sprintf("%7.1f", x) - fmt_delta <- function(now, prev) { - if (is.na(now) || is.na(prev)) return("") - sprintf(" (%+.1f MB)", now - prev) +# checkpoints: named numeric vector of RSS snapshots (MB). +# worker_peak_mb: numeric vector of per-worker peak RSS values (may be NA). +# elapsed: total wall-clock seconds (NULL to omit). +.print_memory_report <- function(function_name, checkpoints, worker_peak_mb = NULL, + elapsed = NULL) { + rule_width <- 65L + rule <- strrep("─", rule_width) + format_mb <- function(x) if (is.na(x)) " n/a" else sprintf("%7.1f", x) + format_delta <- function(now, previous_value) { + if (is.na(now) || is.na(previous_value)) return("") + sprintf(" (%+.1f MB)", now - previous_value) } lines <- c( - hr, - sprintf(" MSstats Memory Report — %s", fn_name), - hr, + rule, + sprintf(" MSstats Memory Report — %s", function_name), + rule, sprintf(" %-36s %7s %s", "Checkpoint", "RSS MB", "Delta") ) - prev <- NA_real_ - for (nm in names(checkpoints)) { - val <- checkpoints[[nm]] + previous_value <- NA_real_ + for (checkpoint_name in names(checkpoints)) { + checkpoint_value <- checkpoints[[checkpoint_name]] lines <- c(lines, - sprintf(" %-36s %s%s", nm, fmt_mb(val), fmt_delta(val, prev))) - prev <- val + sprintf(" %-36s %s%s", checkpoint_name, + format_mb(checkpoint_value), + format_delta(checkpoint_value, previous_value))) + previous_value <- checkpoint_value } - if (!is.null(worker_mems)) { - wm <- worker_mems[!is.na(worker_mems)] - if (length(wm)) { + if (!is.null(worker_peak_mb)) { + observed_peaks <- worker_peak_mb[!is.na(worker_peak_mb)] + if (length(observed_peaks)) { lines <- c(lines, "", sprintf(" Worker RSS min / mean / max : %.1f / %.1f / %.1f MB", - min(wm), mean(wm), max(wm))) + min(observed_peaks), mean(observed_peaks), max(observed_peaks))) } } if (!is.null(elapsed)) lines <- c(lines, sprintf(" Total elapsed: %.1f s", as.numeric(elapsed))) - lines <- c(lines, hr) + lines <- c(lines, rule) message(paste(lines, collapse = "\n")) } -## ── V3 internal helpers ─────────────────────────────────────────────────────── +## ── Protein-slot pack/unpack helpers ─────────────────────────────────────────── ## -## Packed double-vector layout for one protein slot (all column-major matrices): +## Packed double-vector layout for one protein slot (all column-major matrices). +## FL = n_feature_labels (unique FEATURE × LABEL pairs), R = n_runs: ## -## pos 1 : slot_k (protein index, cast to double) +## pos 1 : slot_index (protein index, cast to double) ## pos 2 .. FL*R+1 : newABUNDANCE (FL × R) ## pos FL*R+2 .. 2*FL*R+1 : ABUNDANCE (FL × R; NA for unlabeled/TMP) ## pos 2*FL*R+2 .. 3*FL*R+1 : censored (FL × R; 0.0/1.0) @@ -100,200 +103,201 @@ #' Build the packed double vector and lightweight metadata for one protein slot #' -#' @param dt data.table rows belonging to one protein (or protein × label) slot -#' @param slot_k integer position of this slot in the global protein list +#' @param protein_dt data.table rows belonging to one protein (or protein × label) slot +#' @param slot_index integer position of this slot in the global protein list #' @param all_runs character vector of all run names in global order #' @return list with elements \code{packed} (double vector) and \code{meta} (list) #' @keywords internal -.buildProteinSlotV3 <- function(dt, slot_k, all_runs) { - - R <- length(all_runs) - +.pack_protein_slot <- function(protein_dt, slot_index, all_runs) { + + n_runs <- length(all_runs) + # ── Unique (FEATURE, LABEL) → "effective features", stable ordering ─────── - has_peptide <- "PEPTIDE" %in% colnames(dt) - flp <- unique(dt[, .( + has_peptide <- "PEPTIDE" %in% colnames(protein_dt) + feature_label_dt <- unique(protein_dt[, .( FEATURE = as.character(FEATURE), LABEL = as.character(LABEL), PEPTIDE = if (has_peptide) as.character(PEPTIDE) else as.character(FEATURE) )]) - data.table::setorder(flp, FEATURE, LABEL) - FL <- nrow(flp) - + data.table::setorder(feature_label_dt, FEATURE, LABEL) + n_feature_labels <- nrow(feature_label_dt) + # ── Index maps: feature-label → row index, run name → column index ───────── - fi <- match( - paste(as.character(dt$FEATURE), as.character(dt$LABEL), sep = "\t"), - paste(flp$FEATURE, flp$LABEL, sep = "\t")) - ri_map <- seq_len(R); names(ri_map) <- all_runs - ri <- ri_map[as.character(dt$RUN)] - valid <- !is.na(fi) & !is.na(ri) - - # ── Helper: allocate FL × R matrix and fill from dt column ─────────────── - fill_mat <- function(col_vec, fill = NA_real_) { - m <- matrix(fill, nrow = FL, ncol = R) - m[cbind(fi[valid], ri[valid])] <- as.double(col_vec[valid]) + feature_row_idx <- match( + paste(as.character(protein_dt$FEATURE), as.character(protein_dt$LABEL), sep = "\t"), + paste(feature_label_dt$FEATURE, feature_label_dt$LABEL, sep = "\t")) + run_col_idx_lookup <- seq_len(n_runs); names(run_col_idx_lookup) <- all_runs + run_col_idx <- run_col_idx_lookup[as.character(protein_dt$RUN)] + row_is_valid <- !is.na(feature_row_idx) & !is.na(run_col_idx) + + # ── Helper: allocate n_feature_labels × n_runs matrix, fill from a protein_dt column ── + scatter_into_matrix <- function(col_vec, fill = NA_real_) { + m <- matrix(fill, nrow = n_feature_labels, ncol = n_runs) + m[cbind(feature_row_idx[row_is_valid], run_col_idx[row_is_valid])] <- + as.double(col_vec[row_is_valid]) m } - + # ── Build the five numeric matrices ─────────────────────────────────────── - mat_newABU <- fill_mat(dt$newABUNDANCE) - - has_ABUNDANCE <- "ABUNDANCE" %in% colnames(dt) - mat_ABU <- if (has_ABUNDANCE) fill_mat(dt$ABUNDANCE) else - matrix(NA_real_, FL, R) - - has_censored <- "censored" %in% colnames(dt) - mat_cens <- if (has_censored) fill_mat(as.double(dt$censored)) else - matrix(0.0, FL, R) # treat as non-censored when column absent - - has_cen <- "cen" %in% colnames(dt) - mat_cen <- if (has_cen) fill_mat(dt$cen) else - matrix(NA_real_, FL, R) - - has_anom <- "ANOMALYSCORES" %in% colnames(dt) && - !all(is.na(dt$ANOMALYSCORES)) - mat_anom <- if (has_anom) fill_mat(dt$ANOMALYSCORES) else - matrix(NA_real_, FL, R) - + new_abundance_mat <- scatter_into_matrix(protein_dt$newABUNDANCE) + + has_ABUNDANCE <- "ABUNDANCE" %in% colnames(protein_dt) + abundance_mat <- if (has_ABUNDANCE) scatter_into_matrix(protein_dt$ABUNDANCE) else + matrix(NA_real_, n_feature_labels, n_runs) + + has_censored <- "censored" %in% colnames(protein_dt) + censored_mat <- if (has_censored) scatter_into_matrix(as.double(protein_dt$censored)) else + matrix(0.0, n_feature_labels, n_runs) # treat as non-censored when column absent + + has_cen <- "cen" %in% colnames(protein_dt) + event_mat <- if (has_cen) scatter_into_matrix(protein_dt$cen) else + matrix(NA_real_, n_feature_labels, n_runs) + + has_anom <- "ANOMALYSCORES" %in% colnames(protein_dt) && + !all(is.na(protein_dt$ANOMALYSCORES)) + anomaly_scores_mat <- if (has_anom) scatter_into_matrix(protein_dt$ANOMALYSCORES) else + matrix(NA_real_, n_feature_labels, n_runs) + # ── Per-feature scalar: n_obs (constant within FEATURE × LABEL) ─────────── - n_obs_by_fl <- dt[, .(n_obs = as.double(n_obs[1L])), + n_obs_by_feature_label <- protein_dt[, .(n_obs = as.double(n_obs[1L])), by = .(FEATURE = as.character(FEATURE), LABEL = as.character(LABEL))] - flp_nobs <- n_obs_by_fl[flp, on = c("FEATURE", "LABEL")] - data.table::setorder(flp_nobs, FEATURE, LABEL) - n_obs_vec <- flp_nobs$n_obs - + feature_label_with_nobs <- n_obs_by_feature_label[feature_label_dt, on = c("FEATURE", "LABEL")] + data.table::setorder(feature_label_with_nobs, FEATURE, LABEL) + n_obs_vec <- feature_label_with_nobs$n_obs + # ── Per-run scalars: n_obs_run, prop_features (constant within RUN) ─────── - run_vals <- dt[, .(n_obs_run = as.double(n_obs_run[1L]), + run_level_scalars <- protein_dt[, .(n_obs_run = as.double(n_obs_run[1L]), prop_features = as.double(prop_features[1L])), by = .(RUN = as.character(RUN))] - run_full <- run_vals[data.table::data.table(RUN = all_runs), on = "RUN"] - n_obs_run_vec <- run_full$n_obs_run - prop_features_vec <- run_full$prop_features - + run_scalars_all_runs <- run_level_scalars[data.table::data.table(RUN = all_runs), on = "RUN"] + n_obs_run_vec <- run_scalars_all_runs$n_obs_run + prop_features_vec <- run_scalars_all_runs$prop_features + # ── Pack ───────────────────────────────────────────────────────────────── packed <- c( - as.double(slot_k), - as.vector(mat_newABU), - as.vector(mat_ABU), - as.vector(mat_cens), - as.vector(mat_cen), - as.vector(mat_anom), + as.double(slot_index), + as.vector(new_abundance_mat), + as.vector(abundance_mat), + as.vector(censored_mat), + as.vector(event_mat), + as.vector(anomaly_scores_mat), n_obs_vec, n_obs_run_vec, prop_features_vec ) - + # ── Metadata (string labels, flags; kept in main-process RAM) ──────────── meta <- list( - PROTEIN = as.character(dt$PROTEIN[1L]), - feat_label_pep = as.data.frame(flp), # FL × 3: FEATURE, LABEL, PEPTIDE + PROTEIN = as.character(protein_dt$PROTEIN[1L]), + feature_label_dt = as.data.frame(feature_label_dt), # n_feature_labels × 3: FEATURE, LABEL, PEPTIDE runs = all_runs, - FL = FL, - R = R, - is_labeled_ref = "is_labeled_ref" %in% colnames(dt) && - isTRUE(any(dt$is_labeled_ref, na.rm = TRUE)), + n_feature_labels = n_feature_labels, + n_runs = n_runs, + is_labeled_ref = "is_labeled_ref" %in% colnames(protein_dt) && + isTRUE(any(protein_dt$is_labeled_ref, na.rm = TRUE)), has_ABUNDANCE = has_ABUNDANCE, has_censored = has_censored, has_cen = has_cen, has_anom = has_anom, - add_ref_covariate = "ref_covariate" %in% colnames(dt) + add_ref_covariate = "ref_covariate" %in% colnames(protein_dt) ) - + list(packed = packed, meta = meta) } -#' Reconstruct a per-protein data.table from a V3 packed double vector +#' Reconstruct a per-protein data.table from a packed double vector #' -#' @param packed double vector produced by \code{.buildProteinSlotV3} -#' @param meta metadata list from \code{.buildProteinSlotV3} +#' @param packed double vector produced by \code{.pack_protein_slot} +#' @param meta metadata list from \code{.pack_protein_slot} #' @return data.table compatible with \code{MSstatsSummarizeSingleTMP} / #' \code{MSstatsSummarizeSingleLinear} #' @keywords internal -.reconstructProteinDTV3 <- function(packed, meta) { - - FL <- meta$FL - R <- meta$R - n_mat <- FL * R - - # ── Unpack sections (1-indexed; position 1 is slot_k header) ───────────── - i <- 2L - extract_mat <- function() { - m <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R) - i <<- i + n_mat +.unpack_protein_slot <- function(packed, meta) { + + n_feature_labels <- meta$n_feature_labels + n_runs <- meta$n_runs + matrix_len <- n_feature_labels * n_runs + + # ── Unpack sections (1-indexed; position 1 is slot_index header) ───────── + cursor <- 2L + read_next_matrix <- function() { + m <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs) + cursor <<- cursor + matrix_len m } - mat_newABU <- extract_mat() - mat_ABU <- extract_mat() - mat_cens <- extract_mat() - mat_cen <- extract_mat() - mat_anom <- extract_mat() - - n_obs_vec <- packed[i:(i + FL - 1L)]; i <- i + FL - n_obs_run_vec <- packed[i:(i + R - 1L)]; i <- i + R - prop_feat_vec <- packed[i:(i + R - 1L)] - + new_abundance_mat <- read_next_matrix() + abundance_mat <- read_next_matrix() + censored_mat <- read_next_matrix() + event_mat <- read_next_matrix() + anomaly_scores_mat <- read_next_matrix() + + n_obs_vec <- packed[cursor:(cursor + n_feature_labels - 1L)]; cursor <- cursor + n_feature_labels + n_obs_run_vec <- packed[cursor:(cursor + n_runs - 1L)]; cursor <- cursor + n_runs + prop_features_vec <- packed[cursor:(cursor + n_runs - 1L)] + # ── Column-major melt: rep/each match the matrix column-major ordering ──── - flp <- meta$feat_label_pep # data.frame: FEATURE, LABEL, PEPTIDE (FL rows) + feature_label_dt <- meta$feature_label_dt # data.frame: FEATURE, LABEL, PEPTIDE (n_feature_labels rows) runs <- meta$runs - n_rows <- FL * R - - dt <- data.table::data.table( + n_rows <- n_feature_labels * n_runs + + protein_dt <- data.table::data.table( PROTEIN = rep(meta$PROTEIN, n_rows), - FEATURE = rep(flp$FEATURE, times = R), - LABEL = rep(flp$LABEL, times = R), - PEPTIDE = rep(flp$PEPTIDE, times = R), - RUN = rep(runs, each = FL), - newABUNDANCE = as.vector(mat_newABU), - n_obs = as.integer(rep(n_obs_vec, times = R)), - n_obs_run = as.integer(rep(n_obs_run_vec, each = FL)), - prop_features = rep(prop_feat_vec, each = FL) + FEATURE = rep(feature_label_dt$FEATURE, times = n_runs), + LABEL = rep(feature_label_dt$LABEL, times = n_runs), + PEPTIDE = rep(feature_label_dt$PEPTIDE, times = n_runs), + RUN = rep(runs, each = n_feature_labels), + newABUNDANCE = as.vector(new_abundance_mat), + n_obs = as.integer(rep(n_obs_vec, times = n_runs)), + n_obs_run = as.integer(rep(n_obs_run_vec, each = n_feature_labels)), + prop_features = rep(prop_features_vec, each = n_feature_labels) ) - + # Optional: ABUNDANCE (labeled linear model uses raw log-intensities) if (meta$has_ABUNDANCE) { - dt[, ABUNDANCE := as.vector(mat_ABU)] + protein_dt[, ABUNDANCE := as.vector(abundance_mat)] } - + # censored: stored as 0.0/1.0 doubles; NA → treat as non-censored - dt[, censored := { - v <- as.vector(mat_cens) + protein_dt[, censored := { + v <- as.vector(censored_mat) if (meta$has_censored) as.logical(v > 0.5) else rep(FALSE, n_rows) }] - + # cen: survival event indicator (1 = observed, 0 = left-censored) if (meta$has_cen) { - dt[, cen := as.vector(mat_cen)] + protein_dt[, cen := as.vector(event_mat)] } - + # ANOMALYSCORES: NA when unused (linear model skips anomaly weighting) - dt[, ANOMALYSCORES := as.vector(mat_anom)] - + protein_dt[, ANOMALYSCORES := as.vector(anomaly_scores_mat)] + # SRM-specific columns derived from LABEL and RUN if (meta$is_labeled_ref) { - dt[, is_labeled_ref := (LABEL == "H")] + protein_dt[, is_labeled_ref := (LABEL == "H")] if (meta$add_ref_covariate) { - dt[, ref_covariate := factor( + protein_dt[, ref_covariate := factor( data.table::fifelse(LABEL == "L", as.character(RUN), "0"))] } } - - dt + + protein_dt } -## ── V4: matrix-native TMP ───────────────────────────────────────────────────── +## ── Matrix-native TMP ───────────────────────────────────────────────────────── ## -## .MSstatsSummarizeSingleTMPV2 takes the V3 packed double vector + meta list +## .summarize_protein_tmp_from_packed takes the packed double vector + meta list ## directly and produces TMP results without ever building the full long-format ## data.table or calling dcast. ## ## Compared with MSstatsSummarizeSingleTMP the hot path is: -## V3: packed → long-DT (melt) → filter → AFT fit → dcast → TMP -## V4: packed → FL×R matrix → mask → (AFT fit if needed) → t(mat) → TMP +## long-format: packed → long-DT (melt) → filter → AFT fit → dcast → TMP +## packed: packed → FL×R matrix → mask → (AFT fit if needed) → t(mat) → TMP ## ## The dcast savings are most significant for high-feature proteins (>100 ## features) and large run counts (>50 runs), where the intermediate FL*R @@ -301,7 +305,7 @@ ## ───────────────────────────────────────────────────────────────────────────── -#' Summarize a single protein with TMP directly from a V3 packed double vector +#' Summarize a single protein with TMP directly from a packed double vector #' #' Bypasses the long-format \code{data.table} reconstruction and the #' \code{dcast} inside \code{.fitTukey}. The FL×R packed matrices are @@ -316,8 +320,8 @@ #' \code{median_polish_summary} — no \code{dcast} round-trip. #' } #' -#' @param packed double vector produced by \code{.buildProteinSlotV3} -#' @param meta metadata list from \code{.buildProteinSlotV3} +#' @param packed double vector produced by \code{.pack_protein_slot} +#' @param meta metadata list from \code{.pack_protein_slot} #' @param impute logical; impute censored values with AFT survival model #' @param censored_symbol \code{"0"}, \code{"NA"}, or \code{NULL} #' @param remove50missing logical; skip proteins where all runs are >50\% @@ -326,61 +330,61 @@ #' @return \code{list(result_dt, survival_dt)} matching the format of #' \code{MSstatsSummarizeSingleTMP} #' @keywords internal -.MSstatsSummarizeSingleTMPV2 <- function(packed, meta, - impute, censored_symbol, remove50missing, aft_iterations = 90L) +.summarize_protein_tmp_from_packed <- function(packed, meta, + impute, censored_symbol, remove50missing, aft_iterations = 90L) { - FL <- meta$FL - R <- meta$R - n_mat <- FL * R - PROTEIN <- meta$PROTEIN - flp <- meta$feat_label_pep # data.frame: FEATURE, LABEL, PEPTIDE - runs <- meta$runs - + n_feature_labels <- meta$n_feature_labels + n_runs <- meta$n_runs + matrix_len <- n_feature_labels * n_runs + PROTEIN <- meta$PROTEIN + feature_label_dt <- meta$feature_label_dt # data.frame: FEATURE, LABEL, PEPTIDE + runs <- meta$runs + # ── 1. Unpack matrices ───────────────────────────────────────────────────── # Layout (see top-of-file comment block): - # pos 1 : slot_k (skip) + # pos 1 : slot_index (skip) # pos 2..FL*R+1 : newABUNDANCE (FL × R, column-major) # pos FL*R+2.. : ABUNDANCE (skip — not needed for TMP) # pos 2*FL*R+2.. : censored (FL × R) # pos 3*FL*R+2.. : cen (FL × R) # pos 4*FL*R+2.. : ANOMALYSCORES(skip — not needed for TMP) # then: n_obs(FL), n_obs_run(R), prop_features(R) - i <- 2L - mat_newABU <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R); i <- i + n_mat - i <- i + n_mat # skip ABUNDANCE - mat_cens <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R); i <- i + n_mat - mat_cen <- matrix(packed[i:(i + n_mat - 1L)], nrow = FL, ncol = R); i <- i + n_mat - i <- i + n_mat # skip ANOMALYSCORES - n_obs_vec <- packed[i:(i + FL - 1L)]; i <- i + FL - n_obs_run_vec <- packed[i:(i + R - 1L)]; i <- i + R - prop_feat_vec <- packed[i:(i + R - 1L)] - + cursor <- 2L + new_abundance_mat <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs); cursor <- cursor + matrix_len + cursor <- cursor + matrix_len # skip ABUNDANCE + censored_mat <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs); cursor <- cursor + matrix_len + event_mat <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs); cursor <- cursor + matrix_len + cursor <- cursor + matrix_len # skip ANOMALYSCORES + n_obs_vec <- packed[cursor:(cursor + n_feature_labels - 1L)]; cursor <- cursor + n_feature_labels + n_obs_run_vec <- packed[cursor:(cursor + n_runs - 1L)]; cursor <- cursor + n_runs + prop_features_vec <- packed[cursor:(cursor + n_runs - 1L)] + # ── 2. Row/column masking (mirrors the n_obs / n_obs_run filter) ─────────── - keep_feat <- !is.na(n_obs_vec) & n_obs_vec > 1 - keep_runs <- !is.na(n_obs_run_vec) & n_obs_run_vec > 0 - - mat_newABU <- mat_newABU[keep_feat, keep_runs, drop = FALSE] - mat_cens <- mat_cens [keep_feat, keep_runs, drop = FALSE] - mat_cen <- mat_cen [keep_feat, keep_runs, drop = FALSE] - - flp_filt <- flp[keep_feat, , drop = FALSE] - runs_filt <- runs[keep_runs] - prop_feat_filt <- prop_feat_vec[keep_runs] - FL_filt <- nrow(mat_newABU) - R_filt <- ncol(mat_newABU) - + feature_is_kept <- !is.na(n_obs_vec) & n_obs_vec > 1 + run_is_kept <- !is.na(n_obs_run_vec) & n_obs_run_vec > 0 + + new_abundance_mat <- new_abundance_mat[feature_is_kept, run_is_kept, drop = FALSE] + censored_mat <- censored_mat [feature_is_kept, run_is_kept, drop = FALSE] + event_mat <- event_mat [feature_is_kept, run_is_kept, drop = FALSE] + + feature_label_dt_kept <- feature_label_dt[feature_is_kept, , drop = FALSE] + runs_kept <- runs[run_is_kept] + prop_features_kept <- prop_features_vec[run_is_kept] + n_feature_labels_kept <- nrow(new_abundance_mat) + n_runs_kept <- ncol(new_abundance_mat) + # ── 3. Empty-matrix early exit ───────────────────────────────────────────── - if (FL_filt == 0L || R_filt == 0L) { + if (n_feature_labels_kept == 0L || n_runs_kept == 0L) { msg <- paste("Can't summarize for protein", PROTEIN, "because all measurements are missing or censored.") try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) return(list(NULL, NULL)) } - + # ── 4. remove50missing check (independent of imputation) ────────────────── if (remove50missing && - all(prop_feat_filt <= 0.5 | is.na(prop_feat_filt))) { + all(prop_features_kept <= 0.5 | is.na(prop_features_kept))) { msg <- paste("Can't summarize for protein", PROTEIN, "because all runs have more than 50% missing values and", "are removed with the option, remove50missing=TRUE.") @@ -388,36 +392,36 @@ try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) return(list(NULL, NULL)) } - + # ── 5. AFT survival imputation (lazy: only if any censored present) ──────── - n_rows_filt <- FL_filt * R_filt - orig_abu_vec <- as.vector(mat_newABU) - - surv_cols <- c("newABUNDANCE", "cen", "RUN", "FEATURE", "ref_covariate") - any_censored <- meta$has_censored && any(mat_cens > 0.5, na.rm = TRUE) - + n_kept_rows <- n_feature_labels_kept * n_runs_kept + original_abundance_vec <- as.vector(new_abundance_mat) + + survival_columns <- c("newABUNDANCE", "cen", "RUN", "FEATURE", "ref_covariate") + any_censored <- meta$has_censored && any(censored_mat > 0.5, na.rm = TRUE) + if (impute && any_censored && meta$has_cen) { # Build minimal long-format DT — only columns needed by .fitSurvival - surv_dt <- data.table::data.table( - newABUNDANCE = orig_abu_vec, - cen = as.vector(mat_cen), - censored = as.logical(as.vector(mat_cens) > 0.5), - FEATURE = factor(rep(flp_filt$FEATURE, times = R_filt)), - RUN = factor(rep(runs_filt, each = FL_filt)), - LABEL = rep(as.character(flp_filt$LABEL), times = R_filt) + survival_dt <- data.table::data.table( + newABUNDANCE = original_abundance_vec, + cen = as.vector(event_mat), + censored = as.logical(as.vector(censored_mat) > 0.5), + FEATURE = factor(rep(feature_label_dt_kept$FEATURE, times = n_runs_kept)), + RUN = factor(rep(runs_kept, each = n_feature_labels_kept)), + LABEL = rep(as.character(feature_label_dt_kept$LABEL), times = n_runs_kept) ) if (meta$add_ref_covariate) { - surv_dt[, ref_covariate := factor( + survival_dt[, ref_covariate := factor( data.table::fifelse(LABEL == "L", as.character(RUN), "0"))] } - - fit_cols <- intersect(surv_cols, colnames(surv_dt)) + + fit_columns <- intersect(survival_columns, colnames(survival_dt)) fit_data <- if (meta$is_labeled_ref) { - surv_dt[LABEL != "H", fit_cols, with = FALSE] + survival_dt[LABEL != "H", fit_columns, with = FALSE] } else { - surv_dt[, fit_cols, with = FALSE] + survival_dt[, fit_columns, with = FALSE] } - + converged <- TRUE survival_fit <- withCallingHandlers({ MSstats:::.fitSurvival(fit_data, aft_iterations) @@ -427,133 +431,133 @@ converged <<- FALSE } }) - - predicted_all <- if (converged) { - predict(survival_fit, newdata = surv_dt) + + predicted_values <- if (converged) { + predict(survival_fit, newdata = survival_dt) } else { - rep(NA_real_, n_rows_filt) + rep(NA_real_, n_kept_rows) } - - use_imputed <- if (meta$is_labeled_ref) { - surv_dt$censored & surv_dt$LABEL != "H" + + should_impute <- if (meta$is_labeled_ref) { + survival_dt$censored & survival_dt$LABEL != "H" } else { - surv_dt$censored + survival_dt$censored } - - imputed_abu <- ifelse(use_imputed, predicted_all, orig_abu_vec) - surv_dt[, predicted := ifelse(use_imputed, predicted_all, NA_real_)] - surv_dt[, newABUNDANCE := imputed_abu] - mat_newABU <- matrix(imputed_abu, nrow = FL_filt, ncol = R_filt) - survival <- surv_dt[, intersect(c(surv_cols, "LABEL", "predicted"), - colnames(surv_dt)), with = FALSE] - + + imputed_abundance_vec <- ifelse(should_impute, predicted_values, original_abundance_vec) + survival_dt[, predicted := ifelse(should_impute, predicted_values, NA_real_)] + survival_dt[, newABUNDANCE := imputed_abundance_vec] + new_abundance_mat <- matrix(imputed_abundance_vec, nrow = n_feature_labels_kept, ncol = n_runs_kept) + survival_output <- survival_dt[, intersect(c(survival_columns, "LABEL", "predicted"), + colnames(survival_dt)), with = FALSE] + } else { # No imputation: build survival DT from original values - surv_dt <- data.table::data.table( - newABUNDANCE = orig_abu_vec, - cen = if (meta$has_cen) as.vector(mat_cen) - else rep(NA_real_, n_rows_filt), - FEATURE = rep(as.character(flp_filt$FEATURE), times = R_filt), - RUN = rep(runs_filt, each = FL_filt), - LABEL = rep(as.character(flp_filt$LABEL), times = R_filt), + survival_dt <- data.table::data.table( + newABUNDANCE = original_abundance_vec, + cen = if (meta$has_cen) as.vector(event_mat) + else rep(NA_real_, n_kept_rows), + FEATURE = rep(as.character(feature_label_dt_kept$FEATURE), times = n_runs_kept), + RUN = rep(runs_kept, each = n_feature_labels_kept), + LABEL = rep(as.character(feature_label_dt_kept$LABEL), times = n_runs_kept), predicted = NA # logical NA, matching MSstatsSummarizeSingleTMP's # bare `survival[, predicted := NA]` in its # equivalent non-imputed branch (dataProcess.R:586) ) - survival <- surv_dt[, intersect(c(surv_cols, "LABEL", "predicted"), - colnames(surv_dt)), with = FALSE] + survival_output <- survival_dt[, intersect(c(survival_columns, "LABEL", "predicted"), + colnames(survival_dt)), with = FALSE] } - + # ── 6. Post-imputation .isSummarizable check ─────────────────────────────── - if (all(is.na(mat_newABU) | mat_newABU == 0)) { + if (all(is.na(new_abundance_mat) | new_abundance_mat == 0)) { msg <- paste("Can't summarize for protein", PROTEIN, "because all measurements are missing or censored.") try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) return(list(NULL, NULL)) } - + # ── 7. TMP via vectorized scatter into (LABEL×RUN) × FEATURE wide matrix ─── # # Mirrors .fitTukey: dcast(LABEL + RUN ~ FEATURE, value.var="newABUNDANCE") # Row order in wide_mat: (lbl1,run1),(lbl1,run2),...,(lbl2,run1),... # where labels are sorted — matches dcast alphabetical ordering for # character LABEL values. - if (FL_filt == 1L) { + if (n_feature_labels_kept == 1L) { # ── Single (FEATURE, LABEL) row: skip TMP, use values directly ──────── if (meta$is_labeled_ref) { - h_sel <- as.character(flp_filt$LABEL) == "H" - l_sel <- as.character(flp_filt$LABEL) == "L" - if (any(h_sel) && any(l_sel)) { - h_vals <- as.vector(mat_newABU[h_sel, , drop = FALSE]) - l_vals <- as.vector(mat_newABU[l_sel, , drop = FALSE]) - h_median <- stats::median(h_vals, na.rm = TRUE) + is_heavy_row <- as.character(feature_label_dt_kept$LABEL) == "H" + is_light_row <- as.character(feature_label_dt_kept$LABEL) == "L" + if (any(is_heavy_row) && any(is_light_row)) { + heavy_values <- as.vector(new_abundance_mat[is_heavy_row, , drop = FALSE]) + light_values <- as.vector(new_abundance_mat[is_light_row, , drop = FALSE]) + heavy_median <- stats::median(heavy_values, na.rm = TRUE) result <- data.table::data.table( Protein = PROTEIN, LABEL = "L", - RUN = runs_filt, - LogIntensities = l_vals - h_vals + h_median + RUN = runs_kept, + LogIntensities = light_values - heavy_values + heavy_median ) } else { result <- data.table::data.table( Protein = PROTEIN, - LABEL = as.character(flp_filt$LABEL[1L]), - RUN = runs_filt, - LogIntensities = as.vector(mat_newABU) + LABEL = as.character(feature_label_dt_kept$LABEL[1L]), + RUN = runs_kept, + LogIntensities = as.vector(new_abundance_mat) ) } } else { result <- data.table::data.table( Protein = PROTEIN, - LABEL = rep(as.character(flp_filt$LABEL), times = R_filt), - RUN = rep(runs_filt, each = FL_filt), - LogIntensities = as.vector(mat_newABU) + LABEL = rep(as.character(feature_label_dt_kept$LABEL), times = n_runs_kept), + RUN = rep(runs_kept, each = n_feature_labels_kept), + LogIntensities = as.vector(new_abundance_mat) ) } } else { # ── Multi-feature: scatter into wide matrix, apply TMP ───────────────── - unique_labels <- sort(unique(as.character(flp_filt$LABEL))) - unique_feats <- sort(unique(as.character(flp_filt$FEATURE))) - n_lbl <- length(unique_labels) - n_feat <- length(unique_feats) - - lbl_idx_map <- match(as.character(flp_filt$LABEL), unique_labels) - feat_idx_map <- match(as.character(flp_filt$FEATURE), unique_feats) - - # Vectorized scatter: mat[fl, run] → wide[(lbl_idx-1)*R + run, feat_idx] - fl_rep <- rep(seq_len(FL_filt), times = R_filt) - run_rep <- rep(seq_len(R_filt), each = FL_filt) - wide_row <- (lbl_idx_map[fl_rep] - 1L) * R_filt + run_rep - wide_col <- feat_idx_map[fl_rep] - - wide_mat <- matrix(NA_real_, nrow = n_lbl * R_filt, ncol = n_feat) - wide_mat[cbind(wide_row, wide_col)] <- as.vector(mat_newABU) - - tmp_fitted <- MSstats:::median_polish_summary(wide_mat) - - # Row → (LABEL, RUN) mapping: label index cycles every R_filt rows - result_labels <- rep(unique_labels, each = R_filt) - result_runs <- rep(runs_filt, times = n_lbl) - + unique_labels <- sort(unique(as.character(feature_label_dt_kept$LABEL))) + unique_features <- sort(unique(as.character(feature_label_dt_kept$FEATURE))) + n_labels <- length(unique_labels) + n_features <- length(unique_features) + + label_idx_map <- match(as.character(feature_label_dt_kept$LABEL), unique_labels) + feature_idx_map <- match(as.character(feature_label_dt_kept$FEATURE), unique_features) + + # Vectorized scatter: mat[fl, run] → wide[(lbl_idx-1)*n_runs_kept + run, feat_idx] + feature_idx_of_cell <- rep(seq_len(n_feature_labels_kept), times = n_runs_kept) + run_idx_of_cell <- rep(seq_len(n_runs_kept), each = n_feature_labels_kept) + wide_mat_row <- (label_idx_map[feature_idx_of_cell] - 1L) * n_runs_kept + run_idx_of_cell + wide_mat_col <- feature_idx_map[feature_idx_of_cell] + + wide_mat <- matrix(NA_real_, nrow = n_labels * n_runs_kept, ncol = n_features) + wide_mat[cbind(wide_mat_row, wide_mat_col)] <- as.vector(new_abundance_mat) + + tmp_fitted_values <- MSstats:::median_polish_summary(wide_mat) + + # Row → (LABEL, RUN) mapping: label index cycles every n_runs_kept rows + result_labels <- rep(unique_labels, each = n_runs_kept) + result_runs <- rep(runs_kept, times = n_labels) + if (meta$is_labeled_ref) { - h_li <- match("H", unique_labels) - l_li <- match("L", unique_labels) - if (!is.na(h_li) && !is.na(l_li)) { - h_rows <- (h_li - 1L) * R_filt + seq_len(R_filt) - l_rows <- (l_li - 1L) * R_filt + seq_len(R_filt) - h_med <- stats::median(tmp_fitted[h_rows], na.rm = TRUE) + heavy_label_idx <- match("H", unique_labels) + light_label_idx <- match("L", unique_labels) + if (!is.na(heavy_label_idx) && !is.na(light_label_idx)) { + heavy_rows <- (heavy_label_idx - 1L) * n_runs_kept + seq_len(n_runs_kept) + light_rows <- (light_label_idx - 1L) * n_runs_kept + seq_len(n_runs_kept) + heavy_median <- stats::median(tmp_fitted_values[heavy_rows], na.rm = TRUE) result <- data.table::data.table( Protein = PROTEIN, LABEL = "L", - RUN = runs_filt, - LogIntensities = tmp_fitted[l_rows] - tmp_fitted[h_rows] + h_med + RUN = runs_kept, + LogIntensities = tmp_fitted_values[light_rows] - tmp_fitted_values[heavy_rows] + heavy_median ) } else { result <- data.table::data.table( Protein = PROTEIN, LABEL = result_labels, RUN = result_runs, - LogIntensities = tmp_fitted + LogIntensities = tmp_fitted_values ) } } else { @@ -561,25 +565,25 @@ Protein = PROTEIN, LABEL = result_labels, RUN = result_runs, - LogIntensities = tmp_fitted + LogIntensities = tmp_fitted_values ) } } - - list(result, survival) + + list(result, survival_output) } #' Build the per-record worker closure for -#' \code{MSstatsSummarizeWithMultipleCoresV6} +#' \code{MSstatsSummarizeWithMultipleCores} #' #' Defined at package top level — not nested via \code{local()} inside -#' \code{MSstatsSummarizeWithMultipleCoresV6} — so the returned closure's +#' \code{MSstatsSummarizeWithMultipleCores} — so the returned closure's #' enclosing environment chain is this factory's own (small) call frame plus #' the package namespace. \code{BiocParallel}/\code{matter} serialize a #' closure's entire enclosing environment chain to ship it to each socket #' worker, not just the variables the closure body actually references. A #' closure built via \code{local()} inside -#' \code{MSstatsSummarizeWithMultipleCoresV6} would have that function's own +#' \code{MSstatsSummarizeWithMultipleCores} would have that function's own #' evaluation frame in its chain — which holds \code{input}, #' \code{protein_records}, and other run-scale objects — so every worker #' would receive and retain a serialized copy of them even though @@ -587,69 +591,72 @@ #' its captured state to only the scalar run parameters. #' #' @keywords internal -.buildSummarizeWorkerV6 <- function( +.build_summarize_worker <- function( use_TMP, impute, censored_symbol, remove50missing, aft_iterations, equal_variance ) { - reconstruct_ <- .reconstructProteinDTV3 + unpack_fn <- .unpack_protein_slot use_TMP_ <- use_TMP impute_ <- impute censored_symbol_ <- censored_symbol remove50missing_ <- remove50missing aft_iterations_ <- aft_iterations equal_variance_ <- equal_variance - + function(record) { packed <- record$packed meta <- record$meta result <- if (use_TMP_) { - .MSstatsSummarizeSingleTMPV2( + .summarize_protein_tmp_from_packed( packed, meta, impute_, censored_symbol_, remove50missing_, aft_iterations_) } else { - dt <- reconstruct_(packed, meta) + protein_dt <- unpack_fn(packed, meta) MSstatsSummarizeSingleLinear( - dt, impute_, censored_symbol_, + protein_dt, impute_, censored_symbol_, remove50missing_, aft_iterations_, equal_variances = equal_variance_) } - # Normalize column types/levels at the V6 worker boundary rather than - # inside .MSstatsSummarizeSingleTMPV2/MSstatsSummarizeSingleLinear: - # those are shared with V3-V5 and dataProcess.R respectively, and the - # TMP path in particular carries RUN as plain character (from - # meta$runs) and cen as double (from the packed matrix). Re-leveling - # against meta$runs (rather than a bare factor(RUN)) preserves the - # numeric run order that levels(input$RUN) already established - # upstream — a bare factor() call would instead re-sort alphabetically - # ("1","10","11","2",... for >=10 runs). as.character() first strips - # any existing level order before re-leveling, since .MSstatsSummarize - # SingleTMPV2's survival table (result[[2L]]) already wraps RUN in a - # bare factor() for the imputed branch — that call has the exact same - # alphabetical-resort bug even though it looks pre-typed. droplevels() - # then matches MSstatsSummarizeSingleTMP/SingleCore, whose factor(RUN) - # only ever sees — and so only ever keeps — the runs present for that - # protein, since meta$runs carries every run across the whole input. + # Normalize column types/levels at the worker boundary rather than + # inside .summarize_protein_tmp_from_packed/MSstatsSummarizeSingleLinear: + # those are shared with earlier packed-vector iterations and + # dataProcess.R respectively, and the TMP path in particular carries + # RUN as plain character (from meta$runs) and cen as double (from the + # packed matrix). Re-leveling against meta$runs (rather than a bare + # factor(RUN)) preserves the numeric run order that + # levels(input$RUN) already established upstream — a bare factor() + # call would instead re-sort alphabetically ("1","10","11","2",... + # for >=10 runs). as.character() first strips any existing level + # order before re-leveling, since + # .summarize_protein_tmp_from_packed's survival table (result[[2L]]) + # already wraps RUN in a bare factor() for the imputed branch — that + # call has the exact same alphabetical-resort bug even though it + # looks pre-typed. droplevels() then matches + # MSstatsSummarizeSingleTMP/SingleCore, whose factor(RUN) only ever + # sees — and so only ever keeps — the runs present for that protein, + # since meta$runs carries every run across the whole input. # - # FEATURE gets the same treatment: .MSstatsSummarizeSingleTMPV2's + # FEATURE gets the same treatment: .summarize_protein_tmp_from_packed's # non-imputed survival branch leaves FEATURE as plain character (only # the imputed branch wraps it in a bare factor()), so result[[2L]]$ # FEATURE isn't reliably a factor the way SingleCore's is (single_ # protein[, FEATURE := factor(FEATURE)] before survival is sliced off # it). # - # Deliberately NOT leveling against meta$feat_label_pep$FEATURE's - # existing row order here: that order comes from .buildProteinSlotV3's - # data.table::setorder(flp, FEATURE, LABEL), and data.table sorts - # character columns in the C-locale (byte order: all uppercase before - # any lowercase) for platform-independence — whereas SingleCore's - # bare factor(FEATURE) goes through base R's factor()/sort(), which - # use the session's collation locale (e.g. en_US.UTF-8, where case is - # interleaved: "a" < "A" < "b" < "B"). Any FEATURE strings with mixed - # case (e.g. modification tags) sort differently between the two, so - # re-deriving the level order with a base R sort() reproduces - # SingleCore's order instead of inheriting data.table's. - feature_levels <- sort(unique(meta$feat_label_pep$FEATURE)) + # Deliberately NOT leveling against meta$feature_label_dt$FEATURE's + # existing row order here: that order comes from .pack_protein_slot's + # data.table::setorder(feature_label_dt, FEATURE, LABEL), and + # data.table sorts character columns in the C-locale (byte order: + # all uppercase before any lowercase) for platform-independence — + # whereas SingleCore's bare factor(FEATURE) goes through base R's + # factor()/sort(), which use the session's collation locale (e.g. + # en_US.UTF-8, where case is interleaved: "a" < "A" < "b" < "B"). Any + # FEATURE strings with mixed case (e.g. modification tags) sort + # differently between the two, so re-deriving the level order with a + # base R sort() reproduces SingleCore's order instead of inheriting + # data.table's. + feature_levels <- sort(unique(meta$feature_label_dt$FEATURE)) for (idx in 1:2) { if (!is.null(result[[idx]])) { if ("RUN" %in% colnames(result[[idx]])) @@ -664,7 +671,7 @@ } } -#' Per-worker peak-RSS query task for \code{MSstatsSummarizeWithMultipleCoresV6} +#' Per-worker peak-RSS query task for \code{MSstatsSummarizeWithMultipleCores} #' #' Dispatched once per worker (via \code{seq_len(bpnworkers(BPPARAM))}) after #' the main summarization \code{bplapply} call, while the persistent workers @@ -673,41 +680,43 @@ #' #' @keywords internal #' @noRd -.reportWorkerPeakV6 <- function(i) { - list(worker = i, pid = Sys.getpid(), peak_mb = .peakRSS_MB()) +.report_worker_peak <- function(i) { + list(worker = i, pid = Sys.getpid(), peak_mb = .peak_rss_mb()) } -# Per-worker warm-up task for MSstatsSummarizeWithMultipleCoresV6: loads +# Per-worker warm-up task for MSstatsSummarizeWithMultipleCores: loads # MSstats once per persistent worker process and pins data.table to a single # thread. Defined at top level rather than inline inside -# MSstatsSummarizeWithMultipleCoresV6: even though its body never references +# MSstatsSummarizeWithMultipleCores: even though its body never references # `input`/`protein_records`, a closure defined inline there would still carry # that function's evaluation frame in its enclosing environment chain, and # BiocParallel would serialize that whole frame — including those run-scale # objects — to every worker just to ship this no-op task. -.warmupV6Worker <- function(i) { +.warmup_worker <- function(i) { library(MSstats, quietly = TRUE, warn.conflicts = FALSE) data.table::setDTthreads(1) NULL } -#' Feature-level data summarization via socket-dispatched protein records (V6) +#' Feature-level data summarization via socket-dispatched protein records #' -#' Fixes a hidden RAM cost in \code{MSstatsSummarizeWithMultipleCoresV5}: V5's -#' worker closure captured \code{meta_list} — the metadata for \emph{every} -#' protein — by reference, so each worker received and retained the full -#' metadata set for the whole run regardless of how many proteins were -#' actually assigned to it. Bounding the per-task \emph{packed-vector} payload -#' (e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce that, -#' since the metadata was baked into the closure once, not re-sliced per task. +#' Fixes a hidden RAM cost present in an earlier iteration of this function, +#' where the worker closure captured \code{meta_list} — the metadata for +#' \emph{every} protein — by reference, so each worker received and retained +#' the full metadata set for the whole run regardless of how many proteins +#' were actually assigned to it. Bounding the per-task \emph{packed-vector} +#' payload (e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce +#' that, since the metadata was baked into the closure once, not re-sliced +#' per task. #' -#' V6 pairs each protein's packed double vector with only its own metadata -#' into a single list element (\code{list(packed = ..., meta = ...)}) and -#' dispatches that combined list through \code{bplapply}. \code{BiocParallel} -#' then only serializes and sends the records actually assigned to a given -#' task, so a worker never holds metadata for proteins outside its own -#' task(s) — unlike V5, RAM scales with the batch actually being processed. +#' This version pairs each protein's packed double vector with only its own +#' metadata into a single list element (\code{list(packed = ..., meta = ...)}) +#' and dispatches that combined list through \code{bplapply}. +#' \code{BiocParallel} then only serializes and sends the records actually +#' assigned to a given task, so a worker never holds metadata for proteins +#' outside its own task(s) — unlike that earlier iteration, RAM scales with +#' the batch actually being processed. #' #' Progress is reported by turning on \code{SnowfastParam}'s built-in #' \code{progressbar}, not by having workers talk back to the parent. @@ -756,8 +765,8 @@ #' \code{X} is divided as evenly as possible across \code{numberOfCores} #' workers. #' -#' @return A named list with one element per protein slot, identical in -#' structure to \code{MSstatsSummarizeWithMultipleCores}. +#' @return A named list with one element per protein slot, keyed by protein +#' (or protein \eqn{\times} label) identifier. #' #' @importFrom matter SnowfastParam #' @importFrom BiocParallel bplapply bpstart bpstop bpisup bpnworkers bpprogressbar @@ -786,8 +795,8 @@ MSstatsSummarizeWithMultipleCores <- function( remove50missing, equal_variance, aft_iterations)) } - t_start <- proc.time()[["elapsed"]] - mem_log <- list() + start_time <- proc.time()[["elapsed"]] + memory_checkpoints <- list() # ── 1. Split input by protein slot ──────────────────────────────────────── is_labeled_reference <- "is_labeled_ref" %in% colnames(input) && @@ -810,27 +819,26 @@ MSstatsSummarizeWithMultipleCores <- function( as.character(sort(unique(input$RUN))) getOption("MSstatsLog")("INFO", - paste0("V6: packing ", num_proteins, " proteins × ", + paste0("Packing ", num_proteins, " proteins × ", length(all_runs), " runs into per-protein records")) - + # ── 2. Pack each protein into a (packed vector, own metadata) record ────── # - # Unlike V5's parallel packed_list/meta_list arrays, each protein's - # metadata travels bundled with its own packed vector. bplapply/BPPARAM - # slice this combined list into tasks, so a worker's closure never needs - # — and never receives — metadata for proteins outside its assigned - # task(s). + # Unlike earlier iterations' parallel packed_list/meta_list arrays, each + # protein's metadata travels bundled with its own packed vector. + # bplapply/BPPARAM slice this combined list into tasks, so a worker's + # closure never needs — and never receives — metadata for proteins + # outside its assigned task(s). protein_records <- vector("list", num_proteins) - for (k in seq_len(num_proteins)) { - slot <- .buildProteinSlotV3( - input[protein_indices[[k]], ], k, all_runs) - protein_records[[k]] <- list(packed = slot$packed, meta = slot$meta) + for (slot_index in seq_len(num_proteins)) { + protein_records[[slot_index]] <- .pack_protein_slot( + input[protein_indices[[slot_index]], ], slot_index, all_runs) } payload_mb <- sum(vapply(protein_records, function(r) length(r$packed), integer(1))) * 8 / 1024^2 getOption("MSstatsLog")("INFO", - paste0("V6: dispatching via sockets (", + paste0("Dispatching via sockets (", format(round(payload_mb, 1)), " MB total packed payload; metadata sharded per task)")) @@ -839,14 +847,14 @@ MSstatsSummarizeWithMultipleCores <- function( # No captured meta_list: each task's records already carry their own # metadata, so the closure only needs the scalar run parameters. # - # Built via .buildSummarizeWorkerV6() (defined at package top level) + # Built via .build_summarize_worker() (defined at package top level) # rather than local() here, so the closure's enclosing environment chain # never includes this function's own frame — which holds `input`, # `protein_records`, etc. — and BiocParallel doesn't serialize those # run-scale objects to every worker. use_TMP <- identical(method, "TMP") - - .worker <- .buildSummarizeWorkerV6( + + worker_fn <- .build_summarize_worker( use_TMP, impute, censored_symbol, remove50missing, aft_iterations, equal_variance) @@ -866,7 +874,7 @@ MSstatsSummarizeWithMultipleCores <- function( 0L } getOption("MSstatsLog")("INFO", - paste0("V6: dispatching as ", + paste0("Dispatching as ", if (tasks > 0L) tasks else numberOfCores, " task(s) (max_proteins_per_worker = ", max_proteins_per_worker, ")")) @@ -898,27 +906,27 @@ MSstatsSummarizeWithMultipleCores <- function( BiocParallel::bpprogressbar(BPPARAM) <- FALSE BiocParallel::bplapply( seq_len(BiocParallel::bpnworkers(BPPARAM)), - .warmupV6Worker, BPPARAM = BPPARAM) + .warmup_worker, BPPARAM = BPPARAM) BiocParallel::bpprogressbar(BPPARAM) <- show_progress - results <- BiocParallel::bplapply(protein_records, .worker, BPPARAM = BPPARAM) + results <- BiocParallel::bplapply(protein_records, worker_fn, BPPARAM = BPPARAM) names(results) <- protein_ids - + worker_peaks <- NULL if (track_memory) { BiocParallel::bpprogressbar(BPPARAM) <- FALSE worker_peaks <- BiocParallel::bplapply( seq_len(BiocParallel::bpnworkers(BPPARAM)), - .reportWorkerPeakV6, BPPARAM = BPPARAM) # must run BEFORE bpstop() while workers are alive + .report_worker_peak, BPPARAM = BPPARAM) # must run BEFORE bpstop() while workers are alive BiocParallel::bpprogressbar(BPPARAM) <- show_progress - mem_log[["parent peak (main)"]] <- .peakRSS_MB() + memory_checkpoints[["parent peak (main)"]] <- .peak_rss_mb() worker_peak_mb <- vapply(worker_peaks, function(x) x$peak_mb, numeric(1L)) - .printMemReport( - "MSstatsSummarizeWithMultipleCoresV6", - mem_log, worker_peak_mb, - elapsed = proc.time()[["elapsed"]] - t_start) + .print_memory_report( + "MSstatsSummarizeWithMultipleCores", + memory_checkpoints, worker_peak_mb, + elapsed = proc.time()[["elapsed"]] - start_time) } - - getOption("MSstatsLog")("INFO", "V6: summarization complete.") + + getOption("MSstatsLog")("INFO", "Summarization complete.") results } From 29f34dad9cd249e52c293325e53dffd679925d26 Mon Sep 17 00:00:00 2001 From: tonywu1999 Date: Wed, 5 Aug 2026 14:00:48 -0400 Subject: [PATCH 03/18] add memory tracking for windows --- R/MSstatsSummarizeWithMultipleCores.R | 65 ++++++++++++++++++--------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index e4ffd6ce..38dfae1f 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -15,30 +15,53 @@ # Cross-platform peak-RSS reader. Reflects the true lifetime peak of the # calling process, regardless of when you call it — no polling required. +# All three branches report OS-level peak resident/working-set memory, +# so figures are comparable across platforms. .peak_rss_mb <- function() { - if (file.exists("/proc/self/status")) { - # Linux: VmHWM = kernel-maintained peak resident set size - ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) - if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024) - } - # macOS (also works as a Linux fallback): POSIX getrusage() ru_maxrss - # is likewise a lifetime high-water mark, just different units per OS. - if (!exists(".rusage_maxrss_mb_impl", mode = "function")) { - Rcpp::cppFunction( - depends = "Rcpp", - includes = "#include ", - code = " - double rusage_maxrss_mb_impl() { - struct rusage ru; getrusage(RUSAGE_SELF, &ru); - #ifdef __APPLE__ - return (double) ru.ru_maxrss / (1024.0*1024.0); // bytes -> MB - #else - return (double) ru.ru_maxrss / 1024.0; // KB -> MB - #endif + if (.Platform$OS.type == "windows") { + if (!exists(".peakRSS_windows_impl", mode = "function")) { + Rcpp::cppFunction( + depends = "Rcpp", + includes = c( + "#include ", + "#include "), + code = " + double peakRSS_windows_impl() { + PROCESS_MEMORY_COUNTERS pmc; + if (GetProcessMemoryInfo(GetCurrentProcess(), &pmc, sizeof(pmc))) { + return (double) pmc.PeakWorkingSetSize / (1024.0*1024.0); // bytes -> MB + } + return NA_REAL; }") - assign(".rusage_maxrss_mb_impl", rusage_maxrss_mb_impl, envir = .GlobalEnv) + assign(".peakRSS_windows_impl", peakRSS_windows_impl, envir = .GlobalEnv) } - .rusage_maxrss_mb_impl() + return(.peakRSS_windows_impl()) + } + + if (file.exists("/proc/self/status")) { + # Linux: VmHWM = kernel-maintained peak resident set size + ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) + if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024) + } + + # macOS (also works as a Linux fallback): POSIX getrusage() ru_maxrss + # is likewise a lifetime high-water mark, just different units per OS. + if (!exists(".rusage_maxrss_mb_impl", mode = "function")) { + Rcpp::cppFunction( + depends = "Rcpp", + includes = "#include ", + code = " + double rusage_maxrss_mb_impl() { + struct rusage ru; getrusage(RUSAGE_SELF, &ru); + #ifdef __APPLE__ + return (double) ru.ru_maxrss / (1024.0*1024.0); // bytes -> MB + #else + return (double) ru.ru_maxrss / 1024.0; // KB -> MB + #endif + }") + assign(".rusage_maxrss_mb_impl", rusage_maxrss_mb_impl, envir = .GlobalEnv) + } + .rusage_maxrss_mb_impl() } # Print a formatted memory report to stderr via message(). From d232302c6145f6da32ee7d661ac0856a4ab42439 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 14:45:35 -0400 Subject: [PATCH 04/18] remove custom TMP summarization --- R/MSstatsSummarizeWithMultipleCores.R | 356 ++---------------- man/MSstatsSummarizeWithMultipleCores.Rd | 34 +- man/dot-MSstatsSummarizeSingleTMPV2.Rd | 48 --- ...kerV6.Rd => dot-build_summarize_worker.Rd} | 12 +- ...teinSlotV3.Rd => dot-pack_protein_slot.Rd} | 10 +- man/dot-reconstructProteinDTV3.Rd | 21 -- man/dot-unpack_protein_slot.Rd | 21 ++ 7 files changed, 80 insertions(+), 422 deletions(-) delete mode 100644 man/dot-MSstatsSummarizeSingleTMPV2.Rd rename man/{dot-buildSummarizeWorkerV6.Rd => dot-build_summarize_worker.Rd} (78%) rename man/{dot-buildProteinSlotV3.Rd => dot-pack_protein_slot.Rd} (64%) delete mode 100644 man/dot-reconstructProteinDTV3.Rd create mode 100644 man/dot-unpack_protein_slot.Rd diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 38dfae1f..08e3d2a3 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -312,290 +312,6 @@ -## ── Matrix-native TMP ───────────────────────────────────────────────────────── -## -## .summarize_protein_tmp_from_packed takes the packed double vector + meta list -## directly and produces TMP results without ever building the full long-format -## data.table or calling dcast. -## -## Compared with MSstatsSummarizeSingleTMP the hot path is: -## long-format: packed → long-DT (melt) → filter → AFT fit → dcast → TMP -## packed: packed → FL×R matrix → mask → (AFT fit if needed) → t(mat) → TMP -## -## The dcast savings are most significant for high-feature proteins (>100 -## features) and large run counts (>50 runs), where the intermediate FL*R -## data.table allocation dominates. -## ───────────────────────────────────────────────────────────────────────────── - - -#' Summarize a single protein with TMP directly from a packed double vector -#' -#' Bypasses the long-format \code{data.table} reconstruction and the -#' \code{dcast} inside \code{.fitTukey}. The FL×R packed matrices are -#' operated on directly: -#' \enumerate{ -#' \item Row/column masking replaces the \code{n_obs}/\code{n_obs_run} filter. -#' \item AFT survival fitting uses a lazily-built minimal \code{data.table} -#' (constructed only when \code{impute=TRUE} and censored values are -#' present). -#' \item TMP is applied via a vectorized scatter into a -#' \code{(LABEL×RUN) × FEATURE} wide matrix followed by -#' \code{median_polish_summary} — no \code{dcast} round-trip. -#' } -#' -#' @param packed double vector produced by \code{.pack_protein_slot} -#' @param meta metadata list from \code{.pack_protein_slot} -#' @param impute logical; impute censored values with AFT survival model -#' @param censored_symbol \code{"0"}, \code{"NA"}, or \code{NULL} -#' @param remove50missing logical; skip proteins where all runs are >50\% -#' missing -#' @param aft_iterations integer; max AFT iterations -#' @return \code{list(result_dt, survival_dt)} matching the format of -#' \code{MSstatsSummarizeSingleTMP} -#' @keywords internal -.summarize_protein_tmp_from_packed <- function(packed, meta, - impute, censored_symbol, remove50missing, aft_iterations = 90L) -{ - n_feature_labels <- meta$n_feature_labels - n_runs <- meta$n_runs - matrix_len <- n_feature_labels * n_runs - PROTEIN <- meta$PROTEIN - feature_label_dt <- meta$feature_label_dt # data.frame: FEATURE, LABEL, PEPTIDE - runs <- meta$runs - - # ── 1. Unpack matrices ───────────────────────────────────────────────────── - # Layout (see top-of-file comment block): - # pos 1 : slot_index (skip) - # pos 2..FL*R+1 : newABUNDANCE (FL × R, column-major) - # pos FL*R+2.. : ABUNDANCE (skip — not needed for TMP) - # pos 2*FL*R+2.. : censored (FL × R) - # pos 3*FL*R+2.. : cen (FL × R) - # pos 4*FL*R+2.. : ANOMALYSCORES(skip — not needed for TMP) - # then: n_obs(FL), n_obs_run(R), prop_features(R) - cursor <- 2L - new_abundance_mat <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs); cursor <- cursor + matrix_len - cursor <- cursor + matrix_len # skip ABUNDANCE - censored_mat <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs); cursor <- cursor + matrix_len - event_mat <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs); cursor <- cursor + matrix_len - cursor <- cursor + matrix_len # skip ANOMALYSCORES - n_obs_vec <- packed[cursor:(cursor + n_feature_labels - 1L)]; cursor <- cursor + n_feature_labels - n_obs_run_vec <- packed[cursor:(cursor + n_runs - 1L)]; cursor <- cursor + n_runs - prop_features_vec <- packed[cursor:(cursor + n_runs - 1L)] - - # ── 2. Row/column masking (mirrors the n_obs / n_obs_run filter) ─────────── - feature_is_kept <- !is.na(n_obs_vec) & n_obs_vec > 1 - run_is_kept <- !is.na(n_obs_run_vec) & n_obs_run_vec > 0 - - new_abundance_mat <- new_abundance_mat[feature_is_kept, run_is_kept, drop = FALSE] - censored_mat <- censored_mat [feature_is_kept, run_is_kept, drop = FALSE] - event_mat <- event_mat [feature_is_kept, run_is_kept, drop = FALSE] - - feature_label_dt_kept <- feature_label_dt[feature_is_kept, , drop = FALSE] - runs_kept <- runs[run_is_kept] - prop_features_kept <- prop_features_vec[run_is_kept] - n_feature_labels_kept <- nrow(new_abundance_mat) - n_runs_kept <- ncol(new_abundance_mat) - - # ── 3. Empty-matrix early exit ───────────────────────────────────────────── - if (n_feature_labels_kept == 0L || n_runs_kept == 0L) { - msg <- paste("Can't summarize for protein", PROTEIN, - "because all measurements are missing or censored.") - try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) - try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) - return(list(NULL, NULL)) - } - - # ── 4. remove50missing check (independent of imputation) ────────────────── - if (remove50missing && - all(prop_features_kept <= 0.5 | is.na(prop_features_kept))) { - msg <- paste("Can't summarize for protein", PROTEIN, - "because all runs have more than 50% missing values and", - "are removed with the option, remove50missing=TRUE.") - try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) - try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) - return(list(NULL, NULL)) - } - - # ── 5. AFT survival imputation (lazy: only if any censored present) ──────── - n_kept_rows <- n_feature_labels_kept * n_runs_kept - original_abundance_vec <- as.vector(new_abundance_mat) - - survival_columns <- c("newABUNDANCE", "cen", "RUN", "FEATURE", "ref_covariate") - any_censored <- meta$has_censored && any(censored_mat > 0.5, na.rm = TRUE) - - if (impute && any_censored && meta$has_cen) { - # Build minimal long-format DT — only columns needed by .fitSurvival - survival_dt <- data.table::data.table( - newABUNDANCE = original_abundance_vec, - cen = as.vector(event_mat), - censored = as.logical(as.vector(censored_mat) > 0.5), - FEATURE = factor(rep(feature_label_dt_kept$FEATURE, times = n_runs_kept)), - RUN = factor(rep(runs_kept, each = n_feature_labels_kept)), - LABEL = rep(as.character(feature_label_dt_kept$LABEL), times = n_runs_kept) - ) - if (meta$add_ref_covariate) { - survival_dt[, ref_covariate := factor( - data.table::fifelse(LABEL == "L", as.character(RUN), "0"))] - } - - fit_columns <- intersect(survival_columns, colnames(survival_dt)) - fit_data <- if (meta$is_labeled_ref) { - survival_dt[LABEL != "H", fit_columns, with = FALSE] - } else { - survival_dt[, fit_columns, with = FALSE] - } - - converged <- TRUE - survival_fit <- withCallingHandlers({ - MSstats:::.fitSurvival(fit_data, aft_iterations) - }, warning = function(w) { - if (grepl("converge", conditionMessage(w), ignore.case = TRUE)) { - message("Convergence warning caught: ", conditionMessage(w)) - converged <<- FALSE - } - }) - - predicted_values <- if (converged) { - predict(survival_fit, newdata = survival_dt) - } else { - rep(NA_real_, n_kept_rows) - } - - should_impute <- if (meta$is_labeled_ref) { - survival_dt$censored & survival_dt$LABEL != "H" - } else { - survival_dt$censored - } - - imputed_abundance_vec <- ifelse(should_impute, predicted_values, original_abundance_vec) - survival_dt[, predicted := ifelse(should_impute, predicted_values, NA_real_)] - survival_dt[, newABUNDANCE := imputed_abundance_vec] - new_abundance_mat <- matrix(imputed_abundance_vec, nrow = n_feature_labels_kept, ncol = n_runs_kept) - survival_output <- survival_dt[, intersect(c(survival_columns, "LABEL", "predicted"), - colnames(survival_dt)), with = FALSE] - - } else { - # No imputation: build survival DT from original values - survival_dt <- data.table::data.table( - newABUNDANCE = original_abundance_vec, - cen = if (meta$has_cen) as.vector(event_mat) - else rep(NA_real_, n_kept_rows), - FEATURE = rep(as.character(feature_label_dt_kept$FEATURE), times = n_runs_kept), - RUN = rep(runs_kept, each = n_feature_labels_kept), - LABEL = rep(as.character(feature_label_dt_kept$LABEL), times = n_runs_kept), - predicted = NA # logical NA, matching MSstatsSummarizeSingleTMP's - # bare `survival[, predicted := NA]` in its - # equivalent non-imputed branch (dataProcess.R:586) - ) - survival_output <- survival_dt[, intersect(c(survival_columns, "LABEL", "predicted"), - colnames(survival_dt)), with = FALSE] - } - - # ── 6. Post-imputation .isSummarizable check ─────────────────────────────── - if (all(is.na(new_abundance_mat) | new_abundance_mat == 0)) { - msg <- paste("Can't summarize for protein", PROTEIN, - "because all measurements are missing or censored.") - try(getOption("MSstatsMsg")("INFO", msg), silent = TRUE) - try(getOption("MSstatsLog")("INFO", msg), silent = TRUE) - return(list(NULL, NULL)) - } - - # ── 7. TMP via vectorized scatter into (LABEL×RUN) × FEATURE wide matrix ─── - # - # Mirrors .fitTukey: dcast(LABEL + RUN ~ FEATURE, value.var="newABUNDANCE") - # Row order in wide_mat: (lbl1,run1),(lbl1,run2),...,(lbl2,run1),... - # where labels are sorted — matches dcast alphabetical ordering for - # character LABEL values. - if (n_feature_labels_kept == 1L) { - # ── Single (FEATURE, LABEL) row: skip TMP, use values directly ──────── - if (meta$is_labeled_ref) { - is_heavy_row <- as.character(feature_label_dt_kept$LABEL) == "H" - is_light_row <- as.character(feature_label_dt_kept$LABEL) == "L" - if (any(is_heavy_row) && any(is_light_row)) { - heavy_values <- as.vector(new_abundance_mat[is_heavy_row, , drop = FALSE]) - light_values <- as.vector(new_abundance_mat[is_light_row, , drop = FALSE]) - heavy_median <- stats::median(heavy_values, na.rm = TRUE) - result <- data.table::data.table( - Protein = PROTEIN, - LABEL = "L", - RUN = runs_kept, - LogIntensities = light_values - heavy_values + heavy_median - ) - } else { - result <- data.table::data.table( - Protein = PROTEIN, - LABEL = as.character(feature_label_dt_kept$LABEL[1L]), - RUN = runs_kept, - LogIntensities = as.vector(new_abundance_mat) - ) - } - } else { - result <- data.table::data.table( - Protein = PROTEIN, - LABEL = rep(as.character(feature_label_dt_kept$LABEL), times = n_runs_kept), - RUN = rep(runs_kept, each = n_feature_labels_kept), - LogIntensities = as.vector(new_abundance_mat) - ) - } - } else { - # ── Multi-feature: scatter into wide matrix, apply TMP ───────────────── - unique_labels <- sort(unique(as.character(feature_label_dt_kept$LABEL))) - unique_features <- sort(unique(as.character(feature_label_dt_kept$FEATURE))) - n_labels <- length(unique_labels) - n_features <- length(unique_features) - - label_idx_map <- match(as.character(feature_label_dt_kept$LABEL), unique_labels) - feature_idx_map <- match(as.character(feature_label_dt_kept$FEATURE), unique_features) - - # Vectorized scatter: mat[fl, run] → wide[(lbl_idx-1)*n_runs_kept + run, feat_idx] - feature_idx_of_cell <- rep(seq_len(n_feature_labels_kept), times = n_runs_kept) - run_idx_of_cell <- rep(seq_len(n_runs_kept), each = n_feature_labels_kept) - wide_mat_row <- (label_idx_map[feature_idx_of_cell] - 1L) * n_runs_kept + run_idx_of_cell - wide_mat_col <- feature_idx_map[feature_idx_of_cell] - - wide_mat <- matrix(NA_real_, nrow = n_labels * n_runs_kept, ncol = n_features) - wide_mat[cbind(wide_mat_row, wide_mat_col)] <- as.vector(new_abundance_mat) - - tmp_fitted_values <- MSstats:::median_polish_summary(wide_mat) - - # Row → (LABEL, RUN) mapping: label index cycles every n_runs_kept rows - result_labels <- rep(unique_labels, each = n_runs_kept) - result_runs <- rep(runs_kept, times = n_labels) - - if (meta$is_labeled_ref) { - heavy_label_idx <- match("H", unique_labels) - light_label_idx <- match("L", unique_labels) - if (!is.na(heavy_label_idx) && !is.na(light_label_idx)) { - heavy_rows <- (heavy_label_idx - 1L) * n_runs_kept + seq_len(n_runs_kept) - light_rows <- (light_label_idx - 1L) * n_runs_kept + seq_len(n_runs_kept) - heavy_median <- stats::median(tmp_fitted_values[heavy_rows], na.rm = TRUE) - result <- data.table::data.table( - Protein = PROTEIN, - LABEL = "L", - RUN = runs_kept, - LogIntensities = tmp_fitted_values[light_rows] - tmp_fitted_values[heavy_rows] + heavy_median - ) - } else { - result <- data.table::data.table( - Protein = PROTEIN, - LABEL = result_labels, - RUN = result_runs, - LogIntensities = tmp_fitted_values - ) - } - } else { - result <- data.table::data.table( - Protein = PROTEIN, - LABEL = result_labels, - RUN = result_runs, - LogIntensities = tmp_fitted_values - ) - } - } - - list(result, survival_output) -} - #' Build the per-record worker closure for #' \code{MSstatsSummarizeWithMultipleCores} #' @@ -627,58 +343,46 @@ equal_variance_ <- equal_variance function(record) { - packed <- record$packed meta <- record$meta + protein_dt <- unpack_fn(record$packed, meta) result <- if (use_TMP_) { - .summarize_protein_tmp_from_packed( - packed, meta, - impute_, censored_symbol_, + MSstatsSummarizeSingleTMP( + protein_dt, impute_, censored_symbol_, remove50missing_, aft_iterations_) } else { - protein_dt <- unpack_fn(packed, meta) MSstatsSummarizeSingleLinear( protein_dt, impute_, censored_symbol_, remove50missing_, aft_iterations_, equal_variances = equal_variance_) } - # Normalize column types/levels at the worker boundary rather than - # inside .summarize_protein_tmp_from_packed/MSstatsSummarizeSingleLinear: - # those are shared with earlier packed-vector iterations and - # dataProcess.R respectively, and the TMP path in particular carries - # RUN as plain character (from meta$runs) and cen as double (from the - # packed matrix). Re-leveling against meta$runs (rather than a bare - # factor(RUN)) preserves the numeric run order that - # levels(input$RUN) already established upstream — a bare factor() - # call would instead re-sort alphabetically ("1","10","11","2",... - # for >=10 runs). as.character() first strips any existing level - # order before re-leveling, since - # .summarize_protein_tmp_from_packed's survival table (result[[2L]]) - # already wraps RUN in a bare factor() for the imputed branch — that - # call has the exact same alphabetical-resort bug even though it - # looks pre-typed. droplevels() then matches - # MSstatsSummarizeSingleTMP/SingleCore, whose factor(RUN) only ever - # sees — and so only ever keeps — the runs present for that protein, - # since meta$runs carries every run across the whole input. - # - # FEATURE gets the same treatment: .summarize_protein_tmp_from_packed's - # non-imputed survival branch leaves FEATURE as plain character (only - # the imputed branch wraps it in a bare factor()), so result[[2L]]$ - # FEATURE isn't reliably a factor the way SingleCore's is (single_ - # protein[, FEATURE := factor(FEATURE)] before survival is sliced off - # it). + # Normalize column types/levels at the worker boundary. Both + # MSstatsSummarizeSingleTMP and MSstatsSummarizeSingleLinear call a + # bare factor(RUN)/factor(FEATURE) on their input. For + # MSstatsSummarizeWithSingleCore that's a no-op order-wise, since RUN + # and FEATURE arrive already factors there and bare factor() keeps + # existing levels. Here protein_dt comes from .unpack_protein_slot, + # which reconstructs RUN/FEATURE as plain character — so the bare + # factor() call inside those functions re-sorts alphabetically + # instead ("1","10","11","2",... for >=10 runs). Re-leveling against + # meta$runs restores the numeric run order that levels(input$RUN) + # established upstream; as.character() first strips the alphabetical + # level order factor() just introduced. droplevels() then matches + # SingleCore's factor(RUN), which only ever sees — and so only ever + # keeps — the runs present for that protein, since meta$runs carries + # every run across the whole input. # - # Deliberately NOT leveling against meta$feature_label_dt$FEATURE's - # existing row order here: that order comes from .pack_protein_slot's - # data.table::setorder(feature_label_dt, FEATURE, LABEL), and - # data.table sorts character columns in the C-locale (byte order: - # all uppercase before any lowercase) for platform-independence — - # whereas SingleCore's bare factor(FEATURE) goes through base R's - # factor()/sort(), which use the session's collation locale (e.g. - # en_US.UTF-8, where case is interleaved: "a" < "A" < "b" < "B"). Any - # FEATURE strings with mixed case (e.g. modification tags) sort - # differently between the two, so re-deriving the level order with a - # base R sort() reproduces SingleCore's order instead of inheriting - # data.table's. + # FEATURE is re-leveled with a fresh base R sort() rather than + # inheriting meta$feature_label_dt's existing row order: that order + # comes from .pack_protein_slot's data.table::setorder(feature_label_dt, + # FEATURE, LABEL), and data.table sorts character columns in the + # C-locale (byte order: all uppercase before any lowercase) for + # platform-independence — whereas SingleCore's bare factor(FEATURE) + # goes through base R's factor()/sort(), which use the session's + # collation locale (e.g. en_US.UTF-8, where case is interleaved: + # "a" < "A" < "b" < "B"). Any FEATURE strings with mixed case (e.g. + # modification tags) sort differently between the two, so + # re-deriving the level order with a base R sort() reproduces + # SingleCore's order instead of inheriting data.table's. feature_levels <- sort(unique(meta$feature_label_dt$FEATURE)) for (idx in 1:2) { if (!is.null(result[[idx]])) { diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index dcf9ee4c..ebdeb629 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R \name{MSstatsSummarizeWithMultipleCores} \alias{MSstatsSummarizeWithMultipleCores} -\title{Feature-level data summarization via socket-dispatched protein records (V6)} +\title{Feature-level data summarization via socket-dispatched protein records} \usage{ MSstatsSummarizeWithMultipleCores( input, @@ -63,25 +63,27 @@ track progress. Only works for Linux & Mac OS. Default is 1.} workers.} } \value{ -A named list with one element per protein slot, identical in - structure to \code{MSstatsSummarizeWithMultipleCores}. +A named list with one element per protein slot, keyed by protein + (or protein \eqn{\times} label) identifier. } \description{ -Fixes a hidden RAM cost in \code{MSstatsSummarizeWithMultipleCoresV5}: V5's -worker closure captured \code{meta_list} — the metadata for \emph{every} -protein — by reference, so each worker received and retained the full -metadata set for the whole run regardless of how many proteins were -actually assigned to it. Bounding the per-task \emph{packed-vector} payload -(e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce that, -since the metadata was baked into the closure once, not re-sliced per task. +Fixes a hidden RAM cost present in an earlier iteration of this function, +where the worker closure captured \code{meta_list} — the metadata for +\emph{every} protein — by reference, so each worker received and retained +the full metadata set for the whole run regardless of how many proteins +were actually assigned to it. Bounding the per-task \emph{packed-vector} +payload (e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce +that, since the metadata was baked into the closure once, not re-sliced +per task. } \details{ -V6 pairs each protein's packed double vector with only its own metadata -into a single list element (\code{list(packed = ..., meta = ...)}) and -dispatches that combined list through \code{bplapply}. \code{BiocParallel} -then only serializes and sends the records actually assigned to a given -task, so a worker never holds metadata for proteins outside its own -task(s) — unlike V5, RAM scales with the batch actually being processed. +This version pairs each protein's packed double vector with only its own +metadata into a single list element (\code{list(packed = ..., meta = ...)}) +and dispatches that combined list through \code{bplapply}. +\code{BiocParallel} then only serializes and sends the records actually +assigned to a given task, so a worker never holds metadata for proteins +outside its own task(s) — unlike that earlier iteration, RAM scales with +the batch actually being processed. Progress is reported by turning on \code{SnowfastParam}'s built-in \code{progressbar}, not by having workers talk back to the parent. diff --git a/man/dot-MSstatsSummarizeSingleTMPV2.Rd b/man/dot-MSstatsSummarizeSingleTMPV2.Rd deleted file mode 100644 index fb1f9060..00000000 --- a/man/dot-MSstatsSummarizeSingleTMPV2.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.MSstatsSummarizeSingleTMPV2} -\alias{.MSstatsSummarizeSingleTMPV2} -\title{Summarize a single protein with TMP directly from a V3 packed double vector} -\usage{ -.MSstatsSummarizeSingleTMPV2( - packed, - meta, - impute, - censored_symbol, - remove50missing, - aft_iterations = 90L -) -} -\arguments{ -\item{packed}{double vector produced by \code{.buildProteinSlotV3}} - -\item{meta}{metadata list from \code{.buildProteinSlotV3}} - -\item{impute}{logical; impute censored values with AFT survival model} - -\item{censored_symbol}{\code{"0"}, \code{"NA"}, or \code{NULL}} - -\item{remove50missing}{logical; skip proteins where all runs are >50\% -missing} - -\item{aft_iterations}{integer; max AFT iterations} -} -\value{ -\code{list(result_dt, survival_dt)} matching the format of - \code{MSstatsSummarizeSingleTMP} -} -\description{ -Bypasses the long-format \code{data.table} reconstruction and the -\code{dcast} inside \code{.fitTukey}. The FL×R packed matrices are -operated on directly: -\enumerate{ - \item Row/column masking replaces the \code{n_obs}/\code{n_obs_run} filter. - \item AFT survival fitting uses a lazily-built minimal \code{data.table} - (constructed only when \code{impute=TRUE} and censored values are - present). - \item TMP is applied via a vectorized scatter into a - \code{(LABEL×RUN) × FEATURE} wide matrix followed by - \code{median_polish_summary} — no \code{dcast} round-trip. -} -} -\keyword{internal} diff --git a/man/dot-buildSummarizeWorkerV6.Rd b/man/dot-build_summarize_worker.Rd similarity index 78% rename from man/dot-buildSummarizeWorkerV6.Rd rename to man/dot-build_summarize_worker.Rd index acd91a5e..8970d6a9 100644 --- a/man/dot-buildSummarizeWorkerV6.Rd +++ b/man/dot-build_summarize_worker.Rd @@ -1,11 +1,11 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.buildSummarizeWorkerV6} -\alias{.buildSummarizeWorkerV6} +\name{.build_summarize_worker} +\alias{.build_summarize_worker} \title{Build the per-record worker closure for -\code{MSstatsSummarizeWithMultipleCoresV6}} +\code{MSstatsSummarizeWithMultipleCores}} \usage{ -.buildSummarizeWorkerV6( +.build_summarize_worker( use_TMP, impute, censored_symbol, @@ -16,13 +16,13 @@ } \description{ Defined at package top level — not nested via \code{local()} inside -\code{MSstatsSummarizeWithMultipleCoresV6} — so the returned closure's +\code{MSstatsSummarizeWithMultipleCores} — so the returned closure's enclosing environment chain is this factory's own (small) call frame plus the package namespace. \code{BiocParallel}/\code{matter} serialize a closure's entire enclosing environment chain to ship it to each socket worker, not just the variables the closure body actually references. A closure built via \code{local()} inside -\code{MSstatsSummarizeWithMultipleCoresV6} would have that function's own +\code{MSstatsSummarizeWithMultipleCores} would have that function's own evaluation frame in its chain — which holds \code{input}, \code{protein_records}, and other run-scale objects — so every worker would receive and retain a serialized copy of them even though diff --git a/man/dot-buildProteinSlotV3.Rd b/man/dot-pack_protein_slot.Rd similarity index 64% rename from man/dot-buildProteinSlotV3.Rd rename to man/dot-pack_protein_slot.Rd index 2b576d70..70849c0e 100644 --- a/man/dot-buildProteinSlotV3.Rd +++ b/man/dot-pack_protein_slot.Rd @@ -1,15 +1,15 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.buildProteinSlotV3} -\alias{.buildProteinSlotV3} +\name{.pack_protein_slot} +\alias{.pack_protein_slot} \title{Build the packed double vector and lightweight metadata for one protein slot} \usage{ -.buildProteinSlotV3(dt, slot_k, all_runs) +.pack_protein_slot(protein_dt, slot_index, all_runs) } \arguments{ -\item{dt}{data.table rows belonging to one protein (or protein × label) slot} +\item{protein_dt}{data.table rows belonging to one protein (or protein × label) slot} -\item{slot_k}{integer position of this slot in the global protein list} +\item{slot_index}{integer position of this slot in the global protein list} \item{all_runs}{character vector of all run names in global order} } diff --git a/man/dot-reconstructProteinDTV3.Rd b/man/dot-reconstructProteinDTV3.Rd deleted file mode 100644 index 099a6f28..00000000 --- a/man/dot-reconstructProteinDTV3.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.reconstructProteinDTV3} -\alias{.reconstructProteinDTV3} -\title{Reconstruct a per-protein data.table from a V3 packed double vector} -\usage{ -.reconstructProteinDTV3(packed, meta) -} -\arguments{ -\item{packed}{double vector produced by \code{.buildProteinSlotV3}} - -\item{meta}{metadata list from \code{.buildProteinSlotV3}} -} -\value{ -data.table compatible with \code{MSstatsSummarizeSingleTMP} / - \code{MSstatsSummarizeSingleLinear} -} -\description{ -Reconstruct a per-protein data.table from a V3 packed double vector -} -\keyword{internal} diff --git a/man/dot-unpack_protein_slot.Rd b/man/dot-unpack_protein_slot.Rd new file mode 100644 index 00000000..d0a11738 --- /dev/null +++ b/man/dot-unpack_protein_slot.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.unpack_protein_slot} +\alias{.unpack_protein_slot} +\title{Reconstruct a per-protein data.table from a packed double vector} +\usage{ +.unpack_protein_slot(packed, meta) +} +\arguments{ +\item{packed}{double vector produced by \code{.pack_protein_slot}} + +\item{meta}{metadata list from \code{.pack_protein_slot}} +} +\value{ +data.table compatible with \code{MSstatsSummarizeSingleTMP} / + \code{MSstatsSummarizeSingleLinear} +} +\description{ +Reconstruct a per-protein data.table from a packed double vector +} +\keyword{internal} From 9d977259849dbaf7050aa3768e7c9fcabb1cae2d Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 14:47:20 -0400 Subject: [PATCH 05/18] remove unnecessary current memory rss function --- R/MSstatsSummarizeWithMultipleCores.R | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 08e3d2a3..fd099f8d 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -1,18 +1,3 @@ -## ── Memory-monitoring helpers ───────────────────────────────────────────────── - -# RSS of the current process in MB. -# On Linux reads /proc/self/status (VmRSS); elsewhere falls back to gc() counts. -.current_rss_mb <- function() { - if (file.exists("/proc/self/status")) { - ln <- readLines("/proc/self/status", warn = FALSE) - m <- grep("^VmRSS:", ln, value = TRUE) - if (length(m)) - return(as.numeric(gsub("[^0-9]", "", m[1L])) / 1024) - } - g <- gc(reset = FALSE) - (g["Ncells", "used"] * 8L + g["Vcells", "used"] * 8L) / 1024^2 -} - # Cross-platform peak-RSS reader. Reflects the true lifetime peak of the # calling process, regardless of when you call it — no polling required. # All three branches report OS-level peak resident/working-set memory, From 1deddffd2ace9e4820f2c4375a000c448b85df13 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 15:15:46 -0400 Subject: [PATCH 06/18] clean up comments --- R/MSstatsSummarizeWithMultipleCores.R | 268 +++-------------------- man/MSstatsSummarizeSingleTMP.Rd | 17 +- man/MSstatsSummarizeWithMultipleCores.Rd | 77 ++----- man/MSstatsSummarizeWithSingleCore.Rd | 24 +- man/dot-build_summarize_worker.Rd | 18 +- man/dot-getNonMissingFilterStats.Rd | 8 +- man/dot-pack_protein_slot.Rd | 6 +- man/dot-runTukey.Rd | 11 +- 8 files changed, 68 insertions(+), 361 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index fd099f8d..9c694777 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -1,7 +1,3 @@ -# Cross-platform peak-RSS reader. Reflects the true lifetime peak of the -# calling process, regardless of when you call it — no polling required. -# All three branches report OS-level peak resident/working-set memory, -# so figures are comparable across platforms. .peak_rss_mb <- function() { if (.Platform$OS.type == "windows") { if (!exists(".peakRSS_windows_impl", mode = "function")) { @@ -14,7 +10,7 @@ double peakRSS_windows_impl() { PROCESS_MEMORY_COUNTERS pmc; if (GetProcessMemoryInfo(GetCurrentProcess(), &pmc, sizeof(pmc))) { - return (double) pmc.PeakWorkingSetSize / (1024.0*1024.0); // bytes -> MB + return (double) pmc.PeakWorkingSetSize / (1024.0*1024.0); } return NA_REAL; }") @@ -22,15 +18,12 @@ } return(.peakRSS_windows_impl()) } - + if (file.exists("/proc/self/status")) { - # Linux: VmHWM = kernel-maintained peak resident set size ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024) } - - # macOS (also works as a Linux fallback): POSIX getrusage() ru_maxrss - # is likewise a lifetime high-water mark, just different units per OS. + if (!exists(".rusage_maxrss_mb_impl", mode = "function")) { Rcpp::cppFunction( depends = "Rcpp", @@ -39,9 +32,9 @@ double rusage_maxrss_mb_impl() { struct rusage ru; getrusage(RUSAGE_SELF, &ru); #ifdef __APPLE__ - return (double) ru.ru_maxrss / (1024.0*1024.0); // bytes -> MB + return (double) ru.ru_maxrss / (1024.0*1024.0); #else - return (double) ru.ru_maxrss / 1024.0; // KB -> MB + return (double) ru.ru_maxrss / 1024.0; #endif }") assign(".rusage_maxrss_mb_impl", rusage_maxrss_mb_impl, envir = .GlobalEnv) @@ -49,10 +42,6 @@ .rusage_maxrss_mb_impl() } -# Print a formatted memory report to stderr via message(). -# checkpoints: named numeric vector of RSS snapshots (MB). -# worker_peak_mb: numeric vector of per-worker peak RSS values (may be NA). -# elapsed: total wall-clock seconds (NULL to omit). .print_memory_report <- function(function_name, checkpoints, worker_peak_mb = NULL, elapsed = NULL) { rule_width <- 65L @@ -91,27 +80,9 @@ message(paste(lines, collapse = "\n")) } -## ── Protein-slot pack/unpack helpers ─────────────────────────────────────────── -## -## Packed double-vector layout for one protein slot (all column-major matrices). -## FL = n_feature_labels (unique FEATURE × LABEL pairs), R = n_runs: -## -## pos 1 : slot_index (protein index, cast to double) -## pos 2 .. FL*R+1 : newABUNDANCE (FL × R) -## pos FL*R+2 .. 2*FL*R+1 : ABUNDANCE (FL × R; NA for unlabeled/TMP) -## pos 2*FL*R+2 .. 3*FL*R+1 : censored (FL × R; 0.0/1.0) -## pos 3*FL*R+2 .. 4*FL*R+1 : cen (FL × R; 1-censored event flag) -## pos 4*FL*R+2 .. 5*FL*R+1 : ANOMALYSCORES (FL × R; NA if unused) -## pos 5*FL*R+2 .. 5*FL*R+FL+1 : n_obs (FL; per feature-label) -## pos 5*FL*R+FL+2.. 5*FL*R+FL+R+1: n_obs_run (R; per run) -## pos 5*FL*R+FL+R+2..5*FL*R+FL+2R+1: prop_features (R; per run) -## -## Total length: 1 + 5*FL*R + FL + 2*R -## ───────────────────────────────────────────────────────────────────────────── - -#' Build the packed double vector and lightweight metadata for one protein slot +#' Pack one protein slot into a double vector plus metadata #' -#' @param protein_dt data.table rows belonging to one protein (or protein × label) slot +#' @param protein_dt data.table rows for one protein (or protein x label) slot #' @param slot_index integer position of this slot in the global protein list #' @param all_runs character vector of all run names in global order #' @return list with elements \code{packed} (double vector) and \code{meta} (list) @@ -120,7 +91,6 @@ n_runs <- length(all_runs) - # ── Unique (FEATURE, LABEL) → "effective features", stable ordering ─────── has_peptide <- "PEPTIDE" %in% colnames(protein_dt) feature_label_dt <- unique(protein_dt[, .( FEATURE = as.character(FEATURE), @@ -130,7 +100,6 @@ data.table::setorder(feature_label_dt, FEATURE, LABEL) n_feature_labels <- nrow(feature_label_dt) - # ── Index maps: feature-label → row index, run name → column index ───────── feature_row_idx <- match( paste(as.character(protein_dt$FEATURE), as.character(protein_dt$LABEL), sep = "\t"), paste(feature_label_dt$FEATURE, feature_label_dt$LABEL, sep = "\t")) @@ -138,7 +107,6 @@ run_col_idx <- run_col_idx_lookup[as.character(protein_dt$RUN)] row_is_valid <- !is.na(feature_row_idx) & !is.na(run_col_idx) - # ── Helper: allocate n_feature_labels × n_runs matrix, fill from a protein_dt column ── scatter_into_matrix <- function(col_vec, fill = NA_real_) { m <- matrix(fill, nrow = n_feature_labels, ncol = n_runs) m[cbind(feature_row_idx[row_is_valid], run_col_idx[row_is_valid])] <- @@ -146,7 +114,6 @@ m } - # ── Build the five numeric matrices ─────────────────────────────────────── new_abundance_mat <- scatter_into_matrix(protein_dt$newABUNDANCE) has_ABUNDANCE <- "ABUNDANCE" %in% colnames(protein_dt) @@ -155,7 +122,7 @@ has_censored <- "censored" %in% colnames(protein_dt) censored_mat <- if (has_censored) scatter_into_matrix(as.double(protein_dt$censored)) else - matrix(0.0, n_feature_labels, n_runs) # treat as non-censored when column absent + matrix(0.0, n_feature_labels, n_runs) has_cen <- "cen" %in% colnames(protein_dt) event_mat <- if (has_cen) scatter_into_matrix(protein_dt$cen) else @@ -166,7 +133,6 @@ anomaly_scores_mat <- if (has_anom) scatter_into_matrix(protein_dt$ANOMALYSCORES) else matrix(NA_real_, n_feature_labels, n_runs) - # ── Per-feature scalar: n_obs (constant within FEATURE × LABEL) ─────────── n_obs_by_feature_label <- protein_dt[, .(n_obs = as.double(n_obs[1L])), by = .(FEATURE = as.character(FEATURE), LABEL = as.character(LABEL))] @@ -174,7 +140,6 @@ data.table::setorder(feature_label_with_nobs, FEATURE, LABEL) n_obs_vec <- feature_label_with_nobs$n_obs - # ── Per-run scalars: n_obs_run, prop_features (constant within RUN) ─────── run_level_scalars <- protein_dt[, .(n_obs_run = as.double(n_obs_run[1L]), prop_features = as.double(prop_features[1L])), by = .(RUN = as.character(RUN))] @@ -182,7 +147,6 @@ n_obs_run_vec <- run_scalars_all_runs$n_obs_run prop_features_vec <- run_scalars_all_runs$prop_features - # ── Pack ───────────────────────────────────────────────────────────────── packed <- c( as.double(slot_index), as.vector(new_abundance_mat), @@ -195,10 +159,9 @@ prop_features_vec ) - # ── Metadata (string labels, flags; kept in main-process RAM) ──────────── meta <- list( PROTEIN = as.character(protein_dt$PROTEIN[1L]), - feature_label_dt = as.data.frame(feature_label_dt), # n_feature_labels × 3: FEATURE, LABEL, PEPTIDE + feature_label_dt = as.data.frame(feature_label_dt), runs = all_runs, n_feature_labels = n_feature_labels, n_runs = n_runs, @@ -228,7 +191,6 @@ n_runs <- meta$n_runs matrix_len <- n_feature_labels * n_runs - # ── Unpack sections (1-indexed; position 1 is slot_index header) ───────── cursor <- 2L read_next_matrix <- function() { m <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs) @@ -245,8 +207,7 @@ n_obs_run_vec <- packed[cursor:(cursor + n_runs - 1L)]; cursor <- cursor + n_runs prop_features_vec <- packed[cursor:(cursor + n_runs - 1L)] - # ── Column-major melt: rep/each match the matrix column-major ordering ──── - feature_label_dt <- meta$feature_label_dt # data.frame: FEATURE, LABEL, PEPTIDE (n_feature_labels rows) + feature_label_dt <- meta$feature_label_dt runs <- meta$runs n_rows <- n_feature_labels * n_runs @@ -262,26 +223,21 @@ prop_features = rep(prop_features_vec, each = n_feature_labels) ) - # Optional: ABUNDANCE (labeled linear model uses raw log-intensities) if (meta$has_ABUNDANCE) { protein_dt[, ABUNDANCE := as.vector(abundance_mat)] } - # censored: stored as 0.0/1.0 doubles; NA → treat as non-censored protein_dt[, censored := { v <- as.vector(censored_mat) if (meta$has_censored) as.logical(v > 0.5) else rep(FALSE, n_rows) }] - # cen: survival event indicator (1 = observed, 0 = left-censored) if (meta$has_cen) { protein_dt[, cen := as.vector(event_mat)] } - # ANOMALYSCORES: NA when unused (linear model skips anomaly weighting) protein_dt[, ANOMALYSCORES := as.vector(anomaly_scores_mat)] - # SRM-specific columns derived from LABEL and RUN if (meta$is_labeled_ref) { protein_dt[, is_labeled_ref := (LABEL == "H")] if (meta$add_ref_covariate) { @@ -294,25 +250,10 @@ } - - - -#' Build the per-record worker closure for -#' \code{MSstatsSummarizeWithMultipleCores} +#' Build the per-record worker closure for \code{MSstatsSummarizeWithMultipleCores} #' -#' Defined at package top level — not nested via \code{local()} inside -#' \code{MSstatsSummarizeWithMultipleCores} — so the returned closure's -#' enclosing environment chain is this factory's own (small) call frame plus -#' the package namespace. \code{BiocParallel}/\code{matter} serialize a -#' closure's entire enclosing environment chain to ship it to each socket -#' worker, not just the variables the closure body actually references. A -#' closure built via \code{local()} inside -#' \code{MSstatsSummarizeWithMultipleCores} would have that function's own -#' evaluation frame in its chain — which holds \code{input}, -#' \code{protein_records}, and other run-scale objects — so every worker -#' would receive and retain a serialized copy of them even though -#' \code{.worker} never touches them. Building the closure here instead keeps -#' its captured state to only the scalar run parameters. +#' Defined at package top level so the closure only captures the scalar run +#' parameters, not the caller's run-scale objects. #' #' @keywords internal .build_summarize_worker <- function( @@ -340,34 +281,6 @@ remove50missing_, aft_iterations_, equal_variances = equal_variance_) } - # Normalize column types/levels at the worker boundary. Both - # MSstatsSummarizeSingleTMP and MSstatsSummarizeSingleLinear call a - # bare factor(RUN)/factor(FEATURE) on their input. For - # MSstatsSummarizeWithSingleCore that's a no-op order-wise, since RUN - # and FEATURE arrive already factors there and bare factor() keeps - # existing levels. Here protein_dt comes from .unpack_protein_slot, - # which reconstructs RUN/FEATURE as plain character — so the bare - # factor() call inside those functions re-sorts alphabetically - # instead ("1","10","11","2",... for >=10 runs). Re-leveling against - # meta$runs restores the numeric run order that levels(input$RUN) - # established upstream; as.character() first strips the alphabetical - # level order factor() just introduced. droplevels() then matches - # SingleCore's factor(RUN), which only ever sees — and so only ever - # keeps — the runs present for that protein, since meta$runs carries - # every run across the whole input. - # - # FEATURE is re-leveled with a fresh base R sort() rather than - # inheriting meta$feature_label_dt's existing row order: that order - # comes from .pack_protein_slot's data.table::setorder(feature_label_dt, - # FEATURE, LABEL), and data.table sorts character columns in the - # C-locale (byte order: all uppercase before any lowercase) for - # platform-independence — whereas SingleCore's bare factor(FEATURE) - # goes through base R's factor()/sort(), which use the session's - # collation locale (e.g. en_US.UTF-8, where case is interleaved: - # "a" < "A" < "b" < "B"). Any FEATURE strings with mixed case (e.g. - # modification tags) sort differently between the two, so - # re-deriving the level order with a base R sort() reproduces - # SingleCore's order instead of inheriting data.table's. feature_levels <- sort(unique(meta$feature_label_dt$FEATURE)) for (idx in 1:2) { if (!is.null(result[[idx]])) { @@ -385,25 +298,12 @@ #' Per-worker peak-RSS query task for \code{MSstatsSummarizeWithMultipleCores} #' -#' Dispatched once per worker (via \code{seq_len(bpnworkers(BPPARAM))}) after -#' the main summarization \code{bplapply} call, while the persistent workers -#' are still alive, so each worker reports its own true lifetime-peak RSS -#' rather than a snapshot taken mid-run. -#' #' @keywords internal #' @noRd .report_worker_peak <- function(i) { list(worker = i, pid = Sys.getpid(), peak_mb = .peak_rss_mb()) } -# Per-worker warm-up task for MSstatsSummarizeWithMultipleCores: loads -# MSstats once per persistent worker process and pins data.table to a single -# thread. Defined at top level rather than inline inside -# MSstatsSummarizeWithMultipleCores: even though its body never references -# `input`/`protein_records`, a closure defined inline there would still carry -# that function's evaluation frame in its enclosing environment chain, and -# BiocParallel would serialize that whole frame — including those run-scale -# objects — to every worker just to ship this no-op task. .warmup_worker <- function(i) { library(MSstats, quietly = TRUE, warn.conflicts = FALSE) data.table::setDTthreads(1) @@ -413,69 +313,19 @@ #' Feature-level data summarization via socket-dispatched protein records #' -#' Fixes a hidden RAM cost present in an earlier iteration of this function, -#' where the worker closure captured \code{meta_list} — the metadata for -#' \emph{every} protein — by reference, so each worker received and retained -#' the full metadata set for the whole run regardless of how many proteins -#' were actually assigned to it. Bounding the per-task \emph{packed-vector} -#' payload (e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce -#' that, since the metadata was baked into the closure once, not re-sliced -#' per task. -#' -#' This version pairs each protein's packed double vector with only its own -#' metadata into a single list element (\code{list(packed = ..., meta = ...)}) -#' and dispatches that combined list through \code{bplapply}. -#' \code{BiocParallel} then only serializes and sends the records actually -#' assigned to a given task, so a worker never holds metadata for proteins -#' outside its own task(s) — unlike that earlier iteration, RAM scales with -#' the batch actually being processed. -#' -#' Progress is reported by turning on \code{SnowfastParam}'s built-in -#' \code{progressbar}, not by having workers talk back to the parent. -#' \code{BiocParallel} already ticks that progress bar from inside the -#' manager process itself, once per task result it collects — the same -#' receive that would happen regardless, over the same single -#' \code{bplapply} call. Workers never see the flag and send nothing extra -#' because of it, so this adds no serialization and no IPC beyond what an -#' unmonitored run already does. Reporting granularity therefore tracks -#' however many tasks the run is already split into (see -#' \code{max_proteins_per_worker} below): with the default \code{tasks = 0} -#' that's one step per worker. -#' #' @param input feature-level data processed by dataProcess subfunctions -#' @param method summarization method: "linear" or "TMP" -#' @param equal_variance only for summaryMethod = "linear". Default is TRUE. -#' Logical variable for whether the model should account for heterogeneous variation -#' among intensities from different features. Default is TRUE, which assume equal -#' variance among intensities from features. FALSE means that we cannot assume -#' equal variance among intensities from features, then we will account for -#' heterogeneous variation from different features. -#' @param censored_symbol Indicates how censored missing values are encoded in -#' the 'Intensity' column. 'NA' (default) treats all NA intensities as -#' left-censored (i.e., below the limit of detection). '0' treats zero -#' intensities as left-censored; in this case NA intensities are assumed to be -#' missing at random and are not censored. Skyline output should use '0'. NULL -#' assumes that all missing values are missing at random — no values are treated -#' as censored, and imputation is disabled (impute is ignored). -#' @param remove50missing only for summaryMethod = "TMP". TRUE removes the proteins -#' where every run has at least 50\% missing values for each peptide. FALSE is default. -#' @param impute only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'. -#' TRUE (default) imputes censored missing values using an Accelerated Failure -#' Time model. FALSE excludes censored observations from summarization entirely, -#' treating them as missing at random; no imputed values are introduced. -#' Has no effect when censored_symbol = NULL, since no values are considered censored. -#' @param numberOfCores Number of cores for parallel processing. When > 1, -#' a logfile named `MSstats_dataProcess_log_progress.log` is created to -#' track progress. Only works for Linux & Mac OS. Default is 1. -#' @param max_proteins_per_worker integer; caps how many protein records -#' (packed vector + its own metadata) are bundled into a single -#' \code{bplapply} task sent to one worker. Translated into -#' \code{SnowfastParam(tasks = ...)}: with \code{N} proteins this becomes -#' \code{tasks = ceiling(N / max_proteins_per_worker)}. Only applied when -#' \code{BPPARAM} is \code{NULL}; ignored if the caller supplies -#' \code{BPPARAM} directly. Default \code{0} reproduces \code{tasks = 0}: -#' \code{X} is divided as evenly as possible across \code{numberOfCores} -#' workers. +#' @param method summarization method: "linear" or "TMP" +#' @param impute only for method = "TMP"; imputes censored values via AFT model when TRUE +#' @param censored_symbol how censored values are encoded: 'NA', '0', or NULL for none +#' @param remove50missing only for method = "TMP"; drops proteins missing >=50% per peptide in every run +#' @param equal_variance only for method = "linear"; assume equal variance among feature intensities +#' @param numberOfCores number of cores for parallel processing (Linux/Mac only) +#' @param aft_iterations number of AFT model iterations +#' @param verbose whether to print verbose output +#' @param BPPARAM optional \code{BiocParallelParam} instance +#' @param track_memory whether to report per-worker peak RSS memory usage +#' @param max_proteins_per_worker caps protein records per \code{bplapply} task; +#' 0 uses BiocParallel's default split #' #' @return A named list with one element per protein slot, keyed by protein #' (or protein \eqn{\times} label) identifier. @@ -500,17 +350,15 @@ MSstatsSummarizeWithMultipleCores <- function( track_memory = FALSE, max_proteins_per_worker = 0L ) { - # ── 0. Single-core fallback ──────────────────────────────────────────────── if (numberOfCores <= 1L && is.null(BPPARAM)) { return(MSstatsSummarizeWithSingleCore( input, method, impute, censored_symbol, remove50missing, equal_variance, aft_iterations)) } - + start_time <- proc.time()[["elapsed"]] memory_checkpoints <- list() - - # ── 1. Split input by protein slot ──────────────────────────────────────── + is_labeled_reference <- "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE) split_keys <- if (is_labeled_reference) list(input$PROTEIN) else @@ -518,67 +366,33 @@ MSstatsSummarizeWithMultipleCores <- function( protein_indices <- split(seq_len(nrow(input)), split_keys) protein_ids <- names(protein_indices) num_proteins <- length(protein_indices) - - # Sort on input$RUN's own type (integer/numeric, typically) before - # converting to character. Sorting the character form directly — as a - # naive `sort(as.character(...))` would — collates lexicographically - # ("1","10","11","2",...) whenever RUN isn't already a pre-existing - # factor, e.g. when `input` comes straight from fread() and RUN reads in - # as integer. This mirrors what factor() itself does internally for a - # non-factor input: unique() + order() on the original values, then - # as.character() only at the end. + all_runs <- if (is.factor(input$RUN)) levels(input$RUN) else as.character(sort(unique(input$RUN))) - + getOption("MSstatsLog")("INFO", paste0("Packing ", num_proteins, " proteins × ", length(all_runs), " runs into per-protein records")) - # ── 2. Pack each protein into a (packed vector, own metadata) record ────── - # - # Unlike earlier iterations' parallel packed_list/meta_list arrays, each - # protein's metadata travels bundled with its own packed vector. - # bplapply/BPPARAM slice this combined list into tasks, so a worker's - # closure never needs — and never receives — metadata for proteins - # outside its assigned task(s). protein_records <- vector("list", num_proteins) for (slot_index in seq_len(num_proteins)) { protein_records[[slot_index]] <- .pack_protein_slot( input[protein_indices[[slot_index]], ], slot_index, all_runs) } - + payload_mb <- sum(vapply(protein_records, function(r) length(r$packed), integer(1))) * 8 / 1024^2 getOption("MSstatsLog")("INFO", paste0("Dispatching via sockets (", format(round(payload_mb, 1)), " MB total packed payload; metadata sharded per task)")) - - # ── 3. Worker closure ───────────────────────────────────────────────────── - # - # No captured meta_list: each task's records already carry their own - # metadata, so the closure only needs the scalar run parameters. - # - # Built via .build_summarize_worker() (defined at package top level) - # rather than local() here, so the closure's enclosing environment chain - # never includes this function's own frame — which holds `input`, - # `protein_records`, etc. — and BiocParallel doesn't serialize those - # run-scale objects to every worker. + use_TMP <- identical(method, "TMP") worker_fn <- .build_summarize_worker( use_TMP, impute, censored_symbol, remove50missing, aft_iterations, equal_variance) - - # ── 4. Dispatch ─────────────────────────────────────────────────────────── - # - # tasks controls how many protein records are bundled into one bplapply - # task (i.e. one message sent to one worker). tasks == 0 (default) leaves - # BiocParallel's own behavior in place: X divided as evenly as possible - # over numberOfCores workers. When max_proteins_per_worker > 0, tasks is - # sized so no task exceeds that many records — and because metadata now - # travels with each record instead of being fully captured by the - # closure, this actually bounds peak worker RAM. + if (is.null(BPPARAM)) { tasks <- if (max_proteins_per_worker > 0L) { as.integer(ceiling(num_proteins / max_proteins_per_worker)) @@ -597,30 +411,20 @@ MSstatsSummarizeWithMultipleCores <- function( force.GC = TRUE, stop.on.error = FALSE) } - - # ── Cluster setup ───────────────────────────────────────────────────────── - # Load MSstats once per persistent worker process instead of once per - # protein record, and pin data.table to a single thread per worker — - # otherwise each of the numberOfCores workers would independently - # auto-detect its own DT thread pool, oversubscribing the node's cores. + started_here <- !BiocParallel::bpisup(BPPARAM) if (started_here) { BiocParallel::bpstart(BPPARAM) on.exit(BiocParallel::bpstop(BPPARAM), add = TRUE) } - # Progress bar is on BPPARAM itself, so it would otherwise also print for - # these two trivial one-task-per-worker calls. Toggle it off around them - # and restore whatever it was (set above) so only the main summarization - # bplapply — the one call whose progress is actually informative — shows - # a bar. show_progress <- BiocParallel::bpprogressbar(BPPARAM) - + BiocParallel::bpprogressbar(BPPARAM) <- FALSE BiocParallel::bplapply( seq_len(BiocParallel::bpnworkers(BPPARAM)), .warmup_worker, BPPARAM = BPPARAM) BiocParallel::bpprogressbar(BPPARAM) <- show_progress - + results <- BiocParallel::bplapply(protein_records, worker_fn, BPPARAM = BPPARAM) names(results) <- protein_ids @@ -629,7 +433,7 @@ MSstatsSummarizeWithMultipleCores <- function( BiocParallel::bpprogressbar(BPPARAM) <- FALSE worker_peaks <- BiocParallel::bplapply( seq_len(BiocParallel::bpnworkers(BPPARAM)), - .report_worker_peak, BPPARAM = BPPARAM) # must run BEFORE bpstop() while workers are alive + .report_worker_peak, BPPARAM = BPPARAM) BiocParallel::bpprogressbar(BPPARAM) <- show_progress memory_checkpoints[["parent peak (main)"]] <- .peak_rss_mb() worker_peak_mb <- vapply(worker_peaks, function(x) x$peak_mb, numeric(1L)) diff --git a/man/MSstatsSummarizeSingleTMP.Rd b/man/MSstatsSummarizeSingleTMP.Rd index 041ef76e..26628f2e 100644 --- a/man/MSstatsSummarizeSingleTMP.Rd +++ b/man/MSstatsSummarizeSingleTMP.Rd @@ -15,22 +15,11 @@ MSstatsSummarizeSingleTMP( \arguments{ \item{single_protein}{feature-level data for a single protein} -\item{impute}{only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'. -TRUE (default) imputes censored missing values using an Accelerated Failure -Time model. FALSE excludes censored observations from summarization entirely, -treating them as missing at random; no imputed values are introduced. -Has no effect when censored_symbol = NULL, since no values are considered censored.} +\item{impute}{only for method = "TMP"; imputes censored values via AFT model when TRUE} -\item{censored_symbol}{Indicates how censored missing values are encoded in -the 'Intensity' column. 'NA' (default) treats all NA intensities as -left-censored (i.e., below the limit of detection). '0' treats zero -intensities as left-censored; in this case NA intensities are assumed to be -missing at random and are not censored. Skyline output should use '0'. NULL -assumes that all missing values are missing at random — no values are treated -as censored, and imputation is disabled (impute is ignored).} +\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins -where every run has at least 50\% missing values for each peptide. FALSE is default.} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} \item{aft_iterations}{number of iterations for AFT model fitting} } diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index ebdeb629..f77134c1 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -24,76 +24,31 @@ MSstatsSummarizeWithMultipleCores( \item{method}{summarization method: "linear" or "TMP"} -\item{impute}{only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'. -TRUE (default) imputes censored missing values using an Accelerated Failure -Time model. FALSE excludes censored observations from summarization entirely, -treating them as missing at random; no imputed values are introduced. -Has no effect when censored_symbol = NULL, since no values are considered censored.} +\item{impute}{only for method = "TMP"; imputes censored values via AFT model when TRUE} -\item{censored_symbol}{Indicates how censored missing values are encoded in -the 'Intensity' column. 'NA' (default) treats all NA intensities as -left-censored (i.e., below the limit of detection). '0' treats zero -intensities as left-censored; in this case NA intensities are assumed to be -missing at random and are not censored. Skyline output should use '0'. NULL -assumes that all missing values are missing at random — no values are treated -as censored, and imputation is disabled (impute is ignored).} +\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins -where every run has at least 50\% missing values for each peptide. FALSE is default.} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} -\item{equal_variance}{only for summaryMethod = "linear". Default is TRUE. -Logical variable for whether the model should account for heterogeneous variation -among intensities from different features. Default is TRUE, which assume equal -variance among intensities from features. FALSE means that we cannot assume -equal variance among intensities from features, then we will account for -heterogeneous variation from different features.} +\item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities} -\item{numberOfCores}{Number of cores for parallel processing. When > 1, -a logfile named `MSstats_dataProcess_log_progress.log` is created to -track progress. Only works for Linux & Mac OS. Default is 1.} +\item{numberOfCores}{number of cores for parallel processing (Linux/Mac only)} -\item{max_proteins_per_worker}{integer; caps how many protein records -(packed vector + its own metadata) are bundled into a single -\code{bplapply} task sent to one worker. Translated into -\code{SnowfastParam(tasks = ...)}: with \code{N} proteins this becomes -\code{tasks = ceiling(N / max_proteins_per_worker)}. Only applied when -\code{BPPARAM} is \code{NULL}; ignored if the caller supplies -\code{BPPARAM} directly. Default \code{0} reproduces \code{tasks = 0}: -\code{X} is divided as evenly as possible across \code{numberOfCores} -workers.} +\item{aft_iterations}{number of AFT model iterations} + +\item{verbose}{whether to print verbose output} + +\item{BPPARAM}{optional \code{BiocParallelParam} instance} + +\item{track_memory}{whether to report per-worker peak RSS memory usage} + +\item{max_proteins_per_worker}{caps protein records per \code{bplapply} task; +0 uses BiocParallel's default split} } \value{ A named list with one element per protein slot, keyed by protein (or protein \eqn{\times} label) identifier. } \description{ -Fixes a hidden RAM cost present in an earlier iteration of this function, -where the worker closure captured \code{meta_list} — the metadata for -\emph{every} protein — by reference, so each worker received and retained -the full metadata set for the whole run regardless of how many proteins -were actually assigned to it. Bounding the per-task \emph{packed-vector} -payload (e.g. via \code{SnowfastParam(tasks = ...)}) did nothing to reduce -that, since the metadata was baked into the closure once, not re-sliced -per task. -} -\details{ -This version pairs each protein's packed double vector with only its own -metadata into a single list element (\code{list(packed = ..., meta = ...)}) -and dispatches that combined list through \code{bplapply}. -\code{BiocParallel} then only serializes and sends the records actually -assigned to a given task, so a worker never holds metadata for proteins -outside its own task(s) — unlike that earlier iteration, RAM scales with -the batch actually being processed. - -Progress is reported by turning on \code{SnowfastParam}'s built-in -\code{progressbar}, not by having workers talk back to the parent. -\code{BiocParallel} already ticks that progress bar from inside the -manager process itself, once per task result it collects — the same -receive that would happen regardless, over the same single -\code{bplapply} call. Workers never see the flag and send nothing extra -because of it, so this adds no serialization and no IPC beyond what an -unmonitored run already does. Reporting granularity therefore tracks -however many tasks the run is already split into (see -\code{max_proteins_per_worker} below): with the default \code{tasks = 0} -that's one step per worker. +Feature-level data summarization via socket-dispatched protein records } diff --git a/man/MSstatsSummarizeWithSingleCore.Rd b/man/MSstatsSummarizeWithSingleCore.Rd index 4892dccb..2c3f23ef 100644 --- a/man/MSstatsSummarizeWithSingleCore.Rd +++ b/man/MSstatsSummarizeWithSingleCore.Rd @@ -19,29 +19,13 @@ MSstatsSummarizeWithSingleCore( \item{method}{summarization method: "linear" or "TMP"} -\item{impute}{only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'. -TRUE (default) imputes censored missing values using an Accelerated Failure -Time model. FALSE excludes censored observations from summarization entirely, -treating them as missing at random; no imputed values are introduced. -Has no effect when censored_symbol = NULL, since no values are considered censored.} +\item{impute}{only for method = "TMP"; imputes censored values via AFT model when TRUE} -\item{censored_symbol}{Indicates how censored missing values are encoded in -the 'Intensity' column. 'NA' (default) treats all NA intensities as -left-censored (i.e., below the limit of detection). '0' treats zero -intensities as left-censored; in this case NA intensities are assumed to be -missing at random and are not censored. Skyline output should use '0'. NULL -assumes that all missing values are missing at random — no values are treated -as censored, and imputation is disabled (impute is ignored).} +\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins -where every run has at least 50\% missing values for each peptide. FALSE is default.} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} -\item{equal_variance}{only for summaryMethod = "linear". Default is TRUE. -Logical variable for whether the model should account for heterogeneous variation -among intensities from different features. Default is TRUE, which assume equal -variance among intensities from features. FALSE means that we cannot assume -equal variance among intensities from features, then we will account for -heterogeneous variation from different features.} +\item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities} \item{aft_iterations}{Number of iterations for AFT model fitting. Default is 90.} } diff --git a/man/dot-build_summarize_worker.Rd b/man/dot-build_summarize_worker.Rd index 8970d6a9..f8a526bf 100644 --- a/man/dot-build_summarize_worker.Rd +++ b/man/dot-build_summarize_worker.Rd @@ -2,8 +2,7 @@ % Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R \name{.build_summarize_worker} \alias{.build_summarize_worker} -\title{Build the per-record worker closure for -\code{MSstatsSummarizeWithMultipleCores}} +\title{Build the per-record worker closure for \code{MSstatsSummarizeWithMultipleCores}} \usage{ .build_summarize_worker( use_TMP, @@ -15,18 +14,7 @@ ) } \description{ -Defined at package top level — not nested via \code{local()} inside -\code{MSstatsSummarizeWithMultipleCores} — so the returned closure's -enclosing environment chain is this factory's own (small) call frame plus -the package namespace. \code{BiocParallel}/\code{matter} serialize a -closure's entire enclosing environment chain to ship it to each socket -worker, not just the variables the closure body actually references. A -closure built via \code{local()} inside -\code{MSstatsSummarizeWithMultipleCores} would have that function's own -evaluation frame in its chain — which holds \code{input}, -\code{protein_records}, and other run-scale objects — so every worker -would receive and retain a serialized copy of them even though -\code{.worker} never touches them. Building the closure here instead keeps -its captured state to only the scalar run parameters. +Defined at package top level so the closure only captures the scalar run +parameters, not the caller's run-scale objects. } \keyword{internal} diff --git a/man/dot-getNonMissingFilterStats.Rd b/man/dot-getNonMissingFilterStats.Rd index 7f8b8970..b4ca502f 100644 --- a/man/dot-getNonMissingFilterStats.Rd +++ b/man/dot-getNonMissingFilterStats.Rd @@ -9,13 +9,7 @@ \arguments{ \item{input}{data.table with data for a single protein} -\item{censored_symbol}{Indicates how censored missing values are encoded in -the 'Intensity' column. 'NA' (default) treats all NA intensities as -left-censored (i.e., below the limit of detection). '0' treats zero -intensities as left-censored; in this case NA intensities are assumed to be -missing at random and are not censored. Skyline output should use '0'. NULL -assumes that all missing values are missing at random — no values are treated -as censored, and imputation is disabled (impute is ignored).} +\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} } \value{ data.table diff --git a/man/dot-pack_protein_slot.Rd b/man/dot-pack_protein_slot.Rd index 70849c0e..a9b74d13 100644 --- a/man/dot-pack_protein_slot.Rd +++ b/man/dot-pack_protein_slot.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R \name{.pack_protein_slot} \alias{.pack_protein_slot} -\title{Build the packed double vector and lightweight metadata for one protein slot} +\title{Pack one protein slot into a double vector plus metadata} \usage{ .pack_protein_slot(protein_dt, slot_index, all_runs) } \arguments{ -\item{protein_dt}{data.table rows belonging to one protein (or protein × label) slot} +\item{protein_dt}{data.table rows for one protein (or protein x label) slot} \item{slot_index}{integer position of this slot in the global protein list} @@ -17,6 +17,6 @@ list with elements \code{packed} (double vector) and \code{meta} (list) } \description{ -Build the packed double vector and lightweight metadata for one protein slot +Pack one protein slot into a double vector plus metadata } \keyword{internal} diff --git a/man/dot-runTukey.Rd b/man/dot-runTukey.Rd index b6425a04..0e712023 100644 --- a/man/dot-runTukey.Rd +++ b/man/dot-runTukey.Rd @@ -15,16 +15,9 @@ subtracting the H value and adding back the H median, and only L results are returned. If FALSE (e.g. protein turnover), each label is summarized independently and results for all labels are returned.} -\item{censored_symbol}{Indicates how censored missing values are encoded in -the 'Intensity' column. 'NA' (default) treats all NA intensities as -left-censored (i.e., below the limit of detection). '0' treats zero -intensities as left-censored; in this case NA intensities are assumed to be -missing at random and are not censored. Skyline output should use '0'. NULL -assumes that all missing values are missing at random — no values are treated -as censored, and imputation is disabled (impute is ignored).} +\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins -where every run has at least 50\% missing values for each peptide. FALSE is default.} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} } \value{ data.table From 70113877239039a3975e6b80a94709955dc6d1d2 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 15:29:00 -0400 Subject: [PATCH 07/18] remove cen := as.integer --- R/MSstatsSummarizeWithMultipleCores.R | 2 -- man/reexports.Rd | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 9c694777..4ad59cec 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -290,8 +290,6 @@ result[[idx]][, FEATURE := droplevels(factor(as.character(FEATURE), levels = feature_levels))] } } - if (!is.null(result[[2L]]) && "cen" %in% colnames(result[[2L]])) - result[[2L]][, cen := as.integer(cen)] result } } diff --git a/man/reexports.Rd b/man/reexports.Rd index 04f47fc4..eeac4273 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -21,6 +21,6 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} + \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} }} From 3e91b417c5a615ecc5e36792bc4f5c39a5c16383 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 15:41:29 -0400 Subject: [PATCH 08/18] clean up factoring --- R/MSstatsSummarizeWithMultipleCores.R | 11 +---------- man/reexports.Rd | 2 +- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 4ad59cec..aa19cc64 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -281,15 +281,6 @@ remove50missing_, aft_iterations_, equal_variances = equal_variance_) } - feature_levels <- sort(unique(meta$feature_label_dt$FEATURE)) - for (idx in 1:2) { - if (!is.null(result[[idx]])) { - if ("RUN" %in% colnames(result[[idx]])) - result[[idx]][, RUN := droplevels(factor(as.character(RUN), levels = meta$runs))] - if ("FEATURE" %in% colnames(result[[idx]])) - result[[idx]][, FEATURE := droplevels(factor(as.character(FEATURE), levels = feature_levels))] - } - } result } } @@ -331,7 +322,7 @@ #' @importFrom matter SnowfastParam #' @importFrom BiocParallel bplapply bpstart bpstop bpisup bpnworkers bpprogressbar #' @importFrom data.table data.table fifelse setDTthreads -#' @importFrom stats median predict +#' @importFrom stats median #' #' @export MSstatsSummarizeWithMultipleCores <- function( diff --git a/man/reexports.Rd b/man/reexports.Rd index eeac4273..04f47fc4 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -21,6 +21,6 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} + \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} }} From 5da0bf9b17e4791f31fc9b78ee2b559946f96909 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 15:45:48 -0400 Subject: [PATCH 09/18] fix docs --- R/MSstatsSummarizeWithMultipleCores.R | 2 +- R/utils_censored.R | 2 +- R/utils_summarization.R | 2 +- man/MSstatsSummarizeSingleTMP.Rd | 2 +- man/MSstatsSummarizeWithMultipleCores.Rd | 2 +- man/MSstatsSummarizeWithSingleCore.Rd | 2 +- man/dot-isSummarizable.Rd | 2 +- man/dot-runTukey.Rd | 2 +- man/dot-setCensoredByThreshold.Rd | 2 +- man/reexports.Rd | 2 +- src/RcppExports-1156a8df.o.tmp | 0 11 files changed, 10 insertions(+), 10 deletions(-) create mode 100644 src/RcppExports-1156a8df.o.tmp diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index aa19cc64..c33a4cba 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -306,7 +306,7 @@ #' @param method summarization method: "linear" or "TMP" #' @param impute only for method = "TMP"; imputes censored values via AFT model when TRUE #' @param censored_symbol how censored values are encoded: 'NA', '0', or NULL for none -#' @param remove50missing only for method = "TMP"; drops proteins missing >=50% per peptide in every run +#' @param remove50missing only for method = "TMP"; drops proteins missing >=50\% per peptide in every run #' @param equal_variance only for method = "linear"; assume equal variance among feature intensities #' @param numberOfCores number of cores for parallel processing (Linux/Mac only) #' @param aft_iterations number of AFT model iterations diff --git a/R/utils_censored.R b/R/utils_censored.R index 73b81cb2..252708fe 100644 --- a/R/utils_censored.R +++ b/R/utils_censored.R @@ -84,7 +84,7 @@ MSstatsHandleMissing = function(input, summary_method, impute, #' possible values for left-censored data as the `time` input to the Surv function. #' @param input `data.table` in MSstats format #' @param censored_symbol censoredInt parameter to `dataProcess` -#' @param remove50missing if TRUE, features with at least 50% missing values +#' @param remove50missing if TRUE, features with at least 50\% missing values #' will be removed #' @keywords internal .setCensoredByThreshold = function(input, censored_symbol, remove50missing) { diff --git a/R/utils_summarization.R b/R/utils_summarization.R index 47bbef21..26225316 100644 --- a/R/utils_summarization.R +++ b/R/utils_summarization.R @@ -1,6 +1,6 @@ #' Check if a protein can be summarized with TMP #' @param input data.table -#' @param remove50missing if TRUE, proteins with more than 50% missing values +#' @param remove50missing if TRUE, proteins with more than 50\% missing values #' in all runs will not be summarized #' @return data.table #' @keywords internal diff --git a/man/MSstatsSummarizeSingleTMP.Rd b/man/MSstatsSummarizeSingleTMP.Rd index 26628f2e..cd115723 100644 --- a/man/MSstatsSummarizeSingleTMP.Rd +++ b/man/MSstatsSummarizeSingleTMP.Rd @@ -19,7 +19,7 @@ MSstatsSummarizeSingleTMP( \item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run} \item{aft_iterations}{number of iterations for AFT model fitting} } diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index f77134c1..d4130f11 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -28,7 +28,7 @@ MSstatsSummarizeWithMultipleCores( \item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run} \item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities} diff --git a/man/MSstatsSummarizeWithSingleCore.Rd b/man/MSstatsSummarizeWithSingleCore.Rd index 2c3f23ef..fc711d42 100644 --- a/man/MSstatsSummarizeWithSingleCore.Rd +++ b/man/MSstatsSummarizeWithSingleCore.Rd @@ -23,7 +23,7 @@ MSstatsSummarizeWithSingleCore( \item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run} \item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities} diff --git a/man/dot-isSummarizable.Rd b/man/dot-isSummarizable.Rd index b30b49cf..b13e68fb 100644 --- a/man/dot-isSummarizable.Rd +++ b/man/dot-isSummarizable.Rd @@ -9,7 +9,7 @@ \arguments{ \item{input}{data.table} -\item{remove50missing}{if TRUE, proteins with more than 50% missing values +\item{remove50missing}{if TRUE, proteins with more than 50\% missing values in all runs will not be summarized} } \value{ diff --git a/man/dot-runTukey.Rd b/man/dot-runTukey.Rd index 0e712023..c1049973 100644 --- a/man/dot-runTukey.Rd +++ b/man/dot-runTukey.Rd @@ -17,7 +17,7 @@ independently and results for all labels are returned.} \item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none} -\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50% per peptide in every run} +\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run} } \value{ data.table diff --git a/man/dot-setCensoredByThreshold.Rd b/man/dot-setCensoredByThreshold.Rd index b3ca6e63..45e66237 100644 --- a/man/dot-setCensoredByThreshold.Rd +++ b/man/dot-setCensoredByThreshold.Rd @@ -13,7 +13,7 @@ possible values for left-censored data as the `time` input to the Surv function. \item{censored_symbol}{censoredInt parameter to `dataProcess`} -\item{remove50missing}{if TRUE, features with at least 50% missing values +\item{remove50missing}{if TRUE, features with at least 50\% missing values will be removed} } \description{ diff --git a/man/reexports.Rd b/man/reexports.Rd index 04f47fc4..eeac4273 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -21,6 +21,6 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} + \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} }} diff --git a/src/RcppExports-1156a8df.o.tmp b/src/RcppExports-1156a8df.o.tmp new file mode 100644 index 00000000..e69de29b From 1bf540783068d71fde4e837bff22d9fa7cf1921e Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 15:46:24 -0400 Subject: [PATCH 10/18] remove tmp file from src --- src/RcppExports-1156a8df.o.tmp | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/RcppExports-1156a8df.o.tmp diff --git a/src/RcppExports-1156a8df.o.tmp b/src/RcppExports-1156a8df.o.tmp deleted file mode 100644 index e69de29b..00000000 From a23ac83ee92fc5ab9b1bdba3c99878aacfa85179 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 15:56:07 -0400 Subject: [PATCH 11/18] migrate c++ to cpp files --- R/MSstatsSummarizeWithMultipleCores.R | 40 ++------------------------- R/RcppExports.R | 4 +++ src/Makevars.win | 2 +- src/RcppExports.cpp | 11 ++++++++ src/peak_rss.cpp | 30 ++++++++++++++++++++ 5 files changed, 48 insertions(+), 39 deletions(-) create mode 100644 src/peak_rss.cpp diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index c33a4cba..fe47e642 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -1,45 +1,9 @@ .peak_rss_mb <- function() { - if (.Platform$OS.type == "windows") { - if (!exists(".peakRSS_windows_impl", mode = "function")) { - Rcpp::cppFunction( - depends = "Rcpp", - includes = c( - "#include ", - "#include "), - code = " - double peakRSS_windows_impl() { - PROCESS_MEMORY_COUNTERS pmc; - if (GetProcessMemoryInfo(GetCurrentProcess(), &pmc, sizeof(pmc))) { - return (double) pmc.PeakWorkingSetSize / (1024.0*1024.0); - } - return NA_REAL; - }") - assign(".peakRSS_windows_impl", peakRSS_windows_impl, envir = .GlobalEnv) - } - return(.peakRSS_windows_impl()) - } - - if (file.exists("/proc/self/status")) { + if (.Platform$OS.type != "windows" && file.exists("/proc/self/status")) { ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024) } - - if (!exists(".rusage_maxrss_mb_impl", mode = "function")) { - Rcpp::cppFunction( - depends = "Rcpp", - includes = "#include ", - code = " - double rusage_maxrss_mb_impl() { - struct rusage ru; getrusage(RUSAGE_SELF, &ru); - #ifdef __APPLE__ - return (double) ru.ru_maxrss / (1024.0*1024.0); - #else - return (double) ru.ru_maxrss / 1024.0; - #endif - }") - assign(".rusage_maxrss_mb_impl", rusage_maxrss_mb_impl, envir = .GlobalEnv) - } - .rusage_maxrss_mb_impl() + peak_rss_mb() } .print_memory_report <- function(function_name, checkpoints, worker_peak_mb = NULL, diff --git a/R/RcppExports.R b/R/RcppExports.R index 65f2c4fb..26dbc834 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -17,3 +17,7 @@ median_polish_summary <- function(x, eps = 0.01, maxiter = 10L) { .Call(`_MSstats_median_polish_summary`, x, eps, maxiter) } +peak_rss_mb <- function() { + .Call(`_MSstats_peak_rss_mb`) +} + diff --git a/src/Makevars.win b/src/Makevars.win index 6660c7f9..eac6e524 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,3 +1,3 @@ CXX_STD = CXX14 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lpsapi diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5275017b..a27a4301 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -67,12 +67,23 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// peak_rss_mb +double peak_rss_mb(); +RcppExport SEXP _MSstats_peak_rss_mb() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(peak_rss_mb()); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_MSstats_get_estimable_fixed_random", (DL_FUNC) &_MSstats_get_estimable_fixed_random, 2}, {"_MSstats_make_contrast_run_quant", (DL_FUNC) &_MSstats_make_contrast_run_quant, 6}, {"_MSstats_get_linear_summary", (DL_FUNC) &_MSstats_get_linear_summary, 5}, {"_MSstats_median_polish_summary", (DL_FUNC) &_MSstats_median_polish_summary, 3}, + {"_MSstats_peak_rss_mb", (DL_FUNC) &_MSstats_peak_rss_mb, 0}, {NULL, NULL, 0} }; diff --git a/src/peak_rss.cpp b/src/peak_rss.cpp new file mode 100644 index 00000000..1dcd9277 --- /dev/null +++ b/src/peak_rss.cpp @@ -0,0 +1,30 @@ +#include +using namespace Rcpp; + +#if defined(_WIN32) +#include +#include +#elif defined(__unix__) || defined(__APPLE__) +#include +#endif + +// [[Rcpp::export]] +double peak_rss_mb() { +#if defined(_WIN32) + PROCESS_MEMORY_COUNTERS pmc; + if (GetProcessMemoryInfo(GetCurrentProcess(), &pmc, sizeof(pmc))) { + return (double) pmc.PeakWorkingSetSize / (1024.0 * 1024.0); + } + return NA_REAL; +#elif defined(__unix__) || defined(__APPLE__) + struct rusage ru; + getrusage(RUSAGE_SELF, &ru); +#ifdef __APPLE__ + return (double) ru.ru_maxrss / (1024.0 * 1024.0); +#else + return (double) ru.ru_maxrss / 1024.0; +#endif +#else + return NA_REAL; +#endif +} From cfb812ea9c72c088ac4b76630a51f4eb86018d21 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 16:32:50 -0400 Subject: [PATCH 12/18] add hpc script for benchmarking after merging --- benchmark/benchmark_summarize_perf_selevsek.R | 34 +++++++++++++++++++ benchmark/config.slurm | 2 +- 2 files changed, 35 insertions(+), 1 deletion(-) create mode 100644 benchmark/benchmark_summarize_perf_selevsek.R diff --git a/benchmark/benchmark_summarize_perf_selevsek.R b/benchmark/benchmark_summarize_perf_selevsek.R new file mode 100644 index 00000000..cb628b33 --- /dev/null +++ b/benchmark/benchmark_summarize_perf_selevsek.R @@ -0,0 +1,34 @@ +library(MSstats) + +# Performance expectations for MSstatsSummarizeWithMultipleCores on this dataset +# (not enforced automatically here -- read the job output and compare by hand): +# - 4 cores should reduce wall time by at least 25% vs. 1 core +# - peak RSS per worker should stay under ~1GB (see the memory report below) + +input <- data.table::fread("/projects/VitekLab/Data/MS/selevsek/before_summarization.csv") + +cat("=== MSstatsSummarizeWithMultipleCores performance check ===\n") +cat("Expectation: >=25% reduction in wall time on 4 cores vs. 1 core.\n") +cat("Expectation: peak RSS per worker should stay under ~1GB (see memory report below).\n\n") + +time_1core <- system.time( + result_1core <- MSstatsSummarizeWithSingleCore( + input, "TMP", TRUE, "NA", FALSE, TRUE + ) +) + +time_4core <- system.time( + result_4core <- MSstatsSummarizeWithMultipleCores( + input, "TMP", TRUE, "NA", FALSE, TRUE, 4, track_memory = TRUE, + max_proteins_per_worker = 200 + ) +) + +elapsed_1core <- time_1core[["elapsed"]] +elapsed_4core <- time_4core[["elapsed"]] +speedup_pct <- (elapsed_1core - elapsed_4core) / elapsed_1core * 100 + +cat(sprintf("1 core wall time: %.2f s\n", elapsed_1core)) +cat(sprintf("4 core wall time: %.2f s\n", elapsed_4core)) +cat(sprintf("Observed speedup: %.1f%% (expectation: >= 25%%)\n", speedup_pct)) +cat("Compare the 'Worker RSS min / mean / max' line above against the ~1GB expectation.\n") diff --git a/benchmark/config.slurm b/benchmark/config.slurm index 28332e1a..28355762 100644 --- a/benchmark/config.slurm +++ b/benchmark/config.slurm @@ -27,7 +27,7 @@ remotes::install_github('Vitek-Lab/MSstats', ref = 'devel', lib = Sys.getenv('R_ remotes::install_github('Vitek-Lab/MSstatsConvert', ref = 'master', lib = Sys.getenv('R_LIBS_USER')); \ install.packages(c('dplyr', 'stringr', 'ggplot2'), lib = Sys.getenv('R_LIBS_USER'), repos = 'https://cloud.r-project.org')" -R_SCRIPTS=("benchmark_Dowell2021-HEqe408_LFQ.R" "benchmark_Puyvelde2022-HYE5600735_LFQ.R") +R_SCRIPTS=("benchmark_Dowell2021-HEqe408_LFQ.R" "benchmark_Puyvelde2022-HYE5600735_LFQ.R" "benchmark_summarize_perf_selevsek.R") for script in "${R_SCRIPTS[@]}"; do echo "Executing script: $script" >> job_output.txt From e33403adcbcf5b6c5b5ceddf4282366786090d94 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 16:49:04 -0400 Subject: [PATCH 13/18] add comparison infrastructure for benchmarking script for selevsek --- benchmark/benchmark_summarize_perf_selevsek.R | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/benchmark/benchmark_summarize_perf_selevsek.R b/benchmark/benchmark_summarize_perf_selevsek.R index cb628b33..c36cc2f8 100644 --- a/benchmark/benchmark_summarize_perf_selevsek.R +++ b/benchmark/benchmark_summarize_perf_selevsek.R @@ -32,3 +32,57 @@ cat(sprintf("1 core wall time: %.2f s\n", elapsed_1core)) cat(sprintf("4 core wall time: %.2f s\n", elapsed_4core)) cat(sprintf("Observed speedup: %.1f%% (expectation: >= 25%%)\n", speedup_pct)) cat("Compare the 'Worker RSS min / mean / max' line above against the ~1GB expectation.\n") + +cat("\n=== Verifying 1-core and 4-core results are identical ===\n") + +compare_df <- function(df1, df2) { + df1 <- as.data.frame(df1) + df2 <- as.data.frame(df2) + + if (!setequal(names(df1), names(df2))) return(FALSE) + + df2 <- df2[, names(df1), drop = FALSE] + + # Convert factors to character so ordering/comparison is based on + # content, not on (possibly differing) factor level order + df1[] <- lapply(df1, function(x) if (is.factor(x)) as.character(x) else x) + df2[] <- lapply(df2, function(x) if (is.factor(x)) as.character(x) else x) + + ord1 <- do.call(order, as.list(df1)) + ord2 <- do.call(order, as.list(df2)) + + df1 <- df1[ord1, , drop = FALSE] + df2 <- df2[ord2, , drop = FALSE] + + rownames(df1) <- NULL + rownames(df2) <- NULL + + isTRUE(all.equal(df1, df2, tolerance = 1e-14)) +} + +n_proteins <- length(result_1core) + +protein_level_matches <- logical(n_proteins) +for (i in seq_len(n_proteins)) { + protein_level_matches[i] <- compare_df(result_1core[[i]][[1]], result_4core[[i]][[1]]) +} + +survival_data_matches <- logical(n_proteins) +for (i in seq_len(n_proteins)) { + survival_data_matches[i] <- compare_df(result_1core[[i]][[2]], result_4core[[i]][[2]]) +} + +n_protein_level_match <- sum(protein_level_matches) +n_survival_data_match <- sum(survival_data_matches) + +cat(sprintf("Protein-level results matching: %d / %d\n", n_protein_level_match, n_proteins)) +cat(sprintf("Survival data matching: %d / %d\n", n_survival_data_match, n_proteins)) + +stopifnot( + "Protein-level results differ between 1 core and 4 cores" = + n_protein_level_match == n_proteins, + "Survival data differs between 1 core and 4 cores" = + n_survival_data_match == n_proteins +) + +cat("1-core and 4-core results are identical.\n") From 42e0b07d8553447c720b96e27d93f519a1c44446 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Wed, 5 Aug 2026 16:51:05 -0400 Subject: [PATCH 14/18] fix benchmarking script to say feature level instead of survival --- benchmark/benchmark_summarize_perf_selevsek.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/benchmark/benchmark_summarize_perf_selevsek.R b/benchmark/benchmark_summarize_perf_selevsek.R index c36cc2f8..16180ac1 100644 --- a/benchmark/benchmark_summarize_perf_selevsek.R +++ b/benchmark/benchmark_summarize_perf_selevsek.R @@ -67,22 +67,22 @@ for (i in seq_len(n_proteins)) { protein_level_matches[i] <- compare_df(result_1core[[i]][[1]], result_4core[[i]][[1]]) } -survival_data_matches <- logical(n_proteins) +feature_level_data_matches <- logical(n_proteins) for (i in seq_len(n_proteins)) { - survival_data_matches[i] <- compare_df(result_1core[[i]][[2]], result_4core[[i]][[2]]) + feature_level_data_matches[i] <- compare_df(result_1core[[i]][[2]], result_4core[[i]][[2]]) } n_protein_level_match <- sum(protein_level_matches) -n_survival_data_match <- sum(survival_data_matches) +n_feature_level_data_match <- sum(feature_level_data_matches) cat(sprintf("Protein-level results matching: %d / %d\n", n_protein_level_match, n_proteins)) -cat(sprintf("Survival data matching: %d / %d\n", n_survival_data_match, n_proteins)) +cat(sprintf("Feature_level data matching: %d / %d\n", n_feature_level_data_match, n_proteins)) stopifnot( "Protein-level results differ between 1 core and 4 cores" = n_protein_level_match == n_proteins, - "Survival data differs between 1 core and 4 cores" = - n_survival_data_match == n_proteins + "Feature_level data differs between 1 core and 4 cores" = + n_feature_level_data_match == n_proteins ) cat("1-core and 4-core results are identical.\n") From 549279e2a5ea4bc1f68e153dfb9db810d8085a12 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Thu, 6 Aug 2026 12:11:24 -0400 Subject: [PATCH 15/18] change default for MSstatsSummarizeMultipleCores to 50 --- R/MSstatsSummarizeWithMultipleCores.R | 4 ++-- man/MSstatsSummarizeWithMultipleCores.Rd | 4 ++-- man/reexports.Rd | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index fe47e642..4674449a 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -278,7 +278,7 @@ #' @param BPPARAM optional \code{BiocParallelParam} instance #' @param track_memory whether to report per-worker peak RSS memory usage #' @param max_proteins_per_worker caps protein records per \code{bplapply} task; -#' 0 uses BiocParallel's default split +#' 0 uses BiocParallel's default split, default is 50. #' #' @return A named list with one element per protein slot, keyed by protein #' (or protein \eqn{\times} label) identifier. @@ -301,7 +301,7 @@ MSstatsSummarizeWithMultipleCores <- function( verbose = FALSE, BPPARAM = NULL, track_memory = FALSE, - max_proteins_per_worker = 0L + max_proteins_per_worker = 50L ) { if (numberOfCores <= 1L && is.null(BPPARAM)) { return(MSstatsSummarizeWithSingleCore( diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index d4130f11..9f0360fb 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -16,7 +16,7 @@ MSstatsSummarizeWithMultipleCores( verbose = FALSE, BPPARAM = NULL, track_memory = FALSE, - max_proteins_per_worker = 0L + max_proteins_per_worker = 50L ) } \arguments{ @@ -43,7 +43,7 @@ MSstatsSummarizeWithMultipleCores( \item{track_memory}{whether to report per-worker peak RSS memory usage} \item{max_proteins_per_worker}{caps protein records per \code{bplapply} task; -0 uses BiocParallel's default split} +0 uses BiocParallel's default split, default is 50.} } \value{ A named list with one element per protein slot, keyed by protein diff --git a/man/reexports.Rd b/man/reexports.Rd index eeac4273..04f47fc4 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -21,6 +21,6 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} + \item{MSstatsConvert}{\code{\link[MSstatsConvert:DIANNtoMSstatsFormat]{DIANNtoMSstatsFormat()}}, \code{\link[MSstatsConvert:DIAUmpiretoMSstatsFormat]{DIAUmpiretoMSstatsFormat()}}, \code{\link[MSstatsConvert:FragPipetoMSstatsFormat]{FragPipetoMSstatsFormat()}}, \code{\link[MSstatsConvert:MaxQtoMSstatsFormat]{MaxQtoMSstatsFormat()}}, \code{\link[MSstatsConvert:MZMinetoMSstatsFormat]{MZMinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenMStoMSstatsFormat]{OpenMStoMSstatsFormat()}}, \code{\link[MSstatsConvert:OpenSWATHtoMSstatsFormat]{OpenSWATHtoMSstatsFormat()}}, \code{\link[MSstatsConvert:PDtoMSstatsFormat]{PDtoMSstatsFormat()}}, \code{\link[MSstatsConvert:ProgenesistoMSstatsFormat]{ProgenesistoMSstatsFormat()}}, \code{\link[MSstatsConvert:SkylinetoMSstatsFormat]{SkylinetoMSstatsFormat()}}, \code{\link[MSstatsConvert:SpectronauttoMSstatsFormat]{SpectronauttoMSstatsFormat()}}} }} From 882e360c9f52098d7d003988468fbc99fe7cdfc0 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 17 Aug 2026 16:30:43 -0400 Subject: [PATCH 16/18] reset rss mb for linux machines for memory benchmarking --- R/MSstatsSummarizeWithMultipleCores.R | 50 ++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 4674449a..d8e9fe73 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -6,6 +6,35 @@ peak_rss_mb() } +#' Reset the process peak-RSS high-water mark to the current RSS +#' +#' The kernel-tracked peak RSS (\code{VmHWM} on Linux, \code{ru_maxrss} on +#' macOS/BSD, \code{PeakWorkingSetSize} on Windows) is a running maximum +#' since process start. Without a reset, an earlier one-off allocation in the +#' same R session (e.g. reading an 8 GB CSV) inflates every peak-memory +#' measurement taken afterwards even if that memory has since been freed, +#' making it useless as a baseline for benchmarking this function's own +#' memory use. +#' +#' Linux only: writing \code{"5"} to \code{/proc/self/clear_refs} resets +#' \code{VmHWM} to the current RSS (kernel >= 4.0). macOS and Windows have no +#' equivalent public facility to reset their peak-RSS counters, so this is a +#' no-op there, and \code{track_memory} reports on those platforms may still +#' include memory used before this function was entered. +#' +#' @return invisible \code{TRUE} if the peak was reset, \code{FALSE} otherwise +#' @keywords internal +.reset_peak_rss <- function() { + if (Sys.info()[["sysname"]] == "Linux" && file.exists("/proc/self/clear_refs")) { + reset_ok <- tryCatch({ + cat("5", file = "/proc/self/clear_refs") + TRUE + }, error = function(e) FALSE, warning = function(w) FALSE) + return(invisible(reset_ok)) + } + invisible(FALSE) +} + .print_memory_report <- function(function_name, checkpoints, worker_peak_mb = NULL, elapsed = NULL) { rule_width <- 65L @@ -276,7 +305,14 @@ #' @param aft_iterations number of AFT model iterations #' @param verbose whether to print verbose output #' @param BPPARAM optional \code{BiocParallelParam} instance -#' @param track_memory whether to report per-worker peak RSS memory usage +#' @param track_memory whether to report per-worker peak RSS memory usage. +#' On Linux, the process peak RSS is reset on entry (via +#' \code{/proc/self/clear_refs}) so the reported baseline reflects memory +#' at function entry rather than a historical high-water mark inflated by +#' earlier operations in the session (e.g. reading a large input file). +#' This reset is Linux-only; on macOS and Windows there is no supported +#' facility to reset the peak-RSS counter, so reports there may still +#' include memory used before this function was entered. #' @param max_proteins_per_worker caps protein records per \code{bplapply} task; #' 0 uses BiocParallel's default split, default is 50. #' @@ -312,6 +348,18 @@ MSstatsSummarizeWithMultipleCores <- function( start_time <- proc.time()[["elapsed"]] memory_checkpoints <- list() + if (track_memory) { + peak_reset <- .reset_peak_rss() + if (!peak_reset) { + getOption("MSstatsLog")("INFO", + paste0("Peak RSS reset on entry is only supported on Linux ", + "(via /proc/self/clear_refs); on this platform, the ", + "peak-memory report below may include usage from ", + "before this function was entered.")) + } + memory_checkpoints[["baseline (peak reset on entry)"]] <- .peak_rss_mb() + } + is_labeled_reference <- "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE) split_keys <- if (is_labeled_reference) list(input$PROTEIN) else From e6f2eb05de672faa58c0a36d71616e27d1ffbe69 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Tue, 1 Sep 2026 09:54:21 -0400 Subject: [PATCH 17/18] add documentation, rename peak to max --- R/MSstatsSummarizeWithMultipleCores.R | 80 ++++++++++++------- benchmark/benchmark_summarize_perf_selevsek.R | 5 +- man/MSstatsSummarizeWithMultipleCores.Rd | 9 ++- man/dot-max_rss_mb.Rd | 30 +++++++ man/dot-reset_max_rss.Rd | 29 +++++++ src/peak_rss.cpp | 11 +++ 6 files changed, 133 insertions(+), 31 deletions(-) create mode 100644 man/dot-max_rss_mb.Rd create mode 100644 man/dot-reset_max_rss.Rd diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index d8e9fe73..7c0736fb 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -1,4 +1,24 @@ -.peak_rss_mb <- function() { +#' Maximum RAM (in MB) this R process has ever used +#' +#' There are two paths: +#' +#' \enumerate{ +#' \item \strong{Linux} (non-Windows and \code{/proc/self/status} exists): +#' read the number Linux itself keeps. \code{/proc/self/status} is a +#' plain-text file of \code{Key: value} lines describing this process; the +#' line \code{VmHWM:} holds the maximum RAM ever used, in kB. Preferred +#' because \code{/proc/self/clear_refs} (see \code{\link{.reset_max_rss}}) +#' can reset it, so the reading reflects only recent work. The function +#' pulls the number out of a line that looks like +#' \code{"VmHWM: 123456 kB"}. +#' \item \strong{Everything else} (macOS, Windows, or Linux without +#' \code{/proc}): fall back to the compiled \code{peak_rss_mb()} in +#' \code{src/peak_rss.cpp}, which asks the operating system directly. +#' } +#' +#' @return maximum RAM used, in MB +#' @keywords internal +.max_rss_mb <- function() { if (.Platform$OS.type != "windows" && file.exists("/proc/self/status")) { ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024) @@ -6,25 +26,26 @@ peak_rss_mb() } -#' Reset the process peak-RSS high-water mark to the current RSS +#' Reset the process maximum-RSS high-water mark to the current RSS #' -#' The kernel-tracked peak RSS (\code{VmHWM} on Linux, \code{ru_maxrss} on +#' The kernel-tracked maximum RSS (\code{VmHWM} on Linux, \code{ru_maxrss} on #' macOS/BSD, \code{PeakWorkingSetSize} on Windows) is a running maximum #' since process start. Without a reset, an earlier one-off allocation in the -#' same R session (e.g. reading an 8 GB CSV) inflates every peak-memory +#' same R session (e.g. reading an 8 GB CSV) inflates every maximum-memory #' measurement taken afterwards even if that memory has since been freed, #' making it useless as a baseline for benchmarking this function's own #' memory use. #' #' Linux only: writing \code{"5"} to \code{/proc/self/clear_refs} resets #' \code{VmHWM} to the current RSS (kernel >= 4.0). macOS and Windows have no -#' equivalent public facility to reset their peak-RSS counters, so this is a +#' equivalent public facility to reset their maximum-RSS counters, so this is a #' no-op there, and \code{track_memory} reports on those platforms may still #' include memory used before this function was entered. #' -#' @return invisible \code{TRUE} if the peak was reset, \code{FALSE} otherwise +#' @return invisible \code{TRUE} if the high-water mark was reset, +#' \code{FALSE} otherwise #' @keywords internal -.reset_peak_rss <- function() { +.reset_max_rss <- function() { if (Sys.info()[["sysname"]] == "Linux" && file.exists("/proc/self/clear_refs")) { reset_ok <- tryCatch({ cat("5", file = "/proc/self/clear_refs") @@ -35,7 +56,7 @@ invisible(FALSE) } -.print_memory_report <- function(function_name, checkpoints, worker_peak_mb = NULL, +.print_memory_report <- function(function_name, checkpoints, worker_max_rss_mb = NULL, elapsed = NULL) { rule_width <- 65L rule <- strrep("─", rule_width) @@ -59,12 +80,13 @@ format_delta(checkpoint_value, previous_value))) previous_value <- checkpoint_value } - if (!is.null(worker_peak_mb)) { - observed_peaks <- worker_peak_mb[!is.na(worker_peak_mb)] - if (length(observed_peaks)) { + if (!is.null(worker_max_rss_mb)) { + observed_max_memory <- worker_max_rss_mb[!is.na(worker_max_rss_mb)] + if (length(observed_max_memory)) { lines <- c(lines, "", sprintf(" Worker RSS min / mean / max : %.1f / %.1f / %.1f MB", - min(observed_peaks), mean(observed_peaks), max(observed_peaks))) + min(observed_max_memory), mean(observed_max_memory), + max(observed_max_memory))) } } if (!is.null(elapsed)) @@ -278,12 +300,12 @@ } } -#' Per-worker peak-RSS query task for \code{MSstatsSummarizeWithMultipleCores} +#' Per-worker maximum-RSS query task for \code{MSstatsSummarizeWithMultipleCores} #' #' @keywords internal #' @noRd -.report_worker_peak <- function(i) { - list(worker = i, pid = Sys.getpid(), peak_mb = .peak_rss_mb()) +.report_worker_max_rss <- function(i) { + list(worker = i, pid = Sys.getpid(), max_rss_mb = .max_rss_mb()) } .warmup_worker <- function(i) { @@ -305,13 +327,13 @@ #' @param aft_iterations number of AFT model iterations #' @param verbose whether to print verbose output #' @param BPPARAM optional \code{BiocParallelParam} instance -#' @param track_memory whether to report per-worker peak RSS memory usage. -#' On Linux, the process peak RSS is reset on entry (via +#' @param track_memory whether to report per-worker maximum RSS memory usage. +#' On Linux, the process maximum RSS is reset on entry (via #' \code{/proc/self/clear_refs}) so the reported baseline reflects memory #' at function entry rather than a historical high-water mark inflated by #' earlier operations in the session (e.g. reading a large input file). #' This reset is Linux-only; on macOS and Windows there is no supported -#' facility to reset the peak-RSS counter, so reports there may still +#' facility to reset the maximum-RSS counter, so reports there may still #' include memory used before this function was entered. #' @param max_proteins_per_worker caps protein records per \code{bplapply} task; #' 0 uses BiocParallel's default split, default is 50. @@ -349,15 +371,15 @@ MSstatsSummarizeWithMultipleCores <- function( memory_checkpoints <- list() if (track_memory) { - peak_reset <- .reset_peak_rss() - if (!peak_reset) { + max_rss_reset <- .reset_max_rss() + if (!max_rss_reset) { getOption("MSstatsLog")("INFO", - paste0("Peak RSS reset on entry is only supported on Linux ", + paste0("Maximum RSS reset on entry is only supported on Linux ", "(via /proc/self/clear_refs); on this platform, the ", - "peak-memory report below may include usage from ", + "maximum-memory report below may include usage from ", "before this function was entered.")) } - memory_checkpoints[["baseline (peak reset on entry)"]] <- .peak_rss_mb() + memory_checkpoints[["baseline (max RSS reset on entry)"]] <- .max_rss_mb() } is_labeled_reference <- "is_labeled_ref" %in% colnames(input) && @@ -429,18 +451,18 @@ MSstatsSummarizeWithMultipleCores <- function( results <- BiocParallel::bplapply(protein_records, worker_fn, BPPARAM = BPPARAM) names(results) <- protein_ids - worker_peaks <- NULL + worker_max_rss <- NULL if (track_memory) { BiocParallel::bpprogressbar(BPPARAM) <- FALSE - worker_peaks <- BiocParallel::bplapply( + worker_max_rss <- BiocParallel::bplapply( seq_len(BiocParallel::bpnworkers(BPPARAM)), - .report_worker_peak, BPPARAM = BPPARAM) + .report_worker_max_rss, BPPARAM = BPPARAM) BiocParallel::bpprogressbar(BPPARAM) <- show_progress - memory_checkpoints[["parent peak (main)"]] <- .peak_rss_mb() - worker_peak_mb <- vapply(worker_peaks, function(x) x$peak_mb, numeric(1L)) + memory_checkpoints[["parent maximum (main)"]] <- .max_rss_mb() + worker_max_rss_mb <- vapply(worker_max_rss, function(x) x$max_rss_mb, numeric(1L)) .print_memory_report( "MSstatsSummarizeWithMultipleCores", - memory_checkpoints, worker_peak_mb, + memory_checkpoints, worker_max_rss_mb, elapsed = proc.time()[["elapsed"]] - start_time) } diff --git a/benchmark/benchmark_summarize_perf_selevsek.R b/benchmark/benchmark_summarize_perf_selevsek.R index 16180ac1..5ea301aa 100644 --- a/benchmark/benchmark_summarize_perf_selevsek.R +++ b/benchmark/benchmark_summarize_perf_selevsek.R @@ -5,7 +5,10 @@ library(MSstats) # - 4 cores should reduce wall time by at least 25% vs. 1 core # - peak RSS per worker should stay under ~1GB (see the memory report below) -input <- data.table::fread("/projects/VitekLab/Data/MS/selevsek/before_summarization.csv") +northeastern_high_performance_cluster_file_path <- + "/projects/VitekLab/Data/MS/selevsek/before_summarization.csv" + +input <- data.table::fread(northeastern_high_performance_cluster_file_path) cat("=== MSstatsSummarizeWithMultipleCores performance check ===\n") cat("Expectation: >=25% reduction in wall time on 4 cores vs. 1 core.\n") diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index 9f0360fb..b95b1e81 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -40,7 +40,14 @@ MSstatsSummarizeWithMultipleCores( \item{BPPARAM}{optional \code{BiocParallelParam} instance} -\item{track_memory}{whether to report per-worker peak RSS memory usage} +\item{track_memory}{whether to report per-worker maximum RSS memory usage. +On Linux, the process maximum RSS is reset on entry (via +\code{/proc/self/clear_refs}) so the reported baseline reflects memory +at function entry rather than a historical high-water mark inflated by +earlier operations in the session (e.g. reading a large input file). +This reset is Linux-only; on macOS and Windows there is no supported +facility to reset the maximum-RSS counter, so reports there may still +include memory used before this function was entered.} \item{max_proteins_per_worker}{caps protein records per \code{bplapply} task; 0 uses BiocParallel's default split, default is 50.} diff --git a/man/dot-max_rss_mb.Rd b/man/dot-max_rss_mb.Rd new file mode 100644 index 00000000..db21e5ef --- /dev/null +++ b/man/dot-max_rss_mb.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.max_rss_mb} +\alias{.max_rss_mb} +\title{Maximum RAM (in MB) this R process has ever used} +\usage{ +.max_rss_mb() +} +\value{ +maximum RAM used, in MB +} +\description{ +There are two paths: +} +\details{ +\enumerate{ + \item \strong{Linux} (non-Windows and \code{/proc/self/status} exists): + read the number Linux itself keeps. \code{/proc/self/status} is a + plain-text file of \code{Key: value} lines describing this process; the + line \code{VmHWM:} holds the maximum RAM ever used, in kB. Preferred + because \code{/proc/self/clear_refs} (see \code{\link{.reset_max_rss}}) + can reset it, so the reading reflects only recent work. The function + pulls the number out of a line that looks like + \code{"VmHWM: 123456 kB"}. + \item \strong{Everything else} (macOS, Windows, or Linux without + \code{/proc}): fall back to the compiled \code{peak_rss_mb()} in + \code{src/peak_rss.cpp}, which asks the operating system directly. +} +} +\keyword{internal} diff --git a/man/dot-reset_max_rss.Rd b/man/dot-reset_max_rss.Rd new file mode 100644 index 00000000..41826fd2 --- /dev/null +++ b/man/dot-reset_max_rss.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R +\name{.reset_max_rss} +\alias{.reset_max_rss} +\title{Reset the process maximum-RSS high-water mark to the current RSS} +\usage{ +.reset_max_rss() +} +\value{ +invisible \code{TRUE} if the high-water mark was reset, + \code{FALSE} otherwise +} +\description{ +The kernel-tracked maximum RSS (\code{VmHWM} on Linux, \code{ru_maxrss} on +macOS/BSD, \code{PeakWorkingSetSize} on Windows) is a running maximum +since process start. Without a reset, an earlier one-off allocation in the +same R session (e.g. reading an 8 GB CSV) inflates every maximum-memory +measurement taken afterwards even if that memory has since been freed, +making it useless as a baseline for benchmarking this function's own +memory use. +} +\details{ +Linux only: writing \code{"5"} to \code{/proc/self/clear_refs} resets +\code{VmHWM} to the current RSS (kernel >= 4.0). macOS and Windows have no +equivalent public facility to reset their maximum-RSS counters, so this is a +no-op there, and \code{track_memory} reports on those platforms may still +include memory used before this function was entered. +} +\keyword{internal} diff --git a/src/peak_rss.cpp b/src/peak_rss.cpp index 1dcd9277..35bcf7fc 100644 --- a/src/peak_rss.cpp +++ b/src/peak_rss.cpp @@ -8,6 +8,17 @@ using namespace Rcpp; #include #endif +// Goal: Ask the operating system for the most memory (RAM) this process has +// ever used. It is a running maximum, not the amount in use right now, +// so freeing memory does not lower it. Takes no arguments. +// +// Output: One number: maximum memory used, in MB. NA if the platform is not +// supported, or if the Windows query fails. +// +// The three branches below do the same thing with different OS calls; the unit +// each one returns differs (Windows and macOS give bytes, Linux gives kB), +// which is why the divisors are not identical. +// // [[Rcpp::export]] double peak_rss_mb() { #if defined(_WIN32) From 20ff1079ae9154f312d9b52fedeadb5f54601413 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Tue, 1 Sep 2026 09:56:41 -0400 Subject: [PATCH 18/18] remove RD files for dot functions --- R/MSstatsSummarizeWithMultipleCores.R | 5 +++++ man/dot-build_summarize_worker.Rd | 20 ------------------ man/dot-max_rss_mb.Rd | 30 --------------------------- man/dot-pack_protein_slot.Rd | 22 -------------------- man/dot-reset_max_rss.Rd | 29 -------------------------- man/dot-unpack_protein_slot.Rd | 21 ------------------- 6 files changed, 5 insertions(+), 122 deletions(-) delete mode 100644 man/dot-build_summarize_worker.Rd delete mode 100644 man/dot-max_rss_mb.Rd delete mode 100644 man/dot-pack_protein_slot.Rd delete mode 100644 man/dot-reset_max_rss.Rd delete mode 100644 man/dot-unpack_protein_slot.Rd diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 7c0736fb..1d4678fa 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -18,6 +18,7 @@ #' #' @return maximum RAM used, in MB #' @keywords internal +#' @noRd .max_rss_mb <- function() { if (.Platform$OS.type != "windows" && file.exists("/proc/self/status")) { ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE) @@ -45,6 +46,7 @@ #' @return invisible \code{TRUE} if the high-water mark was reset, #' \code{FALSE} otherwise #' @keywords internal +#' @noRd .reset_max_rss <- function() { if (Sys.info()[["sysname"]] == "Linux" && file.exists("/proc/self/clear_refs")) { reset_ok <- tryCatch({ @@ -102,6 +104,7 @@ #' @param all_runs character vector of all run names in global order #' @return list with elements \code{packed} (double vector) and \code{meta} (list) #' @keywords internal +#' @noRd .pack_protein_slot <- function(protein_dt, slot_index, all_runs) { n_runs <- length(all_runs) @@ -200,6 +203,7 @@ #' @return data.table compatible with \code{MSstatsSummarizeSingleTMP} / #' \code{MSstatsSummarizeSingleLinear} #' @keywords internal +#' @noRd .unpack_protein_slot <- function(packed, meta) { n_feature_labels <- meta$n_feature_labels @@ -271,6 +275,7 @@ #' parameters, not the caller's run-scale objects. #' #' @keywords internal +#' @noRd .build_summarize_worker <- function( use_TMP, impute, censored_symbol, remove50missing, aft_iterations, equal_variance diff --git a/man/dot-build_summarize_worker.Rd b/man/dot-build_summarize_worker.Rd deleted file mode 100644 index f8a526bf..00000000 --- a/man/dot-build_summarize_worker.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.build_summarize_worker} -\alias{.build_summarize_worker} -\title{Build the per-record worker closure for \code{MSstatsSummarizeWithMultipleCores}} -\usage{ -.build_summarize_worker( - use_TMP, - impute, - censored_symbol, - remove50missing, - aft_iterations, - equal_variance -) -} -\description{ -Defined at package top level so the closure only captures the scalar run -parameters, not the caller's run-scale objects. -} -\keyword{internal} diff --git a/man/dot-max_rss_mb.Rd b/man/dot-max_rss_mb.Rd deleted file mode 100644 index db21e5ef..00000000 --- a/man/dot-max_rss_mb.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.max_rss_mb} -\alias{.max_rss_mb} -\title{Maximum RAM (in MB) this R process has ever used} -\usage{ -.max_rss_mb() -} -\value{ -maximum RAM used, in MB -} -\description{ -There are two paths: -} -\details{ -\enumerate{ - \item \strong{Linux} (non-Windows and \code{/proc/self/status} exists): - read the number Linux itself keeps. \code{/proc/self/status} is a - plain-text file of \code{Key: value} lines describing this process; the - line \code{VmHWM:} holds the maximum RAM ever used, in kB. Preferred - because \code{/proc/self/clear_refs} (see \code{\link{.reset_max_rss}}) - can reset it, so the reading reflects only recent work. The function - pulls the number out of a line that looks like - \code{"VmHWM: 123456 kB"}. - \item \strong{Everything else} (macOS, Windows, or Linux without - \code{/proc}): fall back to the compiled \code{peak_rss_mb()} in - \code{src/peak_rss.cpp}, which asks the operating system directly. -} -} -\keyword{internal} diff --git a/man/dot-pack_protein_slot.Rd b/man/dot-pack_protein_slot.Rd deleted file mode 100644 index a9b74d13..00000000 --- a/man/dot-pack_protein_slot.Rd +++ /dev/null @@ -1,22 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.pack_protein_slot} -\alias{.pack_protein_slot} -\title{Pack one protein slot into a double vector plus metadata} -\usage{ -.pack_protein_slot(protein_dt, slot_index, all_runs) -} -\arguments{ -\item{protein_dt}{data.table rows for one protein (or protein x label) slot} - -\item{slot_index}{integer position of this slot in the global protein list} - -\item{all_runs}{character vector of all run names in global order} -} -\value{ -list with elements \code{packed} (double vector) and \code{meta} (list) -} -\description{ -Pack one protein slot into a double vector plus metadata -} -\keyword{internal} diff --git a/man/dot-reset_max_rss.Rd b/man/dot-reset_max_rss.Rd deleted file mode 100644 index 41826fd2..00000000 --- a/man/dot-reset_max_rss.Rd +++ /dev/null @@ -1,29 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.reset_max_rss} -\alias{.reset_max_rss} -\title{Reset the process maximum-RSS high-water mark to the current RSS} -\usage{ -.reset_max_rss() -} -\value{ -invisible \code{TRUE} if the high-water mark was reset, - \code{FALSE} otherwise -} -\description{ -The kernel-tracked maximum RSS (\code{VmHWM} on Linux, \code{ru_maxrss} on -macOS/BSD, \code{PeakWorkingSetSize} on Windows) is a running maximum -since process start. Without a reset, an earlier one-off allocation in the -same R session (e.g. reading an 8 GB CSV) inflates every maximum-memory -measurement taken afterwards even if that memory has since been freed, -making it useless as a baseline for benchmarking this function's own -memory use. -} -\details{ -Linux only: writing \code{"5"} to \code{/proc/self/clear_refs} resets -\code{VmHWM} to the current RSS (kernel >= 4.0). macOS and Windows have no -equivalent public facility to reset their maximum-RSS counters, so this is a -no-op there, and \code{track_memory} reports on those platforms may still -include memory used before this function was entered. -} -\keyword{internal} diff --git a/man/dot-unpack_protein_slot.Rd b/man/dot-unpack_protein_slot.Rd deleted file mode 100644 index d0a11738..00000000 --- a/man/dot-unpack_protein_slot.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R -\name{.unpack_protein_slot} -\alias{.unpack_protein_slot} -\title{Reconstruct a per-protein data.table from a packed double vector} -\usage{ -.unpack_protein_slot(packed, meta) -} -\arguments{ -\item{packed}{double vector produced by \code{.pack_protein_slot}} - -\item{meta}{metadata list from \code{.pack_protein_slot}} -} -\value{ -data.table compatible with \code{MSstatsSummarizeSingleTMP} / - \code{MSstatsSummarizeSingleLinear} -} -\description{ -Reconstruct a per-protein data.table from a packed double vector -} -\keyword{internal}