Baseline psalsa improvements - #118
Open
zeehio wants to merge 27 commits into
Open
Conversation
Plots a few samples' original signal together with their estimated baseline (from nmr_baseline_estimation()), so the baseline estimate can be visually inspected against the signal. - x = chemshift, y = intensity, colour = NMRExperiment - linetype: solid for the signal, dashed for the baseline - chemshift_range accepts a single length-2 vector, or a named list of them to facet several regions side by side (facet_grid with NMRExperiment in rows and, when more than one region is requested, regions in columns) - validates the requested region(s) actually overlap the dataset's ppm axis, with a clear error instead of a downstream failure
facet_grid() has no built-in pagination, unlike facet_wrap(). Since this function grids NMRExperiment (rows) against chemshift_range regions (columns), nrow/ncol paginate each dimension independently: nrow caps NMRExperiments per page, ncol caps regions per page. page indexes through the resulting grid of pages, NMRExperiment pages varying fastest. Defaults to showing everything on one page, matching the previous behavior.
Drop the NMRExperiment facet dimension: every requested NMRExperiment is now overlaid (by colour) within each region's facet_wrap() panel, instead of getting its own facet_grid() row. nrow/ncol/page now paginate the region facets directly, reusing the same snug-grid default (1xn / 2x2 / 2x3 / 3x3) as nmr_baseline_threshold_plot().
Regions can have very different intensity ranges; sharing one y-axis squashes smaller-signal regions flat next to larger ones. Free y (in addition to the already-free x) makes each region's own signal/ baseline shape legible.
NMR peaks are Lorentzian-shaped, unlike gen_synthetic_1d()'s existing "gaussian" and "gex" options. lorentz_peak_1d() is a height/fwhm- parametrized wrapper around the package's existing canonical lorentzian() (area_estimation.R), the same shape peaklist_fit_lorentzians() already fits, so the synthetic tuning pool can now match what NMR peak detection actually assumes. Default peak_shape stays "gaussian" (first choice), so this is purely additive for existing callers of tune_psalsa()/nmr_baseline_estimation().
AlpsNMR's synthetic tuning pool now matches NMR peaks' actual shape by default. nmr_baseline_estimation() calls tune_psalsa() without specifying peak_shape, so its automatic lambda/p/k tuning now also defaults to lorentzian-shaped synthetic peaks instead of gaussian.
peak_area_errors_1d() used to sum the baseline-corrected signal over
each peak's own lo:hi window and compare to its true area -- valid
only when peaks don't overlap, since a neighbor's tail leaking into
that window contaminates the estimate. compute_peak_size_breaks_1d()
and tune_psalsa_params_1d()'s use_peak_area gate then had to restrict
to isolated peaks only (peak_is_isolated_1d()) to stay valid.
For Lorentzian peaks (now the default), that's nearly all of them:
measured 0% isolated on this project's real NMR data at realistic
peak density (vs. 11.2% for Gaussian, same placement), because a
Lorentzian's 1/x^2 tail crosses the 0.1%-height isolation threshold
~10x further from center than a Gaussian's at the same FWHM (measured
31.6x vs 3.1x FWHM). use_peak_area was always false as a result, so
tuning silently fell back to the cruder floor_rmse_1d objective for
any realistically dense synthetic pool -- exactly the regime real
crowded NMR spectra live in, so parameters tuned this way may not
transfer well to them.
Fix: since these are synthetic peaks, we know each one's own true
contribution at every point even where they overlap (gen_synthetic_1d()
now keeps enough per-peak info -- including gex's per-peak a/b, and the
peak_shape itself -- to reconstruct any single peak's curve on demand via
the new reconstruct_peak_1d()). peak_area_errors_1d() now splits the
corrected signal at each point in proportion to each peak's true relative
contribution there, instead of a hard window:
w_i(x) = pk_i(x) / peaks_total(x)
est_area_i = sum_x w_i(x) * corrected(x)
This reduces exactly to the old sum(corrected[lo:hi]) for a genuinely
isolated peak (verified: exact match), and gives every peak a valid,
uncontaminated area estimate regardless of overlap (verified: two
heavily-overlapping synthetic peaks under perfect baseline correction
recover their true areas independently, both at the same small
truncation-only error, summing to ~99.95% of the corrected signal in
their shared window). compute_peak_size_breaks_1d() and the
use_peak_area gate (now n_total >= 15, not n_iso >= 15) use the whole
peak population accordingly. Confirmed use_peak_area now correctly
flips true for lorentzian on this project's real data (was false).
This prioritizes correctness over speed: peak_area_errors_1d() now
reconstructs each peak's own curve inside the objective function,
adding cost on top of what profiling already found to be the dominant
step (phase 3, the Nelder-Mead search). Revisit if that cost matters
in practice.
reconstruct_peak_1d() and peak_area_errors_1d() are both @nord, so [reconstruct_peak_1d()] and [peak_area_errors_1d()] never resolved to anything (roxygen2 warned on every regeneration). Replaced with plain text.
…sa.R Both are @nord (diff2_penalty undocumented) or nonexistent (psalsa2d) so roxygen2 could not resolve the [] links; replaced with plain code text.
With no faceting to page through, pagination instead subsets which NMRExperiments (in the existing numeric-sorted order) land on the sample axis, samples_per_page at a time; requesting a page beyond the number available errors with the page count.
$data1r doesn't exist (the field is $data_1r), so mat silently became NULL, then NULL-minus-baseline lost its matrix dims via recycling. mat is already restricted to threshold_ind columns before the median3mad loop, so re-indexing by threshold_ind there was a second, redundant subscript that errored once mat regained the right dims. Both methods now verified end-to-end.
…plitting On the full MTBLS242 dataset (134 samples, 47417 peaks), nmr_peak_clustering() tried to allocate a single 47417x47417 dense distance matrix (~18GB) and OOM'd. A valid complete-linkage cluster's ppm span must stay under max_dist_thresh_ppb (how num_clusters is chosen), so in 1D, any two peaks sharing a final cluster are necessarily within that same threshold of each other directly. Splitting the peak list at every ppm gap exceeding that threshold and clustering each resulting group independently (split_peaks_into_ppm_components() + cluster_one_ppm_component()) therefore gives identical results to clustering everything at once, while keeping each group's distance matrix down to the size of one crowded region instead of the whole spectrum. On the real dataset this drops peak memory from ~18GB (OOM) to ~4-9GB and completes in about a minute; verified to reproduce the same partition as the old single-matrix path on both a small worked example and a synthetic multi-region case. Also: - set_peak_distances_within_groups() no longer round-trips through a dense matrix to mask same-sample distances to Inf; it assigns directly into the triangular dist vector via index arithmetic, roughly halving peak memory for whichever group ends up largest. - Fixed a latent bug in get_max_dist_ppb_for_num_clusters(): cutree() returns a plain vector instead of a matrix when only one candidate k is tested (which happens whenever a region's peaks all come from a single sample), and the subsequent matrix-style indexing broke on that. Pre-existing, but only surfaced once region splitting made single-sample regions common. nmr_peak_clustering(peak2peak_dist=...) or nmr_peak_clustering(num_clusters=...) still take the original single global hclust path unchanged, since neither has a distance threshold to safely split by.
accept_inflections = FALSE rejects peaks whose apex isn't a genuine local maximum of the raw signal within its own fitted window (ppm_infl_min - ppm_infl_max) -- typically shoulder/inflection points detected on the flank of a much larger neighbouring peak in crowded regions, rather than peaks of their own. Defaults to TRUE (previous behaviour unchanged). is_peak_local_max()/peak_is_local_max_1() check that the peak's own pos is (within a 1% relative tolerance, to absorb fitting noise on a flat/rounded apex) the window's maximum, and that intensity moves monotonically toward it from both edges. The tolerance is scaled by abs(window max), not the raw max, since a window can sit entirely below zero in flat/noisy baseline regions -- an earlier version without the abs() inverted the tolerance direction there, which a regression test now covers. Verified against the real MTBLS242 dataset (134 samples, 70178 peaks): accept_inflections = FALSE rejects 8829 shoulder peaks (12.6%), matching a prior manual analysis of the same data to within a handful of peaks (the remainder explained by peaklist_accept_peaks()'s pre-existing excluded-region filter, which that manual analysis didn't apply). Runs in ~4.3s on the full dataset.
A fitted peak's area is height * pi * gamma, so it scales with linewidth as much as with height -- area_min disproportionately rejects genuinely sharp, narrow peaks (formate being the clearest example: gamma_ppb ~1-3 vs ~2-45 for a typical crowded-region peak, so its area comes out tiny despite a tall, unambiguous peak). intensity_min/intensity_max filter on peak height directly instead, so they don't carry that width bias. Verified against the real MTBLS242 dataset: intensity_min = 500 keeps 167/175 (95.4%) of formate-apex detections, versus 3/175 (1.7%) with the old area_min = 50 approach on the same peaks. Added 3 new unit tests (including a no-op-by-default check) and fixed 3 existing tests whose fake peak_data fixtures lacked an intensity column, which the new criteria now require. All 752 package tests pass.
Splits the input signal(s) into num_regions contiguous stretches and characterizes peak density/fwhm separately per region (instead of one pooled, signal-wide density/fwhm), so the synthetic tuning pool can reproduce a spectrum's real region-to-region variation (e.g. a crowded aliphatic region vs. a sparse downfield region) rather than blending it into one uniform density. Opt-in via tune_psalsa(..., num_regions = N); default (NULL) keeps the original signal-wide pooling untouched. lambda/p/k are still tuned as single, signal-wide values either way -- this only changes how realistic the synthetic pool is.
Generalizes psalsa() to accept lambda/p/k as position-varying vectors (one value per point) instead of only scalars, and adds tune_psalsa_spatial() to build such profiles automatically: each of num_regions regions is tuned independently against its own region-matched synthetic pool, and the per-region optima are combined into smooth lambda/p/k profiles via a natural cubic spline (in log-space for lambda/k, logit-space for p, so the inverse transform always lands back in a valid range regardless of spline over/undershoot between knots). diff2_penalty_weighted() builds the position-varying smoothing penalty directly as t(D2) %*% diag(lambda) %*% D2 via a sqrt-then-crossprod construction, preserving D2's pentadiagonal sparsity pattern exactly -- no memory blowup from going position-varying. psalsa_core()'s p/k are now broadcast to the signal's own length before the reweighting step, fixing a latent misalignment that a vector p/k would otherwise hit when subsetting by the residual sign mask. The existing scalar-lambda/p/k path is untouched (byte-identical), tune_psalsa()'s default behaviour is unaffected, and tune_psalsa_spatial() falls back to tune_psalsa() when fewer than 2 regions have enough signal to spline between. Verified on the two MTBLS242 samples with the known baseline-bleed issue: this cuts the median baseline-corrected offset in the affected downfield window from ~362-386 down to ~85-90, versus only ~7-8% improvement from region-aware synthetic data alone (tune_psalsa's num_regions), confirming the bleed needed a position-varying penalty rather than just more realistic tuning data.
The theta0 anchor added previously wasn't enough on its own: a region with very few peaks in its own synthetic pool still collapsed to the same "disable peak protection" degenerate optimum (lambda -> tiny, k -> huge, p -> near p_max) already described in tune_psalsa_params_1d()'s own comments, at the SAME regularization weight used for the whole-spectrum search. Reproduced directly on a real MTBLS242 region (20/20, 3 peaks/draw): lambda collapsed to ~400 vs. an anchor of ~2.8e7 even with the anchor in place, and the resulting position- varying baseline interpolated the raw signal almost exactly in that region for every sample, not just the ones the tuning was meant to fix. Adds region_lambda_k_weight/region_p_weight (default 2/10, ~100x/~25x the whole-spectrum defaults) so each region's own search is pulled much more strongly toward the signal-wide anchor -- verified empirically across independent synthetic draws to keep a sparse region's tuned values close to the anchor instead of drifting to the escape hatch. Also corrects course on the earlier real-data validation: the ad hoc comparison window used for the earlier commits (8.20-8.60 ppm) was NOT the vignette's actual nmr_baseline_threshold() reference region (9.5-10 ppm). Redone against the correct window: the region-aware synthetic data alone (tune_psalsa(num_regions = 20), already shipped) turns out to already bring the two flagged samples' threshold from 18847/12447 down to ~921/835 -- essentially into the normal range (420-1335 across other samples) by itself. With this regularization fix, tune_psalsa_spatial() now lands at ~893/819 (matching, not dramatically beating, the simpler region-aware tune) with all 18 tuned regions staying sanely clustered (previously as many as 4-5 orders of magnitude apart).
Addresses the sparse-region collapse at its root cause instead of only fighting it with a stronger ridge: merge_sparse_regions_1d() greedily merges ADJACENT regions (left-to-right, backward-merging a trailing under-informed group into its predecessor -- but never force-merging a genuinely empty trailing region into an already well-populated neighbour, which would only dilute it) until each group has an estimated peak count of at least min_peaks (default 15, matching tune_psalsa_params_1d()'s own peak-tiering threshold). Wired into tune_psalsa_region_params_1d()/tune_psalsa_spatial() via a new min_peaks parameter. This is complementary to, not a replacement for, the region_lambda_k_weight/region_p_weight regularization added previously: merging gives a region real data to constrain its own search, while the regularization stays as a safety net for whatever a region still can't resolve on its own (matching tune_psalsa_params_1d()'s own rationale for why it regularizes even a well-populated whole-spectrum search). Verified on the real MTBLS242 data: with the default num_regions = 20, merging collapses the 20 raw bins down to 12 well-informed groups (the previously-pathological sparse region gets absorbed into a wider, properly-populated neighbour), and every tuned lambda now lands in a sane 7.5e6-1.6e8 range (previously spanning 4-5 orders of magnitude). The two originally-flagged samples' reference-window threshold (correct 9.5-10 ppm window) now lands at 777/713, comfortably inside the normal range (764-890) seen across other samples.
Adds a num_regions argument (default NULL, unchanged behaviour): when set, an "auto" lambda/p/k is tuned with tune_psalsa_spatial(num_regions = num_regions) instead of the plain, signal-wide tune_psalsa(), giving a position-varying baseline profile through the actual exported entry point rather than only via internal AlpsNMR::: calls. psalsa()'s own vector-lambda/p/k support (added earlier) means no change was needed to how the tuned values get applied -- the existing do.call(psalsa, ...) call already handles a vector or scalar lambda/p/k transparently. The psalsa_params attribute on data_1r_baseline now also records num_regions alongside lambda/p/k, so a caller can tell whether a stored baseline came from a scalar or position-varying tune; lambda/p/k stay directly reusable via nmr_baseline_estimation(other_dataset, lambda = psalsa_params$lambda, p = psalsa_params$p, k = psalsa_params$k) without retuning, same as before.
…file tune_psalsa_spatial()'s lambda/p/k profiles are as long as the tuned spectrum, impractical to store verbatim (e.g. as a literal for reuse without re-tuning); the much smaller region_params table (a few dozen rows at most) plus noise_sd already determine them exactly via the same spline interpolation tune_psalsa_spatial() itself uses. Adds psalsa_region_params_to_profile(n, region_params, noise_sd) to rebuild the full profiles from just those two small pieces -- verified to reproduce a tune_psalsa_spatial() result's lambda/p/k exactly. Also adds noise_sd to tune_psalsa_spatial()'s own return value (needed to do this reconstruction at all; previously computed internally but never returned). Left unexported, matching the rest of the psalsa/tune_psalsa/ tune_psalsa_spatial family: this stays internal (AlpsNMR:::) machinery for now rather than public API.
Fixes a real, previously-undiagnosed convergence bug: the reweighting iteration depends on a hard threshold (d_geq <- d >= 0) deciding whether each point sits above or below the current baseline estimate. Near a flat/quiet region, many residuals sit close to zero, so a full undamped step in the fitted curve can flip a large batch of points across that threshold at once, which swings the next fitted curve enough to flip a comparable batch back -- an oscillation that can persist right up to maxit instead of settling, rather than slow but monotonic convergence. Verified directly on the MTBLS242 dataset: instrumented traces showed the flip count dropping toward 0 (near-converged) then suddenly jumping back into the thousands at iteration 23-24, for both the already-tuned scalar parameters AND a normal control sample -- this was not specific to any one sample or tuning approach. damping < 1 under-relaxes the weight update (w <- damping * w_target + (1 - damping) * w) instead of replacing w outright, directly shrinking how much the fitted curve moves per iteration and therefore how many points can cross the threshold in one step. damping = 1 (the default) is an exact, byte-identical replacement of the original update -- zero behavioural change unless a caller opts in. Verified on the full 134-sample dataset with damping = 0.7, maxit = 100 (using the SAME already-tuned scalar lambda/p/k, no retuning): every sample now converges within 100 iterations (min 38, median 54.5, max 75), and the reference-window threshold that was previously an 18847 outlier for one sample drops to 1056, with every other sample's threshold landing in a tight, sane 1000-2200 range (previously up to 18847) -- eliminating the anomaly that motivated this session's region-based tuning investigation, without needing any of that region-based machinery.
Covers: damping = 1 is byte-identical to omitting the argument (both psalsa_core() directly and psalsa()/psalsa_one()); the damped update is a proper convex combination of the target and previous weight; damping = 0 keeps every smoother() call on the original unweighted w; damping changes psalsa()'s output relative to the default; damping is threaded through the matrix (multi-sample) path; and a flat, borderline-noisy region (many residuals near zero, structurally the scenario that can trigger the hard-threshold oscillation this fix targets) converges reliably with damping = 0.5.
…uning gen_synthetic_1d() draws peak heights from a lognormal centered on 0.35*A, but every per-region synthetic pool (both tune_psalsa()'s num_regions path via gen_synthetic_1d_regions(), and tune_psalsa_spatial()'s independent per-region tuning via tune_psalsa_region_params_1d()) was calling it with A left at its default of 1 for every region, regardless of that region's own real peak height scale. Only density and fwhm varied by region; amplitude never did, even though it's arguably the dimension that matters most for k (the peak-height parameter) specifically. Verified on real MTBLS242 data this was a real gap, not cosmetic: the BCAA region's real peaks are ~10,000-700,000 in height, four orders of magnitude above the ~100-700 height of small real features in the "quiet" reference region -- yet both were being synthesized from the exact same unit-amplitude peak-height distribution. Adds a per-region `amplitude` (median detected-peak height / 0.35, same convention as analyze_signal_1d()'s own A_peaks) to analyze_signal_regions_1d()/pool_signal_stats_regions_1d(), propagated through merge_sparse_regions_1d() (peak-count-weighted mean across merged members) and into gen_synthetic_1d_regions()/ tune_psalsa_region_params_1d(). csnr is correspondingly rescaled per region (csnr_region = noise_sd / amplitude) so the ABSOLUTE noise level stays equal to the signal-wide noise_sd regardless of a region's own amplitude -- real (thermal/electronic) noise doesn't scale with local peak height the way peak heights themselves should. Layered, fully backward-compatible fallback when amplitude can't be determined: falls back to the signal-wide amplitude (noise_sd/csnr) when available, and to the ORIGINAL pre-amplitude behaviour (no A override, region_profile's own csnr as-is) when neither amplitude nor noise_sd is available -- verified byte-identical to the old behaviour in that case.
…rough nmr_baseline_estimation() compute_psalsa_weights() replaces the hard d>=0 threshold with a cubic-Hermite-interpolated central noise band [-s, s], avoiding the p-floor bias the threshold method imposes on every above-baseline point regardless of how close it is to the baseline. Guards against c_noise > p, which makes the left-region exponential diverge instead of saturating to w_max. weight_method = "smooth" only supports scalar p/k (not tune_psalsa_spatial()'s position-varying profiles). nmr_baseline_estimation() gains a damping = "auto" parameter mirroring its existing maxit = "auto" pattern. Validated against the full 134-sample MTBLS242 dataset: smooth weights converge for every sample with a finite baseline, and the reference-window (9.5-10 ppm) threshold drops from a median of 1299 (damped threshold method) to 492, with a much tighter spread (425-609 vs 1210-2202 IQR).
…tes, in synthetic tuning summarize_peaks_1d() now reports amplitude_q1/q2/q3 (peak height quartiles) alongside the existing fwhm quartiles, so analyze_signal_1d() and analyze_signal_regions_1d() extract identical features -- one set for the whole signal vs. one set per region, both now including an amplitude distribution instead of a single median/weighted-mean point value. pool_signal_stats_1d()/pool_signal_stats_regions_1d() pool density the same way: instead of collapsing straight to a single across-sample median, density_q1/q2/q3 keeps the across-sample quartiles. amplitude_q1/q2/q3 pools via the same column-wise-median-across-samples treatment fwhm already used. merge_sparse_regions_1d() recombines each quantile independently (span-weighted for density, peak-count-weighted for amplitude), generalizing its existing single-value formulas. gen_synthetic_1d()/gen_synthetic_1d_regions() draw a fresh density value from density_q1/q2/q3 on every call/every region (via the new draw_from_quantiles_1d() -- a distribution-free piecewise-linear quantile function through the 3 known points), so a batch of synthetic draws reproduces real sample-to-sample density variation instead of every draw sharing one fixed value. Each individual peak likewise draws its own height from amplitude_q1/q2/q3, rather than every peak in a draw sharing one lognormal centered on a single A. Defaults preserve the old deterministic/ single-A behaviour when callers don't supply real quantiles. Updated the existing hand-built test fixtures across test-psalsa_region_amplitude.R, test-psalsa_region_anchor.R, test-psalsa_spatial.R and test-psalsa_tune_regions.R to the new density_q1/q2/q3 and amplitude_q1/q2/q3 schema.
…or scoring cap_density existed because a hard lo:hi window can't validly attribute area wherever peaks overlap. peak_area_errors_1d()'s fractional attribution (w_i(x) = pk_i(x) / peaks_total(x), using each peak's own known ground-truth shape) scores overlapping peaks correctly too, so that reason no longer applies. Verified on the real MTBLS242 dataset: with cap_density = TRUE (the old default), the peak-spacing cap bound in every region tried, capping the peak count to 3-7x below what the analyzed density actually implied -- silently overriding the density_q1/q2/q3 draw regardless of its value (every synthetic draw got the same, cap-determined peak count). With the default now FALSE, peak counts genuinely vary draw to draw, tracking the analyzed density as intended.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.