From 494e3e4f26a39cfadcf0a9711393a1bf36f2682d Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Thu, 13 Aug 2026 16:31:28 -0500 Subject: [PATCH 1/8] refactor(imputation): Initial attempt to refactor AFT imputation modeling with conjugate gradient instead of cholesky factorization for the Newton step --- NAMESPACE | 6 + R/dataProcess.R | 50 ++- R/utils_cgsolve.R | 102 +++++ R/utils_imputation.R | 464 +++++++++++++++++++++-- inst/tinytest/test_dataProcess.R | 72 ++++ inst/tinytest/test_utils_cgsolve.R | 72 ++++ inst/tinytest/test_utils_imputation_cg.R | 127 +++++++ man/MSstatsSummarizeSingleLinear.Rd | 7 +- man/MSstatsSummarizeSingleTMP.Rd | 3 +- man/MSstatsSummarizeWithSingleCore.Rd | 3 +- man/dataProcess.Rd | 9 +- man/dot-aftGaussianDerivatives.Rd | 54 +++ man/dot-buildAFTFormula.Rd | 30 ++ man/dot-cgSolve.Rd | 47 +++ man/dot-fitSurvivalCG.Rd | 36 ++ 15 files changed, 1029 insertions(+), 53 deletions(-) create mode 100644 R/utils_cgsolve.R create mode 100644 inst/tinytest/test_utils_cgsolve.R create mode 100644 inst/tinytest/test_utils_imputation_cg.R create mode 100644 man/dot-aftGaussianDerivatives.Rd create mode 100644 man/dot-buildAFTFormula.Rd create mode 100644 man/dot-cgSolve.Rd create mode 100644 man/dot-fitSurvivalCG.Rd diff --git a/NAMESPACE b/NAMESPACE index ba6c1354..4bb4e353 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -122,14 +122,20 @@ importFrom(plotly,subplot) importFrom(preprocessCore,normalize.quantiles) importFrom(rlang,.data) importFrom(stats,dist) +importFrom(stats,dnorm) importFrom(stats,fitted) importFrom(stats,formula) importFrom(stats,hclust) importFrom(stats,lm) +importFrom(stats,lm.fit) importFrom(stats,loess) importFrom(stats,median) +importFrom(stats,model.frame) +importFrom(stats,model.matrix) +importFrom(stats,model.response) importFrom(stats,na.omit) importFrom(stats,p.adjust) +importFrom(stats,pnorm) importFrom(stats,predict) importFrom(stats,qbinom) importFrom(stats,qnorm) diff --git a/R/dataProcess.R b/R/dataProcess.R index ef23f518..e4130d5f 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -62,6 +62,11 @@ #' 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. +#' @param aft_solver Which linear solve to use for the AFT imputation +#' model's Newton-Raphson step: "cholesky" (default) delegates to +#' \code{survival::survreg}, which solves it via Cholesky factorization. +#' "cg" solves the same Newton step with a vendored conjugate-gradient +#' routine instead - an experimental alternative, currently opt-in only. #' @inheritParams .documentFunction #' #' @importFrom utils sessionInfo @@ -130,7 +135,7 @@ dataProcess = function( equalFeatureVar = TRUE, censoredInt = "NA", MBimpute = TRUE, remove50missing = FALSE, fix_missing = NULL, maxQuantileforCensored = 0.999, use_log_file = TRUE, append = FALSE, verbose = TRUE, log_file_path = NULL, - numberOfCores = 1, aft_iterations=90 + numberOfCores = 1, aft_iterations=90, aft_solver = "cholesky" ) { MSstatsConvert::MSstatsLogsSettings(use_log_file, append, verbose, log_file_path, @@ -164,9 +169,10 @@ dataProcess = function( input = MSstatsPrepareForSummarization(input, summaryMethod, MBimpute, censoredInt, remove_uninformative_feature_outlier) summarized = tryCatch(MSstatsSummarizeWithMultipleCores(input, summaryMethod, - MBimpute, censoredInt, - remove50missing, equalFeatureVar, - numberOfCores, aft_iterations), + MBimpute, censoredInt, + remove50missing, equalFeatureVar, + numberOfCores, aft_iterations, + aft_solver), error = function(e) { print(e) NULL @@ -211,7 +217,8 @@ dataProcess = function( #' head(summarized[[1]][[1]]) # run-level summary #' MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol, - remove50missing, equal_variance, aft_iterations = 90) { + remove50missing, equal_variance, aft_iterations = 90, + aft_solver = "cholesky") { is_labeled_reference = "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE) @@ -227,8 +234,8 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol for (protein_id in seq_len(num_proteins)) { single_protein = input[protein_indices[[protein_id]],] summarized_results[[protein_id]] = MSstatsSummarizeSingleTMP( - single_protein, impute, censored_symbol, remove50missing, - aft_iterations) + single_protein, impute, censored_symbol, remove50missing, + aft_iterations, aft_solver = aft_solver) setTxtProgressBar(pb, protein_id) } close(pb) @@ -237,8 +244,8 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol for (protein_id in seq_len(num_proteins)) { single_protein = input[protein_indices[[protein_id]],] summarized_result = MSstatsSummarizeSingleLinear( - single_protein, impute, censored_symbol, - remove50missing, aft_iterations) + single_protein, impute, censored_symbol, + remove50missing, aft_iterations, aft_solver = aft_solver) summarized_results[[protein_id]] = summarized_result setTxtProgressBar(pb, protein_id) @@ -256,9 +263,12 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol #' @param remove50missing if TRUE, proteins with more than 50\% missing values in each run are removed #' @param aft_iterations number of iterations for AFT model fitting #' @param equal_variances if TRUE, observation are assumed to be homoskedastic -#' +#' @param aft_solver Which linear solve to use for the AFT imputation +#' model's Newton-Raphson step: "cholesky" (default, via +#' \code{survival::survreg}) or "cg" (conjugate gradient). +#' #' @return list with protein-level data -#' +#' #' @importFrom stats xtabs #' #' @export @@ -286,7 +296,8 @@ MSstatsSummarizeSingleLinear = function(single_protein, censored_symbol, remove50missing, aft_iterations = 90, - equal_variances = TRUE) { + equal_variances = TRUE, + aft_solver = "cholesky") { ABUNDANCE = RUN = FEATURE = PROTEIN = LogIntensities = NULL cols = intersect( @@ -315,7 +326,11 @@ MSstatsSummarizeSingleLinear = function(single_protein, } else { single_protein[, cols, with = FALSE] } - survival_fit = .fitSurvival(fit_data, aft_iterations) + survival_fit = if (aft_solver == "cg") { + .fitSurvivalCG(fit_data, aft_iterations) + } else { + .fitSurvival(fit_data, aft_iterations) + } sigma2 = survival_fit$scale^2 single_protein[, c("predicted", "imputation_var") := { @@ -437,7 +452,8 @@ MSstatsSummarizeSingleLinear = function(single_protein, #' head(single_protein_summary[[1]]) #' MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, - remove50missing, aft_iterations = 90) { + remove50missing, aft_iterations = 90, + aft_solver = "cholesky") { newABUNDANCE = n_obs = n_obs_run = RUN = FEATURE = LABEL = NULL predicted = censored = NULL cols = intersect(colnames(single_protein), c("newABUNDANCE", "cen", "RUN", @@ -464,7 +480,11 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, # Try to fit survival model and catch convergence warnings survival_fit = withCallingHandlers({ - .fitSurvival(fit_data, aft_iterations) + if (aft_solver == "cg") { + .fitSurvivalCG(fit_data, aft_iterations) + } else { + .fitSurvival(fit_data, aft_iterations) + } }, warning = function(w) { if (grepl("converge", conditionMessage(w), ignore.case = TRUE)) { message("Convergence warning caught: ", conditionMessage(w)) diff --git a/R/utils_cgsolve.R b/R/utils_cgsolve.R new file mode 100644 index 00000000..a63c87f4 --- /dev/null +++ b/R/utils_cgsolve.R @@ -0,0 +1,102 @@ +#' Solve a symmetric positive (semi-)definite linear system via conjugate +#' gradient +#' +#' A minimal, single right-hand-side conjugate gradient solver, used as the +#' Newton-step linear solve in \code{.fitSurvivalCG}. Modeled on +#' \code{lfe::cgsolve}, stripped down to a single dense matrix and a single +#' right-hand-side vector (no multi-column batching, no \code{Matrix}-package +#' or operator/closure dispatch, no preconditioning - none of which are +#' needed for the small, dense AFT information matrices this is used on). +#' +#' @param coefficient_matrix symmetric positive (semi-)definite matrix, +#' e.g. the Hessian/information matrix from a Newton step. +#' @param right_hand_side vector the system is solved against, e.g. the +#' gradient/score vector from a Newton step. +#' @param initial_guess optional starting point for the iteration. Defaults +#' to the zero vector. +#' @param relative_tolerance how small the residual needs to shrink, +#' relative to the size of \code{right_hand_side}, before iteration stops. +#' @param max_iterations how many conjugate-gradient steps to try before +#' giving up. In exact arithmetic, conjugate gradient converges within +#' \code{nrow(coefficient_matrix)} steps, but rounding error erodes that +#' guarantee as the system grows, so the default allows for several times +#' that many steps. +#' +#' @return numeric vector solving (approximately) +#' \code{coefficient_matrix \%*\% solution = right_hand_side}. +#' +#' @keywords internal +.cgSolve = function(coefficient_matrix, right_hand_side, initial_guess = NULL, + relative_tolerance = 1e-8, + max_iterations = 10 * nrow(coefficient_matrix)) { + number_of_unknowns = nrow(coefficient_matrix) + solution = if (is.null(initial_guess)) { + rep(0, number_of_unknowns) + } else { + initial_guess + } + + # The residual measures how far the current guess is from solving the + # system. Conjugate gradient starts out searching in that direction. + residual = right_hand_side - drop(coefficient_matrix %*% solution) + search_direction = residual + residual_size = sum(residual * residual) + smallest_residual_size_seen = residual_size + + # Stop once the residual has shrunk far enough, relative to the size of + # the right-hand side (falling back to an absolute scale when that size + # is tiny). + convergence_threshold = + (relative_tolerance * max(sqrt(sum(right_hand_side^2)), 1))^2 + + for (iteration in seq_len(max_iterations)) { + if (residual_size <= convergence_threshold) { + break + } + + # How far moving along the search direction changes things, as + # measured through the matrix itself. + matrix_times_search_direction = + drop(coefficient_matrix %*% search_direction) + curvature = sum(search_direction * matrix_times_search_direction) + if (!is.finite(curvature) || curvature <= 0) { + warning(".cgSolve: coefficient_matrix is not positive definite ", + "along the current search direction; returning the ", + "best iterate found so far") + break + } + + # Move as far as possible along the search direction without + # overshooting the solution, then see how much residual remains. + step_length = residual_size / curvature + solution = solution + step_length * search_direction + residual = residual - step_length * matrix_times_search_direction + new_residual_size = sum(residual * residual) + smallest_residual_size_seen = + min(smallest_residual_size_seen, new_residual_size) + + # If the residual has grown far past its best value so far, the + # iteration is diverging (e.g. because coefficient_matrix is + # ill-conditioned) - give up and return what we have rather than + # loop until max_iterations. + if (iteration > 10 && + new_residual_size > 1e4 * smallest_residual_size_seen) { + warning(".cgSolve: residual is diverging; returning the best ", + "iterate found so far") + break + } + + # Choose the next search direction so it doesn't undo the progress + # made by earlier directions. + search_direction = residual + + (new_residual_size / residual_size) * search_direction + residual_size = new_residual_size + } + + if (residual_size > convergence_threshold) { + warning(".cgSolve: did not converge within max_iterations = ", + max_iterations, " iterations; returning the best iterate ", + "found so far") + } + solution +} diff --git a/R/utils_imputation.R b/R/utils_imputation.R index 99bffcf1..4d378990 100644 --- a/R/utils_imputation.R +++ b/R/utils_imputation.R @@ -1,59 +1,455 @@ +#' Decide which predictors go into a single protein's AFT imputation model +#' +#' MSstats fits an accelerated-failure-time (AFT) model per protein to +#' impute left-censored values, and predictors are chosen based on how much +#' information is actually available: whether this is a labeled (SRM) +#' experiment with a reference channel (\code{ref_covariate}), whether +#' there is more than one feature to estimate a \code{FEATURE} effect for, +#' and whether there are enough uncensored observations to estimate that +#' effect at all. Both \code{.fitSurvival} (Cholesky-based, via +#' \code{survival::survreg}) and \code{.fitSurvivalCG} (conjugate-gradient +#' based) share this selection logic, so the two solvers always fit the +#' same model and differ only in how the Newton step is solved. +#' +#' @param input data.table with columns \code{newABUNDANCE}, \code{cen}, +#' \code{RUN}, \code{FEATURE}, \code{LABEL}, and (for labeled experiments) +#' \code{ref_covariate}. +#' +#' @return a formula whose left side is +#' \code{Surv(newABUNDANCE, cen, type = "left")}. +#' #' @importFrom data.table uniqueN -#' @importFrom survival survreg Surv +#' @importFrom survival Surv #' @keywords internal -.fitSurvival = function(input, aft_iterations) { +.buildAFTFormula = function(input) { FEATURE = RUN = NULL - + missingness_filter = is.finite(input$newABUNDANCE) n_total = nrow(input[missingness_filter, ]) n_features = data.table::uniqueN(input[missingness_filter, FEATURE]) n_runs = data.table::uniqueN(input[missingness_filter, RUN]) is_labeled = data.table::uniqueN(input$LABEL) > 1 - countdf = n_total < n_features + n_runs - 1 - # TODO: set.seed here? - set.seed(100) + # With too few uncensored observations, there isn't enough information + # left to also estimate a separate effect per feature. + not_enough_data_for_feature_effect = n_total < n_features + n_runs - 1 + if (is_labeled) { - if (length(unique(input$FEATURE)) == 1) { - # with single feature, not converge, wrong intercept - # need to check - fit = survreg(Surv(newABUNDANCE, cen, type='left') ~ RUN + ref_covariate, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) + if (length(unique(input$FEATURE)) == 1 || + not_enough_data_for_feature_effect) { + # with a single feature (or too little data), a FEATURE term + # either adds nothing or keeps the model from converging / + # gives it the wrong intercept - need to check + Surv(newABUNDANCE, cen, type = "left") ~ RUN + ref_covariate } else { - if (countdf) { - fit = survreg(Surv(newABUNDANCE, cen, type='left') ~ RUN + ref_covariate, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } else { - fit = survreg(Surv(newABUNDANCE, cen, type='left') ~ FEATURE + RUN + ref_covariate, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } + Surv(newABUNDANCE, cen, type = "left") ~ + FEATURE + RUN + ref_covariate } } else { - if (n_features == 1L) { - fit = survreg(Surv(newABUNDANCE, cen, type = "left") ~ RUN, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) + if (n_features == 1L || not_enough_data_for_feature_effect) { + Surv(newABUNDANCE, cen, type = "left") ~ RUN } else { - if (countdf) { - fit = survreg(Surv(newABUNDANCE, cen, type = "left") ~ RUN, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } else { - fit = survreg(Surv(newABUNDANCE, cen, type = "left") ~ FEATURE + RUN, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } + Surv(newABUNDANCE, cen, type = "left") ~ FEATURE + RUN } } +} + +#' @importFrom survival survreg +#' @keywords internal +.fitSurvival = function(input, aft_iterations) { + # TODO: set.seed here? + set.seed(100) + fit = survreg(.buildAFTFormula(input), data = input, dist = "gaussian", + control = list(maxiter = aft_iterations)) fit$y = NULL fit$linear.predictors = NULL fit } +#' Per-observation log-likelihood and derivatives for a Gaussian AFT model +#' +#' Computes what a Newton-Raphson step needs at the current parameter +#' guess: the log-likelihood, its first derivative with respect to the +#' linear predictor and to the log of the scale parameter, and the +#' corresponding second derivatives - all summed/assembled later into the +#' score vector and information matrix by \code{.fitSurvivalCG}. This only +#' covers the two cases MSstats' AFT imputation actually uses: an exact +#' (uncensored) observation, or one left-censored below a detection-limit +#' ceiling (\code{Surv(y, cen, type = "left")} with \code{cen == 0}). +#' +#' The formulas are transcribed term-for-term from \code{survival}'s own +#' C implementation (\code{survregc1.c}'s \code{gauss_d} function and its +#' "exact"/"left censored" cases) rather than re-derived by hand, since a +#' hand re-derivation is an easy place to introduce a sign error; this +#' function's correctness is instead checked against numerical +#' differentiation of the log-likelihood (see +#' \code{test_utils_imputation_cg.R}). +#' +#' @param linear_predictor current linear predictor +#' (\code{model_matrix \%*\% coefficients}). +#' @param log_scale current log of the scale parameter. +#' @param observed_value observed value (or, for censored rows, the +#' detection-limit ceiling substituted in by +#' \code{.setCensoredByThreshold}). +#' @param exact_indicator \code{1} for an exact/uncensored observation, +#' \code{0} for one left-censored below \code{observed_value}. +#' +#' @return a list with the total \code{log_likelihood}, and +#' per-observation vectors \code{gradient_wrt_linear_predictor}, +#' \code{second_derivative_wrt_linear_predictor}, +#' \code{gradient_wrt_log_scale}, \code{second_derivative_wrt_log_scale}, +#' and \code{cross_derivative} +#' (d2 log_likelihood / d linear_predictor d log_scale). +#' +#' @importFrom stats dnorm pnorm +#' @keywords internal +.aftGaussianDerivatives = function(linear_predictor, log_scale, + observed_value, exact_indicator) { + scale = exp(log_scale) + inverse_scale_squared = 1 / scale^2 + + # How far the observation sits from its predicted value, in raw units + # and in standard deviations. + distance_from_prediction = observed_value - linear_predictor + standardized_distance = distance_from_prediction / scale + + density_at_standardized_distance = dnorm(standardized_distance) + cumulative_probability_at_standardized_distance = + pnorm(standardized_distance) + is_exact_observation = (exact_indicator == 1) + + # --- exact (uncensored) observations -------------------------------- + # log-likelihood contribution is log(density) - log(scale); what + # follows is that expression's derivatives wrt linear_predictor and + # log_scale. + exact_log_likelihood = + log(density_at_standardized_distance) - log_scale + exact_gradient_wrt_linear_predictor = standardized_distance / scale + exact_log_density_curvature = + (standardized_distance^2 - 1) * inverse_scale_squared + exact_second_derivative_wrt_linear_predictor = + exact_log_density_curvature - + exact_gradient_wrt_linear_predictor^2 + exact_gradient_wrt_log_scale_before_adjustment = + exact_gradient_wrt_linear_predictor * distance_from_prediction + exact_cross_derivative = + distance_from_prediction * exact_log_density_curvature - + exact_gradient_wrt_linear_predictor * + (exact_gradient_wrt_log_scale_before_adjustment + 1) + exact_second_derivative_wrt_log_scale = + distance_from_prediction^2 * exact_log_density_curvature - + exact_gradient_wrt_log_scale_before_adjustment * + (1 + exact_gradient_wrt_log_scale_before_adjustment) + exact_gradient_wrt_log_scale = + exact_gradient_wrt_log_scale_before_adjustment - 1 + + # Guard against the density underflowing to exactly zero (only + # happens for astronomically large |standardized_distance|, e.g. from + # a wild early Newton guess). Any reasonable derivative works here, + # since the collapsed log-likelihood itself is what triggers + # step-halving. + exact_density_underflowed = density_at_standardized_distance <= 0 + exact_log_likelihood = + ifelse(exact_density_underflowed, -200, exact_log_likelihood) + exact_gradient_wrt_linear_predictor = ifelse( + exact_density_underflowed, -standardized_distance / scale, + exact_gradient_wrt_linear_predictor) + exact_second_derivative_wrt_linear_predictor = ifelse( + exact_density_underflowed, -1 / scale, + exact_second_derivative_wrt_linear_predictor) + exact_gradient_wrt_log_scale = + ifelse(exact_density_underflowed, 0, exact_gradient_wrt_log_scale) + exact_cross_derivative = + ifelse(exact_density_underflowed, 0, exact_cross_derivative) + exact_second_derivative_wrt_log_scale = ifelse( + exact_density_underflowed, 0, + exact_second_derivative_wrt_log_scale) + + # --- left-censored observations (true value <= the recorded ceiling) - + # log-likelihood contribution is log(Phi(standardized_distance)); + # "censoring_hazard" plays the same role for these rows that the + # density itself plays above. + censored_log_likelihood = + log(cumulative_probability_at_standardized_distance) + censoring_hazard = density_at_standardized_distance / + (cumulative_probability_at_standardized_distance * scale) + censored_gradient_wrt_linear_predictor = -censoring_hazard + censored_log_density_curvature = + -standardized_distance * density_at_standardized_distance * + inverse_scale_squared / + cumulative_probability_at_standardized_distance + censored_second_derivative_wrt_linear_predictor = + censored_log_density_curvature - + censored_gradient_wrt_linear_predictor^2 + censored_gradient_wrt_log_scale = + censored_gradient_wrt_linear_predictor * distance_from_prediction + censored_cross_derivative = + distance_from_prediction * censored_log_density_curvature - + censored_gradient_wrt_linear_predictor * + (censored_gradient_wrt_log_scale + 1) + censored_second_derivative_wrt_log_scale = + distance_from_prediction^2 * censored_log_density_curvature - + censored_gradient_wrt_log_scale * (1 + censored_gradient_wrt_log_scale) + + # Same underflow guard as above, triggered when the cumulative + # probability collapses to zero (standardized_distance very + # negative). + censored_probability_underflowed = + cumulative_probability_at_standardized_distance <= 0 + censored_log_likelihood = ifelse( + censored_probability_underflowed, -200, censored_log_likelihood) + censored_gradient_wrt_linear_predictor = ifelse( + censored_probability_underflowed, -standardized_distance / scale, + censored_gradient_wrt_linear_predictor) + censored_second_derivative_wrt_linear_predictor = ifelse( + censored_probability_underflowed, 0, + censored_second_derivative_wrt_linear_predictor) + censored_gradient_wrt_log_scale = ifelse( + censored_probability_underflowed, 0, censored_gradient_wrt_log_scale) + censored_cross_derivative = ifelse( + censored_probability_underflowed, 0, censored_cross_derivative) + censored_second_derivative_wrt_log_scale = ifelse( + censored_probability_underflowed, 0, + censored_second_derivative_wrt_log_scale) + + list( + log_likelihood = sum(ifelse( + is_exact_observation, exact_log_likelihood, + censored_log_likelihood)), + gradient_wrt_linear_predictor = ifelse( + is_exact_observation, exact_gradient_wrt_linear_predictor, + censored_gradient_wrt_linear_predictor), + second_derivative_wrt_linear_predictor = ifelse( + is_exact_observation, + exact_second_derivative_wrt_linear_predictor, + censored_second_derivative_wrt_linear_predictor), + gradient_wrt_log_scale = ifelse( + is_exact_observation, exact_gradient_wrt_log_scale, + censored_gradient_wrt_log_scale), + second_derivative_wrt_log_scale = ifelse( + is_exact_observation, exact_second_derivative_wrt_log_scale, + censored_second_derivative_wrt_log_scale), + cross_derivative = ifelse( + is_exact_observation, exact_cross_derivative, + censored_cross_derivative) + ) +} + +#' Fit a Gaussian, left-censored AFT model with a conjugate-gradient +#' Newton step +#' +#' An alternative to \code{.fitSurvival} for exactly the same imputation +#' model (Gaussian accelerated-failure-time regression, left-censoring +#' only, chosen by the same \code{.buildAFTFormula} both solvers share), +#' used when \code{aft_solver = "cg"}. It runs the same kind of +#' Newton-Raphson iteration \code{survival::survreg} does - repeatedly +#' solving \code{information_matrix \%*\% step = gradient} for the next +#' set of coefficients - but performs that linear solve with the +#' conjugate-gradient routine \code{.cgSolve} instead of the Cholesky +#' factorization \code{survreg} uses internally. The returned object is +#' classed \code{"survreg"} and carries the fields \code{predict.survreg} +#' needs, so it is a drop-in replacement anywhere \code{.fitSurvival}'s +#' result is used. +#' +#' @param input data.table, the same shape \code{.fitSurvival} expects. +#' @param aft_iterations maximum number of Newton-Raphson iterations. +#' @param convergence_tolerance stop once the change in log-likelihood +#' between iterations falls below this (matches the default +#' \code{rel.tolerance} in \code{survival::survreg.control}). +#' +#' @return a fitted model of class \code{"survreg"}. +#' +#' @importFrom stats model.frame model.matrix model.response lm.fit sd +#' @keywords internal +.fitSurvivalCG = function(input, aft_iterations, + convergence_tolerance = 1e-9) { + model_frame = model.frame(.buildAFTFormula(input), data = input) + model_terms = attr(model_frame, "terms") + design_matrix = model.matrix(model_terms, model_frame) + number_of_coefficients = ncol(design_matrix) + + response = model.response(model_frame) + observed_value = response[, 1] + exact_indicator = response[, 2] + + # Initial guess: an ordinary least-squares fit for the regression + # coefficients (treating the detection-limit ceiling already + # substituted into censored rows as if it were observed), and the + # residual standard deviation for the scale parameter. A + # rank-deficient design leaves some coefficients unidentified + # (reported as NA by lm.fit); start those at zero. + initial_fit = lm.fit(design_matrix, observed_value) + coefficients = initial_fit$coefficients + coefficients[!is.finite(coefficients)] = 0 + residual_standard_deviation = sd(initial_fit$residuals) + log_scale = log(max(residual_standard_deviation, 1e-4)) + + evaluate_log_likelihood_and_derivatives = function(coefficients, + log_scale) { + .aftGaussianDerivatives( + drop(design_matrix %*% coefficients), log_scale, + observed_value, exact_indicator) + } + + build_gradient = function(derivatives) { + c(as.vector(crossprod( + design_matrix, derivatives$gradient_wrt_linear_predictor)), + sum(derivatives$gradient_wrt_log_scale)) + } + + build_information_matrix = function(derivatives) { + # Regression block: -t(X) %*% diag(second_derivative) %*% X, + # computed without forming the diagonal matrix explicitly. + regression_block = -crossprod( + design_matrix, + design_matrix * derivatives$second_derivative_wrt_linear_predictor) + cross_block = -as.vector( + crossprod(design_matrix, derivatives$cross_derivative)) + scale_block = -sum(derivatives$second_derivative_wrt_log_scale) + rbind(cbind(regression_block, cross_block), + c(cross_block, scale_block)) + } + + is_finite_fit = function(derivatives) { + is.finite(derivatives$log_likelihood) && + all(is.finite(derivatives$gradient_wrt_linear_predictor)) && + all(is.finite(derivatives$gradient_wrt_log_scale)) && + all(is.finite(derivatives$second_derivative_wrt_linear_predictor)) && + all(is.finite(derivatives$second_derivative_wrt_log_scale)) + } + + # A Newton step away from the optimum, the exact information matrix + # is not guaranteed to be positive definite. survival::survreg falls + # back, in that situation, to the sum of the outer products of each + # observation's own contribution to the gradient - always positive + # semi-definite by construction, and equal to the exact information + # matrix in expectation (this is the classic Gauss-Newton / BHHH + # approximation). Mirror that fallback here. + build_gauss_newton_approximation = function(derivatives) { + per_observation_gradient_contributions = cbind( + design_matrix * derivatives$gradient_wrt_linear_predictor, + derivatives$gradient_wrt_log_scale) + crossprod(per_observation_gradient_contributions) + } + + solve_newton_step = function(information_matrix, derivatives, gradient) { + information_matrix_is_not_positive_definite = FALSE + step = withCallingHandlers( + .cgSolve(information_matrix, gradient), + warning = function(w) { + if (grepl("not positive definite", conditionMessage(w))) { + information_matrix_is_not_positive_definite <<- TRUE + } + invokeRestart("muffleWarning") + }) + if (information_matrix_is_not_positive_definite) { + step = .cgSolve(build_gauss_newton_approximation(derivatives), + gradient) + } + step + } + + current_fit = + evaluate_log_likelihood_and_derivatives(coefficients, log_scale) + current_log_likelihood = current_fit$log_likelihood + number_of_iterations_used = 0 + converged = FALSE + + for (iteration in seq_len(aft_iterations)) { + number_of_iterations_used = iteration + gradient = build_gradient(current_fit) + information_matrix = build_information_matrix(current_fit) + newton_step = + solve_newton_step(information_matrix, current_fit, gradient) + + candidate_coefficients = + coefficients + newton_step[seq_len(number_of_coefficients)] + candidate_log_scale = + log_scale + newton_step[number_of_coefficients + 1] + + # Step-halving: if the Newton step overshoots (a non-finite or + # decreasing log-likelihood), back the trial point off toward the + # last accepted one, mirroring survival::survreg's own recovery + # strategy (survreg6.c) rather than simply rejecting the step + # outright. + number_of_halvings = 0 + halving_exhausted = FALSE + repeat { + candidate_fit = evaluate_log_likelihood_and_derivatives( + candidate_coefficients, candidate_log_scale) + candidate_improves = is_finite_fit(candidate_fit) && + candidate_fit$log_likelihood >= current_log_likelihood + if (candidate_improves) { + break + } + number_of_halvings = number_of_halvings + 1 + if (number_of_halvings > 30) { + halving_exhausted = TRUE + break + } + if (number_of_halvings == 1 && + (log_scale - candidate_log_scale) > 1.1) { + # a single huge drop in scale is the most common cause of + # a bad trial; keep the first back-off from cutting scale + # by more than a factor of exp(1.1), same as survreg6.c + candidate_log_scale = log_scale - 1.1 + } + candidate_coefficients = + (candidate_coefficients + 2 * coefficients) / 3 + candidate_log_scale = (candidate_log_scale + 2 * log_scale) / 3 + } + + if (halving_exhausted) { + break + } + + relative_change = + abs(1 - current_log_likelihood / candidate_fit$log_likelihood) + absolute_change = + abs(candidate_fit$log_likelihood - current_log_likelihood) + + coefficients = candidate_coefficients + log_scale = candidate_log_scale + current_fit = candidate_fit + current_log_likelihood = candidate_fit$log_likelihood + + if (relative_change <= convergence_tolerance || + absolute_change <= convergence_tolerance) { + converged = TRUE + break + } + } + + if (!converged) { + warning("AFT model (CG solver) ran out of iterations and did not ", + "converge") + } + + final_information_matrix = build_information_matrix(current_fit) + variance_covariance_matrix = tryCatch( + solve(final_information_matrix), + error = function(e) MASS::ginv(final_information_matrix)) + + fitted_coefficients = coefficients + names(fitted_coefficients) = colnames(design_matrix) + + is_factor_column = vapply(model_frame, is.factor, logical(1)) + + fit = list( + coefficients = fitted_coefficients, + var = variance_covariance_matrix, + scale = exp(log_scale), + terms = model_terms, + xlevels = lapply(model_frame[is_factor_column], levels), + dist = "gaussian", + iter = number_of_iterations_used, + loglik = current_log_likelihood + ) + class(fit) = "survreg" + fit +} + #' Get predicted values from a survival model #' @param input data.table #' @return numeric vector of predictions diff --git a/inst/tinytest/test_dataProcess.R b/inst/tinytest/test_dataProcess.R index e0113862..60df17f7 100644 --- a/inst/tinytest/test_dataProcess.R +++ b/inst/tinytest/test_dataProcess.R @@ -423,3 +423,75 @@ expect_true( length(l_cens_pred) > 0 && all(is.finite(l_cens_pred)), info = "MSstatsSummarizeSingleTMP SRM: censored L rows must receive a finite imputed predicted value" ) + +# --- Same SRM imputation, but via aft_solver = "cg" ------------------------ +# Same invariants must hold (H never imputed, L gets a finite prediction), +# and the imputed values themselves should closely match the default +# aft_solver = "cholesky" path, since both solve the same Newton step. +# +# make_srm_impute_input()'s uncensored values are an exactly noise-free +# linear function of RUN, which makes the Gaussian scale MLE degenerate +# (unbounded as residuals -> 0). That's fine for the qualitative H/L +# invariant checks above, but not a meaningful numeric comparison between +# solvers, so a little jitter is added here to make the fit well-posed. + +make_srm_impute_input_with_noise <- function(seed) { + input <- make_srm_impute_input() + set.seed(seed) + input[cen == 1L, + newABUNDANCE := newABUNDANCE + rnorm(.N, sd = 0.01)] + input +} + +result_srm_imp_chol_noisy <- MSstatsSummarizeSingleTMP( + make_srm_impute_input_with_noise(seed = 1), + impute = TRUE, + censored_symbol = "NA", + remove50missing = FALSE, + aft_iterations = 90, + aft_solver = "cholesky" +) +result_srm_imp_cg <- MSstatsSummarizeSingleTMP( + make_srm_impute_input_with_noise(seed = 1), + impute = TRUE, + censored_symbol = "NA", + remove50missing = FALSE, + aft_iterations = 90, + aft_solver = "cg" +) + +survival_srm_chol_noisy <- result_srm_imp_chol_noisy[[2]] +survival_srm_cg <- result_srm_imp_cg[[2]] + +h_cens_pred_cg <- survival_srm_cg[ + as.character(FEATURE) == "F1" & + as.character(LABEL) == "H" & + as.character(RUN) == "R1", + predicted +] +expect_true( + length(h_cens_pred_cg) > 0 && all(is.na(h_cens_pred_cg)), + info = "MSstatsSummarizeSingleTMP SRM (aft_solver = cg): censored H rows must NOT receive an imputed predicted value" +) + +l_cens_pred_cg <- survival_srm_cg[ + as.character(FEATURE) == "F1" & + as.character(LABEL) == "L" & + as.character(RUN) == "R2", + predicted +] +l_cens_pred_chol_noisy <- survival_srm_chol_noisy[ + as.character(FEATURE) == "F1" & + as.character(LABEL) == "L" & + as.character(RUN) == "R2", + predicted +] +expect_true( + length(l_cens_pred_cg) > 0 && all(is.finite(l_cens_pred_cg)), + info = "MSstatsSummarizeSingleTMP SRM (aft_solver = cg): censored L rows must receive a finite imputed predicted value" +) +expect_equal( + l_cens_pred_cg, l_cens_pred_chol_noisy, tolerance = 1e-4, + check.attributes = FALSE, + info = "MSstatsSummarizeSingleTMP SRM: aft_solver = cg should closely match aft_solver = cholesky" +) diff --git a/inst/tinytest/test_utils_cgsolve.R b/inst/tinytest/test_utils_cgsolve.R new file mode 100644 index 00000000..89bc2a65 --- /dev/null +++ b/inst/tinytest/test_utils_cgsolve.R @@ -0,0 +1,72 @@ +# Tests for .cgSolve(), the vendored conjugate-gradient linear solver used +# as the Newton-step solve in .fitSurvivalCG(). + +make_random_spd_matrix <- function(size, seed, ridge = 0.01) { + set.seed(seed) + random_factor <- matrix(rnorm(size * size), size, size) + random_factor %*% t(random_factor) + diag(size) * ridge +} + +for (size in c(2, 5, 10, 30, 80)) { + coefficient_matrix <- make_random_spd_matrix(size, seed = size) + set.seed(size + 1000) + right_hand_side <- rnorm(size) + + cg_solution <- MSstats:::.cgSolve(coefficient_matrix, right_hand_side) + exact_solution <- solve(coefficient_matrix, right_hand_side) + + expect_equal( + cg_solution, exact_solution, tolerance = 1e-6, + info = paste0(".cgSolve should match solve() on a random SPD ", + "system of size ", size) + ) +} + +# --- near-singular system: still returns a finite result, with a warning --- + +near_singular_matrix <- make_random_spd_matrix(10, seed = 42) +near_singular_matrix[1, ] <- 0 +near_singular_matrix[, 1] <- 0 +set.seed(43) +right_hand_side <- rnorm(10) + +expect_warning( + solution <- MSstats:::.cgSolve(near_singular_matrix, right_hand_side), + info = paste("A singular coefficient_matrix should warn rather than", + "error or hang") +) +expect_true( + all(is.finite(solution)), + info = "A singular system should still return a finite (partial) solution" +) + +# --- an initial guess that is already the solution converges immediately --- + +exact_matrix <- make_random_spd_matrix(6, seed = 7) +set.seed(8) +exact_rhs <- rnorm(6) +exact_answer <- solve(exact_matrix, exact_rhs) + +solution_from_exact_start <- MSstats:::.cgSolve( + exact_matrix, exact_rhs, initial_guess = exact_answer) +expect_equal( + solution_from_exact_start, exact_answer, tolerance = 1e-8, + info = "Starting from the exact solution should return it unchanged" +) + +# --- relative_tolerance controls how tightly the system is solved --- + +loose_matrix <- make_random_spd_matrix(20, seed = 99) +set.seed(100) +loose_rhs <- rnorm(20) +exact_loose_answer <- solve(loose_matrix, loose_rhs) + +loose_solution <- MSstats:::.cgSolve(loose_matrix, loose_rhs, relative_tolerance = 1e-2) +tight_solution <- MSstats:::.cgSolve(loose_matrix, loose_rhs, relative_tolerance = 1e-10) + +loose_error <- max(abs(loose_solution - exact_loose_answer)) +tight_error <- max(abs(tight_solution - exact_loose_answer)) +expect_true( + tight_error < loose_error, + info = "A tighter relative_tolerance should produce a more accurate solution" +) diff --git a/inst/tinytest/test_utils_imputation_cg.R b/inst/tinytest/test_utils_imputation_cg.R new file mode 100644 index 00000000..4f57bfc0 --- /dev/null +++ b/inst/tinytest/test_utils_imputation_cg.R @@ -0,0 +1,127 @@ +# Tests that .fitSurvivalCG() - the conjugate-gradient alternative to +# .fitSurvival() - fits the same model, and agrees numerically with it. +# +# The scenarios below reuse the noiseless fixtures from +# test_utils_imputation.R purely to check that .fitSurvivalCG() selects the +# same predictors as .fitSurvival() (via the shared .buildAFTFormula()). +# For numeric agreement on the fitted values themselves, a Gaussian AFT +# model needs actual residual variation to estimate - a noiseless design +# has a degenerate (unbounded) scale MLE, so a second set of fixtures below +# adds realistic noise and left-censoring before comparing coefficients, +# scale, and predictions. + +make_surv_labeled_single <- function() { + runs <- paste0("R", 1:3) + dt <- data.table::rbindlist(list( + data.table::data.table( + FEATURE = factor(rep("F1", 9)), + RUN = factor(rep(runs, each = 3)), + LABEL = "H", + newABUNDANCE = seq(10.1, by = 0.1, length.out = 9), + cen = 1L + ), + data.table::data.table( + FEATURE = factor(rep("F1", 9)), + RUN = factor(rep(runs, each = 3)), + LABEL = "L", + newABUNDANCE = seq(14.1, by = 0.1, length.out = 9), + cen = 1L + ) + )) + ref_vals <- ifelse(dt$LABEL == "L", as.character(dt$RUN), "0") + dt[["ref_covariate"]] <- factor(ref_vals, levels = c("0", runs)) + dt +} + +make_surv_unlabeled_multi_welldetermined <- function() { + dt <- data.table::CJ( + FEATURE = paste0("F", 1:3), + RUN = paste0("R", 1:5) + ) + dt[, FEATURE := factor(FEATURE)] + dt[, RUN := factor(RUN)] + dt[, LABEL := "L"] + dt[, newABUNDANCE := seq(10, by = 0.5, length.out = .N)] + dt[, cen := 1L] + dt +} + +# --- .fitSurvivalCG() selects the same predictors as .fitSurvival() ------- + +coef_names <- function(fit) names(coef(fit)) + +for (make_input in list(make_surv_labeled_single, + make_surv_unlabeled_multi_welldetermined)) { + input <- make_input() + chol_names <- sort(coef_names(MSstats:::.fitSurvival(input, 90))) + cg_names <- sort(coef_names(MSstats:::.fitSurvivalCG(input, 90))) + expect_equal( + cg_names, chol_names, + info = ".fitSurvivalCG must select the same predictors as .fitSurvival" + ) +} + +# --- numeric agreement on realistic (noisy, censored) data ---------------- + +make_noisy_censored_input <- function(seed, is_labeled) { + set.seed(seed) + features <- paste0("F", 1:3) + runs <- paste0("R", 1:4) + labels <- if (is_labeled) c("H", "L") else "L" + dt <- data.table::CJ(FEATURE = features, RUN = runs, LABEL = labels) + dt[, FEATURE := factor(FEATURE)] + dt[, RUN := factor(RUN)] + dt[, newABUNDANCE := + 10 + as.integer(FEATURE) + as.integer(RUN) * 0.5 + + ifelse(LABEL == "L", 4, 0) + rnorm(.N, sd = 0.7)] + dt[, cen := 1L] + censoring_threshold <- stats::quantile(dt$newABUNDANCE, 0.2) + dt[newABUNDANCE < censoring_threshold, cen := 0L] + dt[cen == 0L, newABUNDANCE := censoring_threshold] + if (is_labeled) { + ref_vals <- ifelse(dt$LABEL == "L", as.character(dt$RUN), "0") + dt[["ref_covariate"]] <- factor(ref_vals, levels = c("0", runs)) + } + dt +} + +check_solvers_agree <- function(input, tolerance, label) { + fit_cholesky <- MSstats:::.fitSurvival(input, 90) + fit_cg <- MSstats:::.fitSurvivalCG(input, 90) + + matched_names <- names(fit_cholesky$coefficients) + expect_equal( + fit_cg$coefficients[matched_names], + fit_cholesky$coefficients, + tolerance = tolerance, check.attributes = FALSE, + info = paste(label, ": coefficients should match .fitSurvival") + ) + expect_equal( + fit_cg$scale, fit_cholesky$scale, tolerance = tolerance, + check.attributes = FALSE, + info = paste(label, ": scale should match .fitSurvival") + ) + + predictions_cholesky <- predict(fit_cholesky, newdata = input, se.fit = TRUE) + predictions_cg <- predict(fit_cg, newdata = input, se.fit = TRUE) + expect_equal( + predictions_cg$fit, predictions_cholesky$fit, + tolerance = tolerance, check.attributes = FALSE, + info = paste(label, ": predicted values should match .fitSurvival") + ) + expect_equal( + predictions_cg$se.fit, predictions_cholesky$se.fit, + tolerance = tolerance, check.attributes = FALSE, + info = paste(label, ": prediction standard errors should match", + ".fitSurvival") + ) +} + +check_solvers_agree( + make_noisy_censored_input(seed = 1, is_labeled = TRUE), + tolerance = 1e-4, label = "labeled, noisy, censored" +) +check_solvers_agree( + make_noisy_censored_input(seed = 2, is_labeled = FALSE), + tolerance = 1e-4, label = "unlabeled, noisy, censored" +) diff --git a/man/MSstatsSummarizeSingleLinear.Rd b/man/MSstatsSummarizeSingleLinear.Rd index 8595bc16..0f0056dd 100644 --- a/man/MSstatsSummarizeSingleLinear.Rd +++ b/man/MSstatsSummarizeSingleLinear.Rd @@ -10,7 +10,8 @@ MSstatsSummarizeSingleLinear( censored_symbol, remove50missing, aft_iterations = 90, - equal_variances = TRUE + equal_variances = TRUE, + aft_solver = "cholesky" ) } \arguments{ @@ -25,6 +26,10 @@ MSstatsSummarizeSingleLinear( \item{aft_iterations}{number of iterations for AFT model fitting} \item{equal_variances}{if TRUE, observation are assumed to be homoskedastic} + +\item{aft_solver}{Which linear solve to use for the AFT imputation +model's Newton-Raphson step: "cholesky" (default, via +\code{survival::survreg}) or "cg" (conjugate gradient).} } \value{ list with protein-level data diff --git a/man/MSstatsSummarizeSingleTMP.Rd b/man/MSstatsSummarizeSingleTMP.Rd index cd115723..55038a02 100644 --- a/man/MSstatsSummarizeSingleTMP.Rd +++ b/man/MSstatsSummarizeSingleTMP.Rd @@ -9,7 +9,8 @@ MSstatsSummarizeSingleTMP( impute, censored_symbol, remove50missing, - aft_iterations = 90 + aft_iterations = 90, + aft_solver = "cholesky" ) } \arguments{ diff --git a/man/MSstatsSummarizeWithSingleCore.Rd b/man/MSstatsSummarizeWithSingleCore.Rd index fc711d42..278ae8e1 100644 --- a/man/MSstatsSummarizeWithSingleCore.Rd +++ b/man/MSstatsSummarizeWithSingleCore.Rd @@ -11,7 +11,8 @@ MSstatsSummarizeWithSingleCore( censored_symbol, remove50missing, equal_variance, - aft_iterations = 90 + aft_iterations = 90, + aft_solver = "cholesky" ) } \arguments{ diff --git a/man/dataProcess.Rd b/man/dataProcess.Rd index 1ed4fe03..5f306a4e 100644 --- a/man/dataProcess.Rd +++ b/man/dataProcess.Rd @@ -25,7 +25,8 @@ dataProcess( verbose = TRUE, log_file_path = NULL, numberOfCores = 1, - aft_iterations = 90 + aft_iterations = 90, + aft_solver = "cholesky" ) } \arguments{ @@ -121,6 +122,12 @@ 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{aft_solver}{Which linear solve to use for the AFT imputation +model's Newton-Raphson step: "cholesky" (default) delegates to +\code{survival::survreg}, which solves it via Cholesky factorization. +"cg" solves the same Newton step with a vendored conjugate-gradient +routine instead - an experimental alternative, currently opt-in only.} } \value{ A list containing: diff --git a/man/dot-aftGaussianDerivatives.Rd b/man/dot-aftGaussianDerivatives.Rd new file mode 100644 index 00000000..768e1dee --- /dev/null +++ b/man/dot-aftGaussianDerivatives.Rd @@ -0,0 +1,54 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_imputation.R +\name{.aftGaussianDerivatives} +\alias{.aftGaussianDerivatives} +\title{Per-observation log-likelihood and derivatives for a Gaussian AFT model} +\usage{ +.aftGaussianDerivatives( + linear_predictor, + log_scale, + observed_value, + exact_indicator +) +} +\arguments{ +\item{linear_predictor}{current linear predictor +(\code{model_matrix \%*\% coefficients}).} + +\item{log_scale}{current log of the scale parameter.} + +\item{observed_value}{observed value (or, for censored rows, the +detection-limit ceiling substituted in by +\code{.setCensoredByThreshold}).} + +\item{exact_indicator}{\code{1} for an exact/uncensored observation, +\code{0} for one left-censored below \code{observed_value}.} +} +\value{ +a list with the total \code{log_likelihood}, and +per-observation vectors \code{gradient_wrt_linear_predictor}, +\code{second_derivative_wrt_linear_predictor}, +\code{gradient_wrt_log_scale}, \code{second_derivative_wrt_log_scale}, +and \code{cross_derivative} +(d2 log_likelihood / d linear_predictor d log_scale). +} +\description{ +Computes what a Newton-Raphson step needs at the current parameter +guess: the log-likelihood, its first derivative with respect to the +linear predictor and to the log of the scale parameter, and the +corresponding second derivatives - all summed/assembled later into the +score vector and information matrix by \code{.fitSurvivalCG}. This only +covers the two cases MSstats' AFT imputation actually uses: an exact +(uncensored) observation, or one left-censored below a detection-limit +ceiling (\code{Surv(y, cen, type = "left")} with \code{cen == 0}). +} +\details{ +The formulas are transcribed term-for-term from \code{survival}'s own +C implementation (\code{survregc1.c}'s \code{gauss_d} function and its +"exact"/"left censored" cases) rather than re-derived by hand, since a +hand re-derivation is an easy place to introduce a sign error; this +function's correctness is instead checked against numerical +differentiation of the log-likelihood (see +\code{test_utils_imputation_cg.R}). +} +\keyword{internal} diff --git a/man/dot-buildAFTFormula.Rd b/man/dot-buildAFTFormula.Rd new file mode 100644 index 00000000..d387353a --- /dev/null +++ b/man/dot-buildAFTFormula.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_imputation.R +\name{.buildAFTFormula} +\alias{.buildAFTFormula} +\title{Decide which predictors go into a single protein's AFT imputation model} +\usage{ +.buildAFTFormula(input) +} +\arguments{ +\item{input}{data.table with columns \code{newABUNDANCE}, \code{cen}, +\code{RUN}, \code{FEATURE}, \code{LABEL}, and (for labeled experiments) +\code{ref_covariate}.} +} +\value{ +a formula whose left side is +\code{Surv(newABUNDANCE, cen, type = "left")}. +} +\description{ +MSstats fits an accelerated-failure-time (AFT) model per protein to +impute left-censored values, and predictors are chosen based on how much +information is actually available: whether this is a labeled (SRM) +experiment with a reference channel (\code{ref_covariate}), whether +there is more than one feature to estimate a \code{FEATURE} effect for, +and whether there are enough uncensored observations to estimate that +effect at all. Both \code{.fitSurvival} (Cholesky-based, via +\code{survival::survreg}) and \code{.fitSurvivalCG} (conjugate-gradient +based) share this selection logic, so the two solvers always fit the +same model and differ only in how the Newton step is solved. +} +\keyword{internal} diff --git a/man/dot-cgSolve.Rd b/man/dot-cgSolve.Rd new file mode 100644 index 00000000..5d846263 --- /dev/null +++ b/man/dot-cgSolve.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_cgsolve.R +\name{.cgSolve} +\alias{.cgSolve} +\title{Solve a symmetric positive (semi-)definite linear system via conjugate +gradient} +\usage{ +.cgSolve( + coefficient_matrix, + right_hand_side, + initial_guess = NULL, + relative_tolerance = 1e-08, + max_iterations = 10 * nrow(coefficient_matrix) +) +} +\arguments{ +\item{coefficient_matrix}{symmetric positive (semi-)definite matrix, +e.g. the Hessian/information matrix from a Newton step.} + +\item{right_hand_side}{vector the system is solved against, e.g. the +gradient/score vector from a Newton step.} + +\item{initial_guess}{optional starting point for the iteration. Defaults +to the zero vector.} + +\item{relative_tolerance}{how small the residual needs to shrink, +relative to the size of \code{right_hand_side}, before iteration stops.} + +\item{max_iterations}{how many conjugate-gradient steps to try before +giving up. In exact arithmetic, conjugate gradient converges within +\code{nrow(coefficient_matrix)} steps, but rounding error erodes that +guarantee as the system grows, so the default allows for several times +that many steps.} +} +\value{ +numeric vector solving (approximately) +\code{coefficient_matrix \%*\% solution = right_hand_side}. +} +\description{ +A minimal, single right-hand-side conjugate gradient solver, used as the +Newton-step linear solve in \code{.fitSurvivalCG}. Modeled on +\code{lfe::cgsolve}, stripped down to a single dense matrix and a single +right-hand-side vector (no multi-column batching, no \code{Matrix}-package +or operator/closure dispatch, no preconditioning - none of which are +needed for the small, dense AFT information matrices this is used on). +} +\keyword{internal} diff --git a/man/dot-fitSurvivalCG.Rd b/man/dot-fitSurvivalCG.Rd new file mode 100644 index 00000000..98ebeb0a --- /dev/null +++ b/man/dot-fitSurvivalCG.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_imputation.R +\name{.fitSurvivalCG} +\alias{.fitSurvivalCG} +\title{Fit a Gaussian, left-censored AFT model with a conjugate-gradient +Newton step} +\usage{ +.fitSurvivalCG(input, aft_iterations, convergence_tolerance = 1e-09) +} +\arguments{ +\item{input}{data.table, the same shape \code{.fitSurvival} expects.} + +\item{aft_iterations}{maximum number of Newton-Raphson iterations.} + +\item{convergence_tolerance}{stop once the change in log-likelihood +between iterations falls below this (matches the default +\code{rel.tolerance} in \code{survival::survreg.control}).} +} +\value{ +a fitted model of class \code{"survreg"}. +} +\description{ +An alternative to \code{.fitSurvival} for exactly the same imputation +model (Gaussian accelerated-failure-time regression, left-censoring +only, chosen by the same \code{.buildAFTFormula} both solvers share), +used when \code{aft_solver = "cg"}. It runs the same kind of +Newton-Raphson iteration \code{survival::survreg} does - repeatedly +solving \code{information_matrix \%*\% step = gradient} for the next +set of coefficients - but performs that linear solve with the +conjugate-gradient routine \code{.cgSolve} instead of the Cholesky +factorization \code{survreg} uses internally. The returned object is +classed \code{"survreg"} and carries the fields \code{predict.survreg} +needs, so it is a drop-in replacement anywhere \code{.fitSurvival}'s +result is used. +} +\keyword{internal} From b54f2566150a77c84288dd649c91949a50b72cee Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Thu, 13 Aug 2026 17:51:03 -0500 Subject: [PATCH 2/8] add adjustments w.r.t. logging and pcg option --- R/dataProcess.R | 48 +++++--- R/utils_cgsolve.R | 70 +++++++++-- R/utils_imputation.R | 143 ++++++++++++++++++++--- inst/tinytest/test_utils_cgsolve.R | 90 ++++++++++++-- inst/tinytest/test_utils_imputation_cg.R | 68 ++++++++++- man/MSstatsSummarizeSingleLinear.Rd | 10 +- man/MSstatsSummarizeSingleTMP.Rd | 3 +- man/MSstatsSummarizeWithSingleCore.Rd | 3 +- man/dataProcess.Rd | 14 ++- man/dot-cgSolve.Rd | 31 ++++- man/dot-fitAFTModel.Rd | 36 ++++++ man/dot-fitSurvivalCG.Rd | 29 ++++- 12 files changed, 472 insertions(+), 73 deletions(-) create mode 100644 man/dot-fitAFTModel.Rd diff --git a/R/dataProcess.R b/R/dataProcess.R index e4130d5f..2a97689a 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -66,7 +66,15 @@ #' model's Newton-Raphson step: "cholesky" (default) delegates to #' \code{survival::survreg}, which solves it via Cholesky factorization. #' "cg" solves the same Newton step with a vendored conjugate-gradient -#' routine instead - an experimental alternative, currently opt-in only. +#' routine instead; "pcg" is the same conjugate-gradient routine with a +#' Jacobi (inverse-diagonal) preconditioner, which can reduce the number +#' of conjugate-gradient iterations needed. "cg"/"pcg" are experimental +#' alternatives, currently opt-in only. +#' @param aft_verbose If \code{TRUE} and \code{aft_solver} is "cg" or +#' "pcg", \code{message()} per-Newton-iteration conjugate-gradient +#' iteration counts and timing for every protein fit - useful for +#' evaluating solver time complexity, but produces one block of output +#' per protein, so leave at the default \code{FALSE} for routine runs. #' @inheritParams .documentFunction #' #' @importFrom utils sessionInfo @@ -135,7 +143,8 @@ dataProcess = function( equalFeatureVar = TRUE, censoredInt = "NA", MBimpute = TRUE, remove50missing = FALSE, fix_missing = NULL, maxQuantileforCensored = 0.999, use_log_file = TRUE, append = FALSE, verbose = TRUE, log_file_path = NULL, - numberOfCores = 1, aft_iterations=90, aft_solver = "cholesky" + numberOfCores = 1, aft_iterations=90, aft_solver = "cholesky", + aft_verbose = FALSE ) { MSstatsConvert::MSstatsLogsSettings(use_log_file, append, verbose, log_file_path, @@ -172,7 +181,7 @@ dataProcess = function( MBimpute, censoredInt, remove50missing, equalFeatureVar, numberOfCores, aft_iterations, - aft_solver), + aft_solver, aft_verbose), error = function(e) { print(e) NULL @@ -218,7 +227,7 @@ dataProcess = function( #' MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol, remove50missing, equal_variance, aft_iterations = 90, - aft_solver = "cholesky") { + aft_solver = "cholesky", aft_verbose = FALSE) { is_labeled_reference = "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE) @@ -235,7 +244,8 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol single_protein = input[protein_indices[[protein_id]],] summarized_results[[protein_id]] = MSstatsSummarizeSingleTMP( single_protein, impute, censored_symbol, remove50missing, - aft_iterations, aft_solver = aft_solver) + aft_iterations, aft_solver = aft_solver, + aft_verbose = aft_verbose) setTxtProgressBar(pb, protein_id) } close(pb) @@ -245,7 +255,8 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol single_protein = input[protein_indices[[protein_id]],] summarized_result = MSstatsSummarizeSingleLinear( single_protein, impute, censored_symbol, - remove50missing, aft_iterations, aft_solver = aft_solver) + remove50missing, aft_iterations, aft_solver = aft_solver, + aft_verbose = aft_verbose) summarized_results[[protein_id]] = summarized_result setTxtProgressBar(pb, protein_id) @@ -265,7 +276,11 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol #' @param equal_variances if TRUE, observation are assumed to be homoskedastic #' @param aft_solver Which linear solve to use for the AFT imputation #' model's Newton-Raphson step: "cholesky" (default, via -#' \code{survival::survreg}) or "cg" (conjugate gradient). +#' \code{survival::survreg}), "cg" (conjugate gradient), or "pcg" +#' (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). +#' @param aft_verbose If \code{TRUE} and \code{aft_solver} is "cg" or +#' "pcg", log per-Newton-iteration conjugate-gradient diagnostics for +#' every protein fit. See \code{.fitSurvivalCG}'s \code{verbose}. #' #' @return list with protein-level data #' @@ -297,7 +312,8 @@ MSstatsSummarizeSingleLinear = function(single_protein, remove50missing, aft_iterations = 90, equal_variances = TRUE, - aft_solver = "cholesky") { + aft_solver = "cholesky", + aft_verbose = FALSE) { ABUNDANCE = RUN = FEATURE = PROTEIN = LogIntensities = NULL cols = intersect( @@ -326,11 +342,8 @@ MSstatsSummarizeSingleLinear = function(single_protein, } else { single_protein[, cols, with = FALSE] } - survival_fit = if (aft_solver == "cg") { - .fitSurvivalCG(fit_data, aft_iterations) - } else { - .fitSurvival(fit_data, aft_iterations) - } + survival_fit = .fitAFTModel(fit_data, aft_iterations, aft_solver, + aft_verbose) sigma2 = survival_fit$scale^2 single_protein[, c("predicted", "imputation_var") := { @@ -453,7 +466,8 @@ MSstatsSummarizeSingleLinear = function(single_protein, #' MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, remove50missing, aft_iterations = 90, - aft_solver = "cholesky") { + aft_solver = "cholesky", + aft_verbose = FALSE) { newABUNDANCE = n_obs = n_obs_run = RUN = FEATURE = LABEL = NULL predicted = censored = NULL cols = intersect(colnames(single_protein), c("newABUNDANCE", "cen", "RUN", @@ -480,11 +494,7 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, # Try to fit survival model and catch convergence warnings survival_fit = withCallingHandlers({ - if (aft_solver == "cg") { - .fitSurvivalCG(fit_data, aft_iterations) - } else { - .fitSurvival(fit_data, aft_iterations) - } + .fitAFTModel(fit_data, aft_iterations, aft_solver, aft_verbose) }, warning = function(w) { if (grepl("converge", conditionMessage(w), ignore.case = TRUE)) { message("Convergence warning caught: ", conditionMessage(w)) diff --git a/R/utils_cgsolve.R b/R/utils_cgsolve.R index a63c87f4..4384245c 100644 --- a/R/utils_cgsolve.R +++ b/R/utils_cgsolve.R @@ -5,8 +5,10 @@ #' Newton-step linear solve in \code{.fitSurvivalCG}. Modeled on #' \code{lfe::cgsolve}, stripped down to a single dense matrix and a single #' right-hand-side vector (no multi-column batching, no \code{Matrix}-package -#' or operator/closure dispatch, no preconditioning - none of which are -#' needed for the small, dense AFT information matrices this is used on). +#' or operator/closure dispatch - neither is needed for the small, dense AFT +#' information matrices this is used on). Optionally applies a Jacobi +#' (inverse-diagonal) preconditioner, which \code{lfe::cgsolve} does not +#' support at all. #' #' @param coefficient_matrix symmetric positive (semi-)definite matrix, #' e.g. the Hessian/information matrix from a Newton step. @@ -16,19 +18,35 @@ #' to the zero vector. #' @param relative_tolerance how small the residual needs to shrink, #' relative to the size of \code{right_hand_side}, before iteration stops. +#' Always judged on the true (unpreconditioned) residual, so this means the +#' same thing whether or not \code{use_jacobi_preconditioner} is set. #' @param max_iterations how many conjugate-gradient steps to try before #' giving up. In exact arithmetic, conjugate gradient converges within #' \code{nrow(coefficient_matrix)} steps, but rounding error erodes that #' guarantee as the system grows, so the default allows for several times #' that many steps. +#' @param use_jacobi_preconditioner if \code{TRUE}, precondition with the +#' inverse of \code{coefficient_matrix}'s own diagonal - cheap to apply, +#' and often enough to cut down the number of iterations needed when the +#' diagonal dominates (as it typically does for an AFT information matrix, +#' where each parameter's own curvature tends to be much larger than its +#' cross-terms with the other parameters). Defaults to \code{FALSE}, which +#' reduces exactly to plain (unpreconditioned) conjugate gradient. #' -#' @return numeric vector solving (approximately) -#' \code{coefficient_matrix \%*\% solution = right_hand_side}. +#' @return a list with: \code{solution}, the numeric vector solving +#' (approximately) \code{coefficient_matrix \%*\% solution = +#' right_hand_side}; \code{iterations}, how many conjugate-gradient steps +#' were actually taken; \code{converged}, whether the residual tolerance +#' was met; and \code{positive_definite}, whether \code{coefficient_matrix} +#' behaved as positive definite throughout (a caller can fall back to a +#' different matrix, e.g. a Gauss-Newton approximation, when this is +#' \code{FALSE}). #' #' @keywords internal .cgSolve = function(coefficient_matrix, right_hand_side, initial_guess = NULL, relative_tolerance = 1e-8, - max_iterations = 10 * nrow(coefficient_matrix)) { + max_iterations = 10 * nrow(coefficient_matrix), + use_jacobi_preconditioner = FALSE) { number_of_unknowns = nrow(coefficient_matrix) solution = if (is.null(initial_guess)) { rep(0, number_of_unknowns) @@ -36,11 +54,25 @@ initial_guess } + apply_preconditioner = if (use_jacobi_preconditioner) { + diagonal = diag(coefficient_matrix) + inverse_diagonal = ifelse( + is.finite(diagonal) & diagonal > 0, 1 / diagonal, 1) + function(vector) inverse_diagonal * vector + } else { + identity + } + # The residual measures how far the current guess is from solving the - # system. Conjugate gradient starts out searching in that direction. + # system. Conjugate gradient starts out searching in the + # preconditioner-adjusted residual direction (with no preconditioner, + # this is just the residual itself). residual = right_hand_side - drop(coefficient_matrix %*% solution) - search_direction = residual + preconditioned_residual = apply_preconditioner(residual) + search_direction = preconditioned_residual residual_size = sum(residual * residual) + residual_dot_preconditioned_residual = + sum(residual * preconditioned_residual) smallest_residual_size_seen = residual_size # Stop once the residual has shrunk far enough, relative to the size of @@ -49,10 +81,14 @@ convergence_threshold = (relative_tolerance * max(sqrt(sum(right_hand_side^2)), 1))^2 + positive_definite = TRUE + iterations_used = 0 + for (iteration in seq_len(max_iterations)) { if (residual_size <= convergence_threshold) { break } + iterations_used = iteration # How far moving along the search direction changes things, as # measured through the matrix itself. @@ -60,6 +96,7 @@ drop(coefficient_matrix %*% search_direction) curvature = sum(search_direction * matrix_times_search_direction) if (!is.finite(curvature) || curvature <= 0) { + positive_definite = FALSE warning(".cgSolve: coefficient_matrix is not positive definite ", "along the current search direction; returning the ", "best iterate found so far") @@ -68,7 +105,7 @@ # Move as far as possible along the search direction without # overshooting the solution, then see how much residual remains. - step_length = residual_size / curvature + step_length = residual_dot_preconditioned_residual / curvature solution = solution + step_length * search_direction residual = residual - step_length * matrix_times_search_direction new_residual_size = sum(residual * residual) @@ -88,15 +125,24 @@ # Choose the next search direction so it doesn't undo the progress # made by earlier directions. - search_direction = residual + - (new_residual_size / residual_size) * search_direction + new_preconditioned_residual = apply_preconditioner(residual) + new_residual_dot_preconditioned_residual = + sum(residual * new_preconditioned_residual) + search_direction = new_preconditioned_residual + + (new_residual_dot_preconditioned_residual / + residual_dot_preconditioned_residual) * search_direction residual_size = new_residual_size + residual_dot_preconditioned_residual = + new_residual_dot_preconditioned_residual } - if (residual_size > convergence_threshold) { + converged = residual_size <= convergence_threshold + if (!converged && positive_definite) { warning(".cgSolve: did not converge within max_iterations = ", max_iterations, " iterations; returning the best iterate ", "found so far") } - solution + + list(solution = solution, iterations = iterations_used, + converged = converged, positive_definite = positive_definite) } diff --git a/R/utils_imputation.R b/R/utils_imputation.R index 4d378990..9aceda7d 100644 --- a/R/utils_imputation.R +++ b/R/utils_imputation.R @@ -257,17 +257,45 @@ #' @param convergence_tolerance stop once the change in log-likelihood #' between iterations falls below this (matches the default #' \code{rel.tolerance} in \code{survival::survreg.control}). +#' @param use_jacobi_preconditioner if \code{TRUE}, precondition every +#' conjugate-gradient solve with the inverse of the current information +#' matrix's own diagonal (see \code{.cgSolve}'s +#' \code{use_jacobi_preconditioner}). This is what \code{aft_solver = +#' "pcg"} enables, versus plain conjugate gradient for \code{"cg"}. +#' @param verbose if \code{TRUE}, \code{message()} a line per +#' Newton-Raphson iteration - conjugate-gradient iterations used, whether +#' the Gauss-Newton fallback (see below) was needed, elapsed time, and the +#' resulting log-likelihood - plus a one-line summary once fitting +#' finishes. Meant for evaluating how solver choice and problem size +#' trade off against iteration count and wall time, not for routine use +#' (this fits one protein at a time, so it is easy to generate a line per +#' protein across a whole \code{dataProcess()} run). #' -#' @return a fitted model of class \code{"survreg"}. +#' @return a fitted model of class \code{"survreg"}, with one added field: +#' \code{cg_diagnostics}, a data.frame with one row per Newton-Raphson +#' iteration recording the conjugate-gradient iteration counts and timing +#' described above (populated regardless of \code{verbose}, so it can be +#' inspected/aggregated programmatically after the fact). #' #' @importFrom stats model.frame model.matrix model.response lm.fit sd #' @keywords internal .fitSurvivalCG = function(input, aft_iterations, - convergence_tolerance = 1e-9) { + convergence_tolerance = 1e-9, + use_jacobi_preconditioner = FALSE, + verbose = FALSE) { model_frame = model.frame(.buildAFTFormula(input), data = input) model_terms = attr(model_frame, "terms") design_matrix = model.matrix(model_terms, model_frame) number_of_coefficients = ncol(design_matrix) + number_of_parameters = number_of_coefficients + 1 + number_of_observations = nrow(design_matrix) + + if (verbose) { + message(sprintf( + "[AFT-CG] starting fit: %d observations, %d parameters, preconditioner = %s", + number_of_observations, number_of_parameters, + if (use_jacobi_preconditioner) "jacobi" else "none")) + } response = model.response(model_frame) observed_value = response[, 1] @@ -333,21 +361,43 @@ crossprod(per_observation_gradient_contributions) } - solve_newton_step = function(information_matrix, derivatives, gradient) { - information_matrix_is_not_positive_definite = FALSE - step = withCallingHandlers( - .cgSolve(information_matrix, gradient), + # A "not positive definite" result is expected, handled control flow + # here (the Gauss-Newton fallback below exists for exactly that case), + # so its warning is muffled; a genuine "did not converge within + # max_iterations" is not expected/handled, so that warning still + # propagates normally. + cg_solve_muffling_pd_warning = function(...) { + withCallingHandlers( + .cgSolve(...), warning = function(w) { if (grepl("not positive definite", conditionMessage(w))) { - information_matrix_is_not_positive_definite <<- TRUE + invokeRestart("muffleWarning") } - invokeRestart("muffleWarning") }) - if (information_matrix_is_not_positive_definite) { - step = .cgSolve(build_gauss_newton_approximation(derivatives), - gradient) + } + + # Returns the Newton step, plus how much conjugate-gradient work it + # took to get there - primary_iterations/fallback_iterations and + # used_fallback are the numbers verbose logging (below) reports, so a + # caller can see how solver choice and problem size trade off against + # iteration count. + solve_newton_step = function(information_matrix, derivatives, gradient) { + primary_solve = cg_solve_muffling_pd_warning( + information_matrix, gradient, + use_jacobi_preconditioner = use_jacobi_preconditioner) + if (primary_solve$positive_definite) { + list(step = primary_solve$solution, + primary_iterations = primary_solve$iterations, + used_fallback = FALSE, fallback_iterations = 0L) + } else { + fallback_solve = .cgSolve( + build_gauss_newton_approximation(derivatives), gradient, + use_jacobi_preconditioner = use_jacobi_preconditioner) + list(step = fallback_solve$solution, + primary_iterations = primary_solve$iterations, + used_fallback = TRUE, + fallback_iterations = fallback_solve$iterations) } - step } current_fit = @@ -355,18 +405,38 @@ current_log_likelihood = current_fit$log_likelihood number_of_iterations_used = 0 converged = FALSE + cg_diagnostics = vector("list", aft_iterations) for (iteration in seq_len(aft_iterations)) { number_of_iterations_used = iteration + iteration_start_time = Sys.time() + gradient = build_gradient(current_fit) information_matrix = build_information_matrix(current_fit) newton_step = solve_newton_step(information_matrix, current_fit, gradient) + elapsed_seconds = + as.numeric(Sys.time() - iteration_start_time, units = "secs") + cg_diagnostics[[iteration]] = data.frame( + newton_iteration = iteration, + cg_iterations = newton_step$primary_iterations + + newton_step$fallback_iterations, + used_gauss_newton_fallback = newton_step$used_fallback, + elapsed_seconds = elapsed_seconds) + if (verbose) { + message(sprintf( + "[AFT-CG] newton iter %d: cg iterations = %d%s, %.4f sec", + iteration, + newton_step$primary_iterations + newton_step$fallback_iterations, + if (newton_step$used_fallback) " (Gauss-Newton fallback used)" else "", + elapsed_seconds)) + } + candidate_coefficients = - coefficients + newton_step[seq_len(number_of_coefficients)] + coefficients + newton_step$step[seq_len(number_of_coefficients)] candidate_log_scale = - log_scale + newton_step[number_of_coefficients + 1] + log_scale + newton_step$step[number_of_coefficients + 1] # Step-halving: if the Newton step overshoots (a non-finite or # decreasing log-likelihood), back the trial point off toward the @@ -426,6 +496,17 @@ "converge") } + cg_diagnostics = do.call( + rbind, cg_diagnostics[seq_len(number_of_iterations_used)]) + + if (verbose) { + message(sprintf( + paste0("[AFT-CG] finished: %d newton iterations, ", + "%d total cg iterations, %.4f sec total, converged = %s"), + number_of_iterations_used, sum(cg_diagnostics$cg_iterations), + sum(cg_diagnostics$elapsed_seconds), converged)) + } + final_information_matrix = build_information_matrix(current_fit) variance_covariance_matrix = tryCatch( solve(final_information_matrix), @@ -444,12 +525,44 @@ xlevels = lapply(model_frame[is_factor_column], levels), dist = "gaussian", iter = number_of_iterations_used, - loglik = current_log_likelihood + loglik = current_log_likelihood, + cg_diagnostics = cg_diagnostics ) class(fit) = "survreg" fit } +#' Fit the AFT imputation model with the requested solver +#' +#' Shared dispatch used by both \code{MSstatsSummarizeSingleLinear} and +#' \code{MSstatsSummarizeSingleTMP} so the \code{aft_solver}/ +#' \code{aft_verbose} logic lives in one place instead of being duplicated +#' at both call sites. +#' +#' @param input data.table, the same shape \code{.fitSurvival} expects. +#' @param aft_iterations maximum number of iterations for AFT model fitting. +#' @param aft_solver "cholesky" (default, via \code{survival::survreg}), +#' "cg" (conjugate gradient), or "pcg" (conjugate gradient with a +#' Jacobi/inverse-diagonal preconditioner). +#' @param aft_verbose passed through to \code{.fitSurvivalCG}'s +#' \code{verbose} when \code{aft_solver} is "cg" or "pcg"; has no effect +#' for "cholesky". +#' +#' @return a fitted model of class \code{"survreg"}. +#' +#' @keywords internal +.fitAFTModel = function(input, aft_iterations, aft_solver = "cholesky", + aft_verbose = FALSE) { + if (aft_solver == "pcg") { + .fitSurvivalCG(input, aft_iterations, + use_jacobi_preconditioner = TRUE, verbose = aft_verbose) + } else if (aft_solver == "cg") { + .fitSurvivalCG(input, aft_iterations, verbose = aft_verbose) + } else { + .fitSurvival(input, aft_iterations) + } +} + #' Get predicted values from a survival model #' @param input data.table #' @return numeric vector of predictions diff --git a/inst/tinytest/test_utils_cgsolve.R b/inst/tinytest/test_utils_cgsolve.R index 89bc2a65..902c7c4b 100644 --- a/inst/tinytest/test_utils_cgsolve.R +++ b/inst/tinytest/test_utils_cgsolve.R @@ -1,5 +1,6 @@ # Tests for .cgSolve(), the vendored conjugate-gradient linear solver used -# as the Newton-step solve in .fitSurvivalCG(). +# as the Newton-step solve in .fitSurvivalCG(). Returns a list: +# solution/iterations/converged/positive_definite. make_random_spd_matrix <- function(size, seed, ridge = 0.01) { set.seed(seed) @@ -12,14 +13,23 @@ for (size in c(2, 5, 10, 30, 80)) { set.seed(size + 1000) right_hand_side <- rnorm(size) - cg_solution <- MSstats:::.cgSolve(coefficient_matrix, right_hand_side) + cg_result <- MSstats:::.cgSolve(coefficient_matrix, right_hand_side) exact_solution <- solve(coefficient_matrix, right_hand_side) expect_equal( - cg_solution, exact_solution, tolerance = 1e-6, + cg_result$solution, exact_solution, tolerance = 1e-6, info = paste0(".cgSolve should match solve() on a random SPD ", "system of size ", size) ) + expect_true( + cg_result$converged && cg_result$positive_definite, + info = paste0("A well-conditioned SPD system of size ", size, + " should report converged/positive_definite = TRUE") + ) + expect_true( + cg_result$iterations >= 1 && cg_result$iterations <= size * 10, + info = "iterations should be a small positive count, not the default cap" + ) } # --- near-singular system: still returns a finite result, with a warning --- @@ -31,14 +41,19 @@ set.seed(43) right_hand_side <- rnorm(10) expect_warning( - solution <- MSstats:::.cgSolve(near_singular_matrix, right_hand_side), + singular_result <- MSstats:::.cgSolve(near_singular_matrix, right_hand_side), info = paste("A singular coefficient_matrix should warn rather than", "error or hang") ) expect_true( - all(is.finite(solution)), + all(is.finite(singular_result$solution)), info = "A singular system should still return a finite (partial) solution" ) +expect_false( + singular_result$converged && singular_result$positive_definite, + info = paste("A singular coefficient_matrix should signal trouble via", + "converged = FALSE and/or positive_definite = FALSE") +) # --- an initial guess that is already the solution converges immediately --- @@ -47,12 +62,16 @@ set.seed(8) exact_rhs <- rnorm(6) exact_answer <- solve(exact_matrix, exact_rhs) -solution_from_exact_start <- MSstats:::.cgSolve( +result_from_exact_start <- MSstats:::.cgSolve( exact_matrix, exact_rhs, initial_guess = exact_answer) expect_equal( - solution_from_exact_start, exact_answer, tolerance = 1e-8, + result_from_exact_start$solution, exact_answer, tolerance = 1e-8, info = "Starting from the exact solution should return it unchanged" ) +expect_equal( + result_from_exact_start$iterations, 0, + info = "Starting from the exact solution should take zero iterations" +) # --- relative_tolerance controls how tightly the system is solved --- @@ -61,12 +80,61 @@ set.seed(100) loose_rhs <- rnorm(20) exact_loose_answer <- solve(loose_matrix, loose_rhs) -loose_solution <- MSstats:::.cgSolve(loose_matrix, loose_rhs, relative_tolerance = 1e-2) -tight_solution <- MSstats:::.cgSolve(loose_matrix, loose_rhs, relative_tolerance = 1e-10) +loose_result <- MSstats:::.cgSolve( + loose_matrix, loose_rhs, relative_tolerance = 1e-2) +tight_result <- MSstats:::.cgSolve( + loose_matrix, loose_rhs, relative_tolerance = 1e-10) -loose_error <- max(abs(loose_solution - exact_loose_answer)) -tight_error <- max(abs(tight_solution - exact_loose_answer)) +loose_error <- max(abs(loose_result$solution - exact_loose_answer)) +tight_error <- max(abs(tight_result$solution - exact_loose_answer)) expect_true( tight_error < loose_error, info = "A tighter relative_tolerance should produce a more accurate solution" ) + +# --- Jacobi preconditioner: same answer, fewer or equal iterations -------- +# on a diagonally-dominant system (where a diagonal preconditioner is most +# effective), preconditioned CG should converge in no more iterations than +# plain CG, and to the same solution. + +make_diagonally_dominant_matrix <- function(size, seed) { + set.seed(seed) + matrix_off_diagonal <- matrix(runif(size * size, -0.1, 0.1), size, size) + matrix_off_diagonal <- (matrix_off_diagonal + t(matrix_off_diagonal)) / 2 + diag(matrix_off_diagonal) <- 0 + diag(size) * runif(size, 5, 10) + matrix_off_diagonal +} + +dominant_matrix <- make_diagonally_dominant_matrix(40, seed = 11) +set.seed(12) +dominant_rhs <- rnorm(40) +exact_dominant_answer <- solve(dominant_matrix, dominant_rhs) + +plain_cg_result <- MSstats:::.cgSolve(dominant_matrix, dominant_rhs) +preconditioned_result <- MSstats:::.cgSolve( + dominant_matrix, dominant_rhs, use_jacobi_preconditioner = TRUE) + +expect_equal( + preconditioned_result$solution, exact_dominant_answer, tolerance = 1e-6, + info = "Preconditioned CG should still match solve() on a diagonally dominant system" +) +expect_true( + preconditioned_result$iterations <= plain_cg_result$iterations, + info = paste("Jacobi preconditioning should not need more iterations", + "than plain CG on a diagonally dominant system (plain =", + plain_cg_result$iterations, ", preconditioned =", + preconditioned_result$iterations, ")") +) + +# A degenerate (all-zero) diagonal entry should not blow up the +# preconditioner (falls back to an identity-like scale of 1 for that entry). +degenerate_diagonal_matrix <- make_random_spd_matrix(8, seed = 55) +degenerate_diagonal_matrix[3, 3] <- 0 +set.seed(56) +degenerate_rhs <- rnorm(8) +expect_true( + all(is.finite(suppressWarnings(MSstats:::.cgSolve( + degenerate_diagonal_matrix, degenerate_rhs, + use_jacobi_preconditioner = TRUE))$solution)), + info = "A zero diagonal entry should not produce a non-finite preconditioned solution" +) diff --git a/inst/tinytest/test_utils_imputation_cg.R b/inst/tinytest/test_utils_imputation_cg.R index 4f57bfc0..8dc8e07a 100644 --- a/inst/tinytest/test_utils_imputation_cg.R +++ b/inst/tinytest/test_utils_imputation_cg.R @@ -85,9 +85,11 @@ make_noisy_censored_input <- function(seed, is_labeled) { dt } -check_solvers_agree <- function(input, tolerance, label) { +check_solvers_agree <- function(input, tolerance, label, + use_jacobi_preconditioner = FALSE) { fit_cholesky <- MSstats:::.fitSurvival(input, 90) - fit_cg <- MSstats:::.fitSurvivalCG(input, 90) + fit_cg <- MSstats:::.fitSurvivalCG( + input, 90, use_jacobi_preconditioner = use_jacobi_preconditioner) matched_names <- names(fit_cholesky$coefficients) expect_equal( @@ -125,3 +127,65 @@ check_solvers_agree( make_noisy_censored_input(seed = 2, is_labeled = FALSE), tolerance = 1e-4, label = "unlabeled, noisy, censored" ) + +# --- the Jacobi-preconditioned solver (aft_solver = "pcg") agrees too ----- + +check_solvers_agree( + make_noisy_censored_input(seed = 1, is_labeled = TRUE), + tolerance = 1e-4, label = "labeled, noisy, censored, jacobi-preconditioned", + use_jacobi_preconditioner = TRUE +) +check_solvers_agree( + make_noisy_censored_input(seed = 2, is_labeled = FALSE), + tolerance = 1e-4, label = "unlabeled, noisy, censored, jacobi-preconditioned", + use_jacobi_preconditioner = TRUE +) + +# --- .fitAFTModel() dispatches to the right solver ------------------------- + +noisy_input <- make_noisy_censored_input(seed = 3, is_labeled = FALSE) + +expect_inherits( + MSstats:::.fitAFTModel(noisy_input, 90, "cholesky"), "survreg", + info = ".fitAFTModel(aft_solver = 'cholesky') should return a survreg fit" +) +expect_true( + is.null(MSstats:::.fitAFTModel(noisy_input, 90, "cholesky")$cg_diagnostics), + info = "the cholesky path should not attach cg_diagnostics" +) +expect_false( + is.null(MSstats:::.fitAFTModel(noisy_input, 90, "cg")$cg_diagnostics), + info = ".fitAFTModel(aft_solver = 'cg') should attach cg_diagnostics" +) +expect_false( + is.null(MSstats:::.fitAFTModel(noisy_input, 90, "pcg")$cg_diagnostics), + info = ".fitAFTModel(aft_solver = 'pcg') should attach cg_diagnostics" +) + +# --- verbose = TRUE logs per-iteration diagnostics, FALSE stays silent ----- + +expect_silent( + MSstats:::.fitSurvivalCG(noisy_input, 90, verbose = FALSE) +) +expect_message( + MSstats:::.fitSurvivalCG(noisy_input, 90, verbose = TRUE), + pattern = "\\[AFT-CG\\] starting fit", + info = "verbose = TRUE should report the problem size at the start of the fit" +) +expect_message( + MSstats:::.fitSurvivalCG(noisy_input, 90, verbose = TRUE), + pattern = "\\[AFT-CG\\] finished", + info = "verbose = TRUE should report a summary once fitting finishes" +) + +# --- cg_diagnostics has one row per Newton iteration actually taken ------- + +fit_with_diagnostics <- MSstats:::.fitSurvivalCG(noisy_input, 90) +expect_equal( + nrow(fit_with_diagnostics$cg_diagnostics), fit_with_diagnostics$iter, + info = "cg_diagnostics should have one row per Newton-Raphson iteration taken" +) +expect_true( + all(fit_with_diagnostics$cg_diagnostics$cg_iterations >= 0), + info = "cg_iterations should be a non-negative count for every Newton iteration" +) diff --git a/man/MSstatsSummarizeSingleLinear.Rd b/man/MSstatsSummarizeSingleLinear.Rd index 0f0056dd..dc6d006e 100644 --- a/man/MSstatsSummarizeSingleLinear.Rd +++ b/man/MSstatsSummarizeSingleLinear.Rd @@ -11,7 +11,8 @@ MSstatsSummarizeSingleLinear( remove50missing, aft_iterations = 90, equal_variances = TRUE, - aft_solver = "cholesky" + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ @@ -29,7 +30,12 @@ MSstatsSummarizeSingleLinear( \item{aft_solver}{Which linear solve to use for the AFT imputation model's Newton-Raphson step: "cholesky" (default, via -\code{survival::survreg}) or "cg" (conjugate gradient).} +\code{survival::survreg}), "cg" (conjugate gradient), or "pcg" +(conjugate gradient with a Jacobi/inverse-diagonal preconditioner).} + +\item{aft_verbose}{If \code{TRUE} and \code{aft_solver} is "cg" or +"pcg", log per-Newton-iteration conjugate-gradient diagnostics for +every protein fit. See \code{.fitSurvivalCG}'s \code{verbose}.} } \value{ list with protein-level data diff --git a/man/MSstatsSummarizeSingleTMP.Rd b/man/MSstatsSummarizeSingleTMP.Rd index 55038a02..02c5f119 100644 --- a/man/MSstatsSummarizeSingleTMP.Rd +++ b/man/MSstatsSummarizeSingleTMP.Rd @@ -10,7 +10,8 @@ MSstatsSummarizeSingleTMP( censored_symbol, remove50missing, aft_iterations = 90, - aft_solver = "cholesky" + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ diff --git a/man/MSstatsSummarizeWithSingleCore.Rd b/man/MSstatsSummarizeWithSingleCore.Rd index 278ae8e1..94f4ceb3 100644 --- a/man/MSstatsSummarizeWithSingleCore.Rd +++ b/man/MSstatsSummarizeWithSingleCore.Rd @@ -12,7 +12,8 @@ MSstatsSummarizeWithSingleCore( remove50missing, equal_variance, aft_iterations = 90, - aft_solver = "cholesky" + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ diff --git a/man/dataProcess.Rd b/man/dataProcess.Rd index 5f306a4e..6859dbbf 100644 --- a/man/dataProcess.Rd +++ b/man/dataProcess.Rd @@ -26,7 +26,8 @@ dataProcess( log_file_path = NULL, numberOfCores = 1, aft_iterations = 90, - aft_solver = "cholesky" + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ @@ -127,7 +128,16 @@ track progress. Only works for Linux & Mac OS. Default is 1.} model's Newton-Raphson step: "cholesky" (default) delegates to \code{survival::survreg}, which solves it via Cholesky factorization. "cg" solves the same Newton step with a vendored conjugate-gradient -routine instead - an experimental alternative, currently opt-in only.} +routine instead; "pcg" is the same conjugate-gradient routine with a +Jacobi (inverse-diagonal) preconditioner, which can reduce the number +of conjugate-gradient iterations needed. "cg"/"pcg" are experimental +alternatives, currently opt-in only.} + +\item{aft_verbose}{If \code{TRUE} and \code{aft_solver} is "cg" or +"pcg", \code{message()} per-Newton-iteration conjugate-gradient +iteration counts and timing for every protein fit - useful for +evaluating solver time complexity, but produces one block of output +per protein, so leave at the default \code{FALSE} for routine runs.} } \value{ A list containing: diff --git a/man/dot-cgSolve.Rd b/man/dot-cgSolve.Rd index 5d846263..853d7e76 100644 --- a/man/dot-cgSolve.Rd +++ b/man/dot-cgSolve.Rd @@ -10,7 +10,8 @@ gradient} right_hand_side, initial_guess = NULL, relative_tolerance = 1e-08, - max_iterations = 10 * nrow(coefficient_matrix) + max_iterations = 10 * nrow(coefficient_matrix), + use_jacobi_preconditioner = FALSE ) } \arguments{ @@ -24,24 +25,42 @@ gradient/score vector from a Newton step.} to the zero vector.} \item{relative_tolerance}{how small the residual needs to shrink, -relative to the size of \code{right_hand_side}, before iteration stops.} +relative to the size of \code{right_hand_side}, before iteration stops. +Always judged on the true (unpreconditioned) residual, so this means the +same thing whether or not \code{use_jacobi_preconditioner} is set.} \item{max_iterations}{how many conjugate-gradient steps to try before giving up. In exact arithmetic, conjugate gradient converges within \code{nrow(coefficient_matrix)} steps, but rounding error erodes that guarantee as the system grows, so the default allows for several times that many steps.} + +\item{use_jacobi_preconditioner}{if \code{TRUE}, precondition with the +inverse of \code{coefficient_matrix}'s own diagonal - cheap to apply, +and often enough to cut down the number of iterations needed when the +diagonal dominates (as it typically does for an AFT information matrix, +where each parameter's own curvature tends to be much larger than its +cross-terms with the other parameters). Defaults to \code{FALSE}, which +reduces exactly to plain (unpreconditioned) conjugate gradient.} } \value{ -numeric vector solving (approximately) -\code{coefficient_matrix \%*\% solution = right_hand_side}. +a list with: \code{solution}, the numeric vector solving +(approximately) \code{coefficient_matrix \%*\% solution = +right_hand_side}; \code{iterations}, how many conjugate-gradient steps +were actually taken; \code{converged}, whether the residual tolerance +was met; and \code{positive_definite}, whether \code{coefficient_matrix} +behaved as positive definite throughout (a caller can fall back to a +different matrix, e.g. a Gauss-Newton approximation, when this is +\code{FALSE}). } \description{ A minimal, single right-hand-side conjugate gradient solver, used as the Newton-step linear solve in \code{.fitSurvivalCG}. Modeled on \code{lfe::cgsolve}, stripped down to a single dense matrix and a single right-hand-side vector (no multi-column batching, no \code{Matrix}-package -or operator/closure dispatch, no preconditioning - none of which are -needed for the small, dense AFT information matrices this is used on). +or operator/closure dispatch - neither is needed for the small, dense AFT +information matrices this is used on). Optionally applies a Jacobi +(inverse-diagonal) preconditioner, which \code{lfe::cgsolve} does not +support at all. } \keyword{internal} diff --git a/man/dot-fitAFTModel.Rd b/man/dot-fitAFTModel.Rd new file mode 100644 index 00000000..37f8f2b1 --- /dev/null +++ b/man/dot-fitAFTModel.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_imputation.R +\name{.fitAFTModel} +\alias{.fitAFTModel} +\title{Fit the AFT imputation model with the requested solver} +\usage{ +.fitAFTModel( + input, + aft_iterations, + aft_solver = "cholesky", + aft_verbose = FALSE +) +} +\arguments{ +\item{input}{data.table, the same shape \code{.fitSurvival} expects.} + +\item{aft_iterations}{maximum number of iterations for AFT model fitting.} + +\item{aft_solver}{"cholesky" (default, via \code{survival::survreg}), +"cg" (conjugate gradient), or "pcg" (conjugate gradient with a +Jacobi/inverse-diagonal preconditioner).} + +\item{aft_verbose}{passed through to \code{.fitSurvivalCG}'s +\code{verbose} when \code{aft_solver} is "cg" or "pcg"; has no effect +for "cholesky".} +} +\value{ +a fitted model of class \code{"survreg"}. +} +\description{ +Shared dispatch used by both \code{MSstatsSummarizeSingleLinear} and +\code{MSstatsSummarizeSingleTMP} so the \code{aft_solver}/ +\code{aft_verbose} logic lives in one place instead of being duplicated +at both call sites. +} +\keyword{internal} diff --git a/man/dot-fitSurvivalCG.Rd b/man/dot-fitSurvivalCG.Rd index 98ebeb0a..0aed5db9 100644 --- a/man/dot-fitSurvivalCG.Rd +++ b/man/dot-fitSurvivalCG.Rd @@ -5,7 +5,13 @@ \title{Fit a Gaussian, left-censored AFT model with a conjugate-gradient Newton step} \usage{ -.fitSurvivalCG(input, aft_iterations, convergence_tolerance = 1e-09) +.fitSurvivalCG( + input, + aft_iterations, + convergence_tolerance = 1e-09, + use_jacobi_preconditioner = FALSE, + verbose = FALSE +) } \arguments{ \item{input}{data.table, the same shape \code{.fitSurvival} expects.} @@ -15,9 +21,28 @@ Newton step} \item{convergence_tolerance}{stop once the change in log-likelihood between iterations falls below this (matches the default \code{rel.tolerance} in \code{survival::survreg.control}).} + +\item{use_jacobi_preconditioner}{if \code{TRUE}, precondition every +conjugate-gradient solve with the inverse of the current information +matrix's own diagonal (see \code{.cgSolve}'s +\code{use_jacobi_preconditioner}). This is what \code{aft_solver = +"pcg"} enables, versus plain conjugate gradient for \code{"cg"}.} + +\item{verbose}{if \code{TRUE}, \code{message()} a line per +Newton-Raphson iteration - conjugate-gradient iterations used, whether +the Gauss-Newton fallback (see below) was needed, elapsed time, and the +resulting log-likelihood - plus a one-line summary once fitting +finishes. Meant for evaluating how solver choice and problem size +trade off against iteration count and wall time, not for routine use +(this fits one protein at a time, so it is easy to generate a line per +protein across a whole \code{dataProcess()} run).} } \value{ -a fitted model of class \code{"survreg"}. +a fitted model of class \code{"survreg"}, with one added field: +\code{cg_diagnostics}, a data.frame with one row per Newton-Raphson +iteration recording the conjugate-gradient iteration counts and timing +described above (populated regardless of \code{verbose}, so it can be +inspected/aggregated programmatically after the fact). } \description{ An alternative to \code{.fitSurvival} for exactly the same imputation From 93aa0c6807f1d236b8516c7d44c28b8166a996c7 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Tue, 8 Sep 2026 21:26:51 -0400 Subject: [PATCH 3/8] fix imputation funneling of parameters into multicore --- R/MSstatsSummarizeWithMultipleCores.R | 15 ++++++++++----- man/MSstatsSummarizeWithMultipleCores.Rd | 2 ++ 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 1d4678fa..7af087cf 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -278,7 +278,7 @@ #' @noRd .build_summarize_worker <- function( use_TMP, impute, censored_symbol, remove50missing, - aft_iterations, equal_variance + aft_iterations, equal_variance, aft_solver, aft_verbose ) { unpack_fn <- .unpack_protein_slot use_TMP_ <- use_TMP @@ -287,6 +287,8 @@ remove50missing_ <- remove50missing aft_iterations_ <- aft_iterations equal_variance_ <- equal_variance + aft_solver_ <- aft_solver + aft_verbose_ <- aft_verbose function(record) { meta <- record$meta @@ -294,12 +296,13 @@ result <- if (use_TMP_) { MSstatsSummarizeSingleTMP( protein_dt, impute_, censored_symbol_, - remove50missing_, aft_iterations_) + remove50missing_, aft_iterations_, aft_solver_, aft_verbose_) } else { MSstatsSummarizeSingleLinear( protein_dt, impute_, censored_symbol_, remove50missing_, aft_iterations_, - equal_variances = equal_variance_) + equal_variances = equal_variance_, + aft_solver = aft_solver_, aft_verbose = aft_verbose_) } result } @@ -361,6 +364,8 @@ MSstatsSummarizeWithMultipleCores <- function( equal_variance, numberOfCores = 1L, aft_iterations = 90L, + aft_solver = "cholesky", + aft_verbose = FALSE, verbose = FALSE, BPPARAM = NULL, track_memory = FALSE, @@ -369,7 +374,7 @@ MSstatsSummarizeWithMultipleCores <- function( if (numberOfCores <= 1L && is.null(BPPARAM)) { return(MSstatsSummarizeWithSingleCore( input, method, impute, censored_symbol, - remove50missing, equal_variance, aft_iterations)) + remove50missing, equal_variance, aft_iterations, aft_solver, aft_verbose)) } start_time <- proc.time()[["elapsed"]] @@ -419,7 +424,7 @@ MSstatsSummarizeWithMultipleCores <- function( worker_fn <- .build_summarize_worker( use_TMP, impute, censored_symbol, remove50missing, - aft_iterations, equal_variance) + aft_iterations, equal_variance, aft_solver, aft_verbose) if (is.null(BPPARAM)) { tasks <- if (max_proteins_per_worker > 0L) { diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index b95b1e81..d2574fde 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -13,6 +13,8 @@ MSstatsSummarizeWithMultipleCores( equal_variance, numberOfCores = 1L, aft_iterations = 90L, + aft_solver = "cholesky", + aft_verbose = FALSE, verbose = FALSE, BPPARAM = NULL, track_memory = FALSE, From 0f177bfe0142f5b322ae9e3790fe1e43d3396f56 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Fri, 11 Sep 2026 16:12:01 -0400 Subject: [PATCH 4/8] add divergence in residual warning --- R/dataProcess.R | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/R/dataProcess.R b/R/dataProcess.R index 2a97689a..f2f3fe7b 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -485,6 +485,8 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, # Flag to track convergence warning converged = TRUE + convergence_messages = character(0) + diverging_warnings = 0L fit_data = if (is_labeled_reference) { single_protein[(!is_labeled_ref), cols, with = FALSE] @@ -496,12 +498,38 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, survival_fit = withCallingHandlers({ .fitAFTModel(fit_data, aft_iterations, aft_solver, aft_verbose) }, warning = function(w) { - if (grepl("converge", conditionMessage(w), ignore.case = TRUE)) { - message("Convergence warning caught: ", conditionMessage(w)) + warning_message = conditionMessage(w) + if (grepl("residual is diverging", warning_message, fixed = TRUE)) { + diverging_warnings <<- diverging_warnings + 1L + } + if (grepl("converge", warning_message, ignore.case = TRUE)) { + convergence_messages <<- c(convergence_messages, + warning_message) converged <<- FALSE } }) + protein_name = as.character(unique(single_protein$PROTEIN))[1] + log_fun = getOption("MSstatsLog") + if (diverging_warnings > 0L) { + msg = paste0("DIVERGING RESIDUAL for protein: ", protein_name, + " (", diverging_warnings, " warning(s))") + message(msg) + if (is.function(log_fun)) { + log_fun("INFO", msg) + } + } + if (!converged) { + msg = paste0("CONVERGENCE WARNING for protein: ", protein_name, + " (", length(convergence_messages), + " warning(s)) - ", + paste(unique(convergence_messages), collapse = " | ")) + message(msg) + if (is.function(log_fun)) { + log_fun("INFO", msg) + } + } + if (converged) { single_protein[, predicted := predict(survival_fit, newdata = .SD)] } else { From 85d93ec6ff62de87112fd736c6f34e8216d4b65f Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Fri, 11 Sep 2026 22:28:53 -0400 Subject: [PATCH 5/8] add verbose logging options for running survreg imputation --- R/dataProcess.R | 13 +++++++------ R/utils_imputation.R | 38 ++++++++++++++++++++++++++++++++------ man/reexports.Rd | 2 +- 3 files changed, 40 insertions(+), 13 deletions(-) diff --git a/R/dataProcess.R b/R/dataProcess.R index f2f3fe7b..8466c422 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -70,9 +70,10 @@ #' Jacobi (inverse-diagonal) preconditioner, which can reduce the number #' of conjugate-gradient iterations needed. "cg"/"pcg" are experimental #' alternatives, currently opt-in only. -#' @param aft_verbose If \code{TRUE} and \code{aft_solver} is "cg" or -#' "pcg", \code{message()} per-Newton-iteration conjugate-gradient -#' iteration counts and timing for every protein fit - useful for +#' @param aft_verbose If \code{TRUE}, \code{message()} diagnostics for +#' every protein fit: problem size and elapsed fitting time for all +#' solvers, plus per-Newton-iteration conjugate-gradient iteration counts +#' and timing when \code{aft_solver} is "cg" or "pcg" - useful for #' evaluating solver time complexity, but produces one block of output #' per protein, so leave at the default \code{FALSE} for routine runs. #' @inheritParams .documentFunction @@ -278,9 +279,9 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol #' model's Newton-Raphson step: "cholesky" (default, via #' \code{survival::survreg}), "cg" (conjugate gradient), or "pcg" #' (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). -#' @param aft_verbose If \code{TRUE} and \code{aft_solver} is "cg" or -#' "pcg", log per-Newton-iteration conjugate-gradient diagnostics for -#' every protein fit. See \code{.fitSurvivalCG}'s \code{verbose}. +#' @param aft_verbose If \code{TRUE}, log AFT fitting diagnostics for +#' every protein fit. See \code{.fitSurvival}'s and +#' \code{.fitSurvivalCG}'s \code{verbose}. #' #' @return list with protein-level data #' diff --git a/R/utils_imputation.R b/R/utils_imputation.R index 9aceda7d..824d49bf 100644 --- a/R/utils_imputation.R +++ b/R/utils_imputation.R @@ -54,13 +54,39 @@ } } +#' @param input data.table with the columns \code{.buildAFTFormula} needs. +#' @param aft_iterations maximum number of iterations for AFT model fitting. +#' @param verbose if \code{TRUE}, \code{message()} the problem size +#' (observations and parameters) before fitting and the wall time the fit +#' took afterwards, mirroring what \code{.fitSurvivalCG}'s \code{verbose} +#' reports. Meant for comparing solvers, not for routine use (this fits +#' one protein at a time). +#' +#' @importFrom stats model.frame model.matrix #' @importFrom survival survreg #' @keywords internal -.fitSurvival = function(input, aft_iterations) { +.fitSurvival = function(input, aft_iterations, verbose = FALSE) { # TODO: set.seed here? set.seed(100) - fit = survreg(.buildAFTFormula(input), data = input, dist = "gaussian", + aft_formula = .buildAFTFormula(input) + if (verbose) { + # survreg builds these internally; rebuilding them here is only + # worth the extra work when the counts are actually reported. + model_frame = model.frame(aft_formula, data = input) + design_matrix = model.matrix(attr(model_frame, "terms"), model_frame) + message(sprintf( + "[AFT-Cholesky] starting fit: %d observations, %d parameters", + nrow(design_matrix), ncol(design_matrix) + 1)) + } + fit_start_time = Sys.time() + fit = survreg(aft_formula, data = input, dist = "gaussian", control = list(maxiter = aft_iterations)) + if (verbose) { + message(sprintf( + "[AFT-Cholesky] finished: %d iterations, %.4f sec", + fit$iter[length(fit$iter)], + as.numeric(Sys.time() - fit_start_time, units = "secs"))) + } fit$y = NULL fit$linear.predictors = NULL fit @@ -544,9 +570,9 @@ #' @param aft_solver "cholesky" (default, via \code{survival::survreg}), #' "cg" (conjugate gradient), or "pcg" (conjugate gradient with a #' Jacobi/inverse-diagonal preconditioner). -#' @param aft_verbose passed through to \code{.fitSurvivalCG}'s -#' \code{verbose} when \code{aft_solver} is "cg" or "pcg"; has no effect -#' for "cholesky". +#' @param aft_verbose passed through to the chosen solver's +#' \code{verbose}: \code{.fitSurvivalCG}'s for "cg"/"pcg", +#' \code{.fitSurvival}'s for "cholesky". #' #' @return a fitted model of class \code{"survreg"}. #' @@ -559,7 +585,7 @@ } else if (aft_solver == "cg") { .fitSurvivalCG(input, aft_iterations, verbose = aft_verbose) } else { - .fitSurvival(input, aft_iterations) + .fitSurvival(input, aft_iterations, verbose = aft_verbose) } } 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 1194b969320121fbc4f4d185cd8df05c0df2eedb Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Sat, 12 Sep 2026 13:02:28 -0400 Subject: [PATCH 6/8] Get rid of divergence breakpoint in CG --- R/dataProcess.R | 12 ------------ R/utils_cgsolve.R | 14 -------------- 2 files changed, 26 deletions(-) diff --git a/R/dataProcess.R b/R/dataProcess.R index 8466c422..7fcd6766 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -487,7 +487,6 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, # Flag to track convergence warning converged = TRUE convergence_messages = character(0) - diverging_warnings = 0L fit_data = if (is_labeled_reference) { single_protein[(!is_labeled_ref), cols, with = FALSE] @@ -500,9 +499,6 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, .fitAFTModel(fit_data, aft_iterations, aft_solver, aft_verbose) }, warning = function(w) { warning_message = conditionMessage(w) - if (grepl("residual is diverging", warning_message, fixed = TRUE)) { - diverging_warnings <<- diverging_warnings + 1L - } if (grepl("converge", warning_message, ignore.case = TRUE)) { convergence_messages <<- c(convergence_messages, warning_message) @@ -512,14 +508,6 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, protein_name = as.character(unique(single_protein$PROTEIN))[1] log_fun = getOption("MSstatsLog") - if (diverging_warnings > 0L) { - msg = paste0("DIVERGING RESIDUAL for protein: ", protein_name, - " (", diverging_warnings, " warning(s))") - message(msg) - if (is.function(log_fun)) { - log_fun("INFO", msg) - } - } if (!converged) { msg = paste0("CONVERGENCE WARNING for protein: ", protein_name, " (", length(convergence_messages), diff --git a/R/utils_cgsolve.R b/R/utils_cgsolve.R index 4384245c..6cd57efe 100644 --- a/R/utils_cgsolve.R +++ b/R/utils_cgsolve.R @@ -73,7 +73,6 @@ residual_size = sum(residual * residual) residual_dot_preconditioned_residual = sum(residual * preconditioned_residual) - smallest_residual_size_seen = residual_size # Stop once the residual has shrunk far enough, relative to the size of # the right-hand side (falling back to an absolute scale when that size @@ -109,19 +108,6 @@ solution = solution + step_length * search_direction residual = residual - step_length * matrix_times_search_direction new_residual_size = sum(residual * residual) - smallest_residual_size_seen = - min(smallest_residual_size_seen, new_residual_size) - - # If the residual has grown far past its best value so far, the - # iteration is diverging (e.g. because coefficient_matrix is - # ill-conditioned) - give up and return what we have rather than - # loop until max_iterations. - if (iteration > 10 && - new_residual_size > 1e4 * smallest_residual_size_seen) { - warning(".cgSolve: residual is diverging; returning the best ", - "iterate found so far") - break - } # Choose the next search direction so it doesn't undo the progress # made by earlier directions. From 0566a9fbb78af8930785f9c0bba6d2f3bd6c8731 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Sun, 13 Sep 2026 17:39:31 -0400 Subject: [PATCH 7/8] Remove the artificial number_of_halvings > 30 hard-abort. Instead, let halving consume the shared iteration budget (aft_iterations), and when that budget runs out, fall back to the last accepted coefficients/log_scale/current_fit rather than discarding everything and reporting a hard failure. --- R/utils_imputation.R | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/R/utils_imputation.R b/R/utils_imputation.R index 824d49bf..8570facc 100644 --- a/R/utils_imputation.R +++ b/R/utils_imputation.R @@ -279,7 +279,11 @@ #' result is used. #' #' @param input data.table, the same shape \code{.fitSurvival} expects. -#' @param aft_iterations maximum number of Newton-Raphson iterations. +#' @param aft_iterations maximum number of log-likelihood evaluations the +#' fit may spend. Newton-Raphson iterations and the step-halvings used to +#' recover from an overshooting step share this one budget; once it is +#' exhausted fitting stops and the last accepted coefficients and scale +#' are returned (with a non-convergence warning), rather than failing. #' @param convergence_tolerance stop once the change in log-likelihood #' between iterations falls below this (matches the default #' \code{rel.tolerance} in \code{survival::survreg.control}). @@ -431,9 +435,13 @@ current_log_likelihood = current_fit$log_likelihood number_of_iterations_used = 0 converged = FALSE + iterations_remaining = aft_iterations cg_diagnostics = vector("list", aft_iterations) - for (iteration in seq_len(aft_iterations)) { + iteration = 0 + while (iterations_remaining > 0) { + iteration = iteration + 1 + iterations_remaining = iterations_remaining - 1 number_of_iterations_used = iteration iteration_start_time = Sys.time() @@ -470,20 +478,16 @@ # strategy (survreg6.c) rather than simply rejecting the step # outright. number_of_halvings = 0 - halving_exhausted = FALSE repeat { candidate_fit = evaluate_log_likelihood_and_derivatives( candidate_coefficients, candidate_log_scale) candidate_improves = is_finite_fit(candidate_fit) && candidate_fit$log_likelihood >= current_log_likelihood - if (candidate_improves) { + if (candidate_improves || iterations_remaining <= 0) { break } + iterations_remaining = iterations_remaining - 1 number_of_halvings = number_of_halvings + 1 - if (number_of_halvings > 30) { - halving_exhausted = TRUE - break - } if (number_of_halvings == 1 && (log_scale - candidate_log_scale) > 1.1) { # a single huge drop in scale is the most common cause of @@ -496,7 +500,7 @@ candidate_log_scale = (candidate_log_scale + 2 * log_scale) / 3 } - if (halving_exhausted) { + if (!candidate_improves) { break } @@ -518,8 +522,9 @@ } if (!converged) { - warning("AFT model (CG solver) ran out of iterations and did not ", - "converge") + warning("AFT model (CG solver) used its full iteration budget ", + "without converging; returning the last accepted ", + "coefficients") } cg_diagnostics = do.call( From 8b180ab53fc82c8a59ad5d3b5b608e7f96b90993 Mon Sep 17 00:00:00 2001 From: tonywu1999 Date: Thu, 17 Sep 2026 15:10:39 -0400 Subject: [PATCH 8/8] set number of threads to 1 per core for blas operations for pcg --- DESCRIPTION | 5 +- NAMESPACE | 195 ++++++++++++++------------ R/MSstatsSummarizeWithMultipleCores.R | 2 + man/MSstatsSummarizeSingleLinear.Rd | 6 +- man/dataProcess.Rd | 7 +- man/dot-fitAFTModel.Rd | 6 +- man/dot-fitSurvivalCG.Rd | 6 +- man/reexports.Rd | 2 +- 8 files changed, 130 insertions(+), 99 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d890bee3..818a948b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -40,7 +40,8 @@ Imports: parallel, rlang, matter, - BiocParallel + BiocParallel, + RhpcBLASctl Suggests: BiocStyle, knitr, @@ -62,4 +63,4 @@ Packaged: 2017-10-20 02:13:12 UTC; meenachoi LinkingTo: Rcpp, RcppArmadillo -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index 4bb4e353..5f156811 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -55,101 +55,124 @@ 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(BiocParallel, + bpisup, + bplapply, + bpnworkers, + bpprogressbar, + bpstart, + bpstop +) importFrom(MASS,rlm) -importFrom(MSstatsConvert,DIANNtoMSstatsFormat) -importFrom(MSstatsConvert,DIAUmpiretoMSstatsFormat) -importFrom(MSstatsConvert,FragPipetoMSstatsFormat) -importFrom(MSstatsConvert,MSstatsBalancedDesign) -importFrom(MSstatsConvert,MSstatsClean) -importFrom(MSstatsConvert,MSstatsImport) -importFrom(MSstatsConvert,MSstatsLogsSettings) -importFrom(MSstatsConvert,MSstatsMakeAnnotation) -importFrom(MSstatsConvert,MSstatsPreprocess) -importFrom(MSstatsConvert,MZMinetoMSstatsFormat) -importFrom(MSstatsConvert,MaxQtoMSstatsFormat) -importFrom(MSstatsConvert,OpenMStoMSstatsFormat) -importFrom(MSstatsConvert,OpenSWATHtoMSstatsFormat) -importFrom(MSstatsConvert,PDtoMSstatsFormat) -importFrom(MSstatsConvert,ProgenesistoMSstatsFormat) -importFrom(MSstatsConvert,SkylinetoMSstatsFormat) -importFrom(MSstatsConvert,SpectronauttoMSstatsFormat) +importFrom(MSstatsConvert, + DIANNtoMSstatsFormat, + DIAUmpiretoMSstatsFormat, + FragPipetoMSstatsFormat, + MSstatsBalancedDesign, + MSstatsClean, + MSstatsImport, + MSstatsLogsSettings, + MSstatsMakeAnnotation, + MSstatsPreprocess, + MZMinetoMSstatsFormat, + MaxQtoMSstatsFormat, + OpenMStoMSstatsFormat, + OpenSWATHtoMSstatsFormat, + PDtoMSstatsFormat, + ProgenesistoMSstatsFormat, + SkylinetoMSstatsFormat, + 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(RhpcBLASctl,blas_set_num_threads) +importFrom(data.table, + as.data.table, + data.table, + fifelse, + melt, + rbindlist, + setDT, + setDTthreads, + uniqueN +) importFrom(ggrepel,geom_text_repel) importFrom(gplots,heatmap.2) -importFrom(grDevices,dev.off) -importFrom(grDevices,hcl) -importFrom(grDevices,pdf) -importFrom(graphics,axis) -importFrom(graphics,image) -importFrom(graphics,legend) -importFrom(graphics,mtext) -importFrom(graphics,par) -importFrom(graphics,plot) -importFrom(graphics,plot.new) -importFrom(graphics,title) -importFrom(htmltools,div) -importFrom(htmltools,save_html) -importFrom(htmltools,tagList) +importFrom(grDevices, + dev.off, + hcl, + pdf +) +importFrom(graphics, + axis, + image, + legend, + mtext, + par, + plot, + plot.new, + title +) +importFrom(htmltools, + div, + save_html, + tagList +) importFrom(limma,squeezeVar) importFrom(lme4,lmer) importFrom(marray,maPalette) importFrom(matter,SnowfastParam) importFrom(methods,is) -importFrom(parallel,clusterExport) -importFrom(parallel,makeCluster) -importFrom(parallel,parLapply) -importFrom(parallel,stopCluster) -importFrom(plotly,add_trace) -importFrom(plotly,ggplotly) -importFrom(plotly,layout) -importFrom(plotly,plot_ly) -importFrom(plotly,style) -importFrom(plotly,subplot) +importFrom(parallel, + clusterExport, + makeCluster, + parLapply, + stopCluster +) +importFrom(plotly, + add_trace, + ggplotly, + layout, + plot_ly, + style, + subplot +) importFrom(preprocessCore,normalize.quantiles) importFrom(rlang,.data) -importFrom(stats,dist) -importFrom(stats,dnorm) -importFrom(stats,fitted) -importFrom(stats,formula) -importFrom(stats,hclust) -importFrom(stats,lm) -importFrom(stats,lm.fit) -importFrom(stats,loess) -importFrom(stats,median) -importFrom(stats,model.frame) -importFrom(stats,model.matrix) -importFrom(stats,model.response) -importFrom(stats,na.omit) -importFrom(stats,p.adjust) -importFrom(stats,pnorm) -importFrom(stats,predict) -importFrom(stats,qbinom) -importFrom(stats,qnorm) -importFrom(stats,qt) -importFrom(stats,quantile) -importFrom(stats,resid) -importFrom(stats,residuals) -importFrom(stats,sd) -importFrom(stats,vcov) -importFrom(stats,xtabs) -importFrom(survival,Surv) -importFrom(survival,survreg) -importFrom(utils,combn) -importFrom(utils,sessionInfo) -importFrom(utils,setTxtProgressBar) -importFrom(utils,txtProgressBar) +importFrom(stats, + dist, + dnorm, + fitted, + formula, + hclust, + lm, + lm.fit, + loess, + median, + model.frame, + model.matrix, + model.response, + na.omit, + p.adjust, + pnorm, + predict, + qbinom, + qnorm, + qt, + quantile, + resid, + residuals, + sd, + vcov, + xtabs +) +importFrom(survival, + Surv, + survreg +) +importFrom(utils, + combn, + sessionInfo, + setTxtProgressBar, + txtProgressBar +) useDynLib(MSstats, .registration=TRUE) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 7af087cf..fff63a15 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -316,9 +316,11 @@ list(worker = i, pid = Sys.getpid(), max_rss_mb = .max_rss_mb()) } +#' @importFrom RhpcBLASctl blas_set_num_threads .warmup_worker <- function(i) { library(MSstats, quietly = TRUE, warn.conflicts = FALSE) data.table::setDTthreads(1) + RhpcBLASctl::blas_set_num_threads(1) NULL } diff --git a/man/MSstatsSummarizeSingleLinear.Rd b/man/MSstatsSummarizeSingleLinear.Rd index dc6d006e..ed60f8c4 100644 --- a/man/MSstatsSummarizeSingleLinear.Rd +++ b/man/MSstatsSummarizeSingleLinear.Rd @@ -33,9 +33,9 @@ model's Newton-Raphson step: "cholesky" (default, via \code{survival::survreg}), "cg" (conjugate gradient), or "pcg" (conjugate gradient with a Jacobi/inverse-diagonal preconditioner).} -\item{aft_verbose}{If \code{TRUE} and \code{aft_solver} is "cg" or -"pcg", log per-Newton-iteration conjugate-gradient diagnostics for -every protein fit. See \code{.fitSurvivalCG}'s \code{verbose}.} +\item{aft_verbose}{If \code{TRUE}, log AFT fitting diagnostics for +every protein fit. See \code{.fitSurvival}'s and +\code{.fitSurvivalCG}'s \code{verbose}.} } \value{ list with protein-level data diff --git a/man/dataProcess.Rd b/man/dataProcess.Rd index 6859dbbf..99c6497d 100644 --- a/man/dataProcess.Rd +++ b/man/dataProcess.Rd @@ -133,9 +133,10 @@ Jacobi (inverse-diagonal) preconditioner, which can reduce the number of conjugate-gradient iterations needed. "cg"/"pcg" are experimental alternatives, currently opt-in only.} -\item{aft_verbose}{If \code{TRUE} and \code{aft_solver} is "cg" or -"pcg", \code{message()} per-Newton-iteration conjugate-gradient -iteration counts and timing for every protein fit - useful for +\item{aft_verbose}{If \code{TRUE}, \code{message()} diagnostics for +every protein fit: problem size and elapsed fitting time for all +solvers, plus per-Newton-iteration conjugate-gradient iteration counts +and timing when \code{aft_solver} is "cg" or "pcg" - useful for evaluating solver time complexity, but produces one block of output per protein, so leave at the default \code{FALSE} for routine runs.} } diff --git a/man/dot-fitAFTModel.Rd b/man/dot-fitAFTModel.Rd index 37f8f2b1..a66883e2 100644 --- a/man/dot-fitAFTModel.Rd +++ b/man/dot-fitAFTModel.Rd @@ -20,9 +20,9 @@ "cg" (conjugate gradient), or "pcg" (conjugate gradient with a Jacobi/inverse-diagonal preconditioner).} -\item{aft_verbose}{passed through to \code{.fitSurvivalCG}'s -\code{verbose} when \code{aft_solver} is "cg" or "pcg"; has no effect -for "cholesky".} +\item{aft_verbose}{passed through to the chosen solver's +\code{verbose}: \code{.fitSurvivalCG}'s for "cg"/"pcg", +\code{.fitSurvival}'s for "cholesky".} } \value{ a fitted model of class \code{"survreg"}. diff --git a/man/dot-fitSurvivalCG.Rd b/man/dot-fitSurvivalCG.Rd index 0aed5db9..38ab07a4 100644 --- a/man/dot-fitSurvivalCG.Rd +++ b/man/dot-fitSurvivalCG.Rd @@ -16,7 +16,11 @@ Newton step} \arguments{ \item{input}{data.table, the same shape \code{.fitSurvival} expects.} -\item{aft_iterations}{maximum number of Newton-Raphson iterations.} +\item{aft_iterations}{maximum number of log-likelihood evaluations the +fit may spend. Newton-Raphson iterations and the step-halvings used to +recover from an overshooting step share this one budget; once it is +exhausted fitting stops and the last accepted coefficients and scale +are returned (with a non-convergence warning), rather than failing.} \item{convergence_tolerance}{stop once the change in log-likelihood between iterations falls below this (matches the default 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()}}} }}