Skip to content

Feat/sensitivity improvements - #46

Merged
RobbinBouwmeester merged 51 commits into
mainfrom
feat/sensitivity-improvements
Jul 30, 2026
Merged

Feat/sensitivity improvements#46
RobbinBouwmeester merged 51 commits into
mainfrom
feat/sensitivity-improvements

Conversation

@RobbinBouwmeester

Copy link
Copy Markdown
Member

No description provided.

RobbinBouwmeester and others added 30 commits July 24, 2026 19:17
Update the guide for commit f8f6c08 and this session's validated findings:

- docs/05, docs/02: N-terminal Met-excision in the digest (DigestConfig.
  n_term_met_excision, default on), behavior, config gate, and the two tests.
- docs/13, docs/README, README: augment_library.py added as a fourth
  imported-library helper; scripts inventory 10 -> 11 (7 sidecars + 4 helpers).
- docs/18, docs/20: record the AIF E. coli findings (Met-excision closes the
  not-in-database gap 209 -> 0; missed-peptide ceiling = extract presence/apex
  then rescore then seed; the rescorer is feature-limited not training-limited;
  DIA-NN 2.2.0 comparison) with explicit single-dataset / decoy-FDR-only /
  entrapment-gated caveats.
- Correct drifted file:line refs shifted by the inserted Met-excision tests
  (digest.rs test block) plus two pre-existing imprecise config.rs citations.

Docs only. Extended feature count unchanged (381); 0 em-dashes.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Surface per-run, per-dimension calibration quality so a run's calibration can be
judged good or biased, using data the fits already compute (AlphaDIA emits an
equivalent stat table; MuMDIA had none).

- search-seed masscal.json: add ppm_residual_median and ppm_residual_mad, the
  bias and spread of the calibrant deviations after the offset correction. A
  residual median far from zero flags mass-axis curvature a single scalar offset
  cannot remove (the case for an m/z-dependent calibration later).
- rt-im-train cal.json: add rt_residual_median_s (bias), rt_residual_abs_median_s
  (typical accuracy), and rt_residual_mad_s (spread) over the training anchors.

Diagnostic only; no extraction behavior changes. Cost is on the confident-anchor
set, so it is flat regardless of variable-modification count.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…guard

Two library-hardening additions from the AlphaDIA gap analysis, both flat-compute
at high variable-modification counts.

- Fragment cardinality: a new `cardinality` column on the fragment library
  counting how many distinct precursors share each fragment m/z (0.01 Da bin).
  A low value marks a clean, quantification-friendly ion; a high value an
  interference-prone non-unique ion. Computed once at library build in both the
  imported path (import_diann_lib.py) and the native path (predict_frag.rs, with
  a deterministic ordered-map helper matched to the import binning). Additive and
  ignored by current readers, so no behavior change; it is the deterministic
  input a later interference-aware feature/quant selector will read instead of a
  runtime heuristic.
- Finite guard at the prediction sidecar boundary: fail loudly on a non-finite
  predicted iRT or fragment intensity from MS2PIP/DeepLC, which would otherwise
  silently corrupt the library and every downstream spectral-similarity feature
  (and misorder the top-N intensity sort). The rescore sidecar boundary already
  guards its feature matrix the same way.

Tests: fragment_cardinality bin/distinct-count unit test. Gate green.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…e envelope

Two AlphaDIA-derived options, both config-gated and behavior-preserving by
default (they only change results when explicitly enabled), so they land ahead of
their benchmark. Both are compute-flat with respect to variable-modification
count: the apex work is per extracted candidate and the envelope is strictly
post-FDR on quantified survivors.

- extract.apex_gaussian_sigma_scans (default 0.0 = unchanged): when > 0, smooth
  the per-scan fragment-count series with a Gaussian matched filter (radius
  3*sigma) before apex selection instead of the truncated rolling sum. Targets
  the RT-misassignment that the selected-apex-is-correct-only-about-half-the-time
  measurement attributes to apex localization. The default rolling-sum path is
  numerically identical (integer counts widened to f64).

- quant.interference_envelope (default false): apply an apex-outward
  interference-correction envelope (center_envelope_1d, a cumulative-minimum walk
  from the apex) to each fragment trace before trapezoidal integration, clipping
  co-eluting interference in the peak wings out of the reported area. Identity on
  a clean monotone peak.

Both changes are opt-in and must be entrapment or quant-accuracy validated (and
evaluated against nn_torch seed noise) before any default flip.

Tests: center_envelope_1d clip + area-shrink unit test. Gate green (130 tests).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
quant-lfq previously emitted only a protein-group by-run matrix. Add
peptide-level (by stripped base sequence) and precursor-level (by peptidoform +
charge) matrices as sibling files next to the protein output
(`<out>.peptide.parquet`, `<out>.precursor.parquet`), rolled up from the same
normalized feature-by-run matrix with the same MaxLFQ/directLFQ engine. The
protein output is byte-unchanged.

Purely additive analysis granularity, strictly post-FDR, so it changes no
identification or FDR result and its cost does not scale with variable-mod count.

Tests: base_sequence strip unit test. Gate green.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…lumn

Reflect the just-landed changes in the tracked guide:
- docs/02: extract.apex_gaussian_sigma_scans and quant.interference_envelope
  (both opt-in, default off, benchmark-gated).
- docs/12: the quant-lfq .peptide.parquet / .precursor.parquet sibling matrices.
- docs/15: the fragment library `cardinality` column.

Docs only.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…-in)

Replace the single scalar fragment mass offset with an optional m/z-dependent
correction. When search_seed.mass_cal_loess is set, search-seed fits a LOESS of
the calibrant ppm deviation versus fragment m/z and samples it onto a 25-Th grid
in <seed>.masscal.json; extract then interpolates a per-peak offset from that grid
(MassOffset) instead of the flat offset, removing any m/z-correlated curvature the
scalar offset leaves.

Default false (scalar offset, byte-identical), opt-in and benchmark-gated. The
grid is applied only when present and well-formed, so old masscal artifacts and
the scalar path are unchanged. Deterministic (calibrant pairs sorted by m/z before
the fit; grid interpolation is an ordered binary search).

Benchmark (E. coli AIF, MUMDIA_NN_SEEDS=5): +8 peptides vs scalar (within noise),
because the Orbitrap mass axis is already well de-biased (post-correction residual
median 0.0, no m/z curvature). The value is on data with real mass-axis structure
(uncalibrated instruments, drift), which the new calibration-quality residual stats
now flag.

Tests: MassOffset scalar + grid-interpolation + clamping unit test. Gate green.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ore + MBR)

Add the experiment-wide orchestrator `mumdia run-experiment`, reifying the
previously-manual multi-run flow into one subcommand and making the rescuable-tier
MBR first-class.

Flow: build/import the library once, run the per-file search chain
(convert -> seed -> optional per-run DeepLC fine-tune -> rt-cal -> extract ->
features -> compete) over N runs keeping each run's chromatograms, then one
experiment-wide rescore over all competed tables (shared classifier + FDR),
optional rescuable MBR transfer (permuted-RT decoy-transfer FDR), a schema-
preserving by-`source` split (arrow row filter), per-run quant with `q_filter=psm_q`
so each run keeps its own confident PSMs, and cross-run MaxLFQ (protein + peptide +
precursor matrices). Writes an experiment_manifest.json.

- CLI: repeat --mzml per run (>= 2), optional --run-names; library-input or --fasta.
- MBR is optional (mbr.strategy); when it accepts no transfers (e.g. every
  precursor is confident in all runs) the augmented scored table is absent and the
  flow falls back to the combined scored table.
- New mumdia-io helper write_batches: write pre-built arrow RecordBatches
  preserving their schema, used for the by-source split.

Validated end-to-end on the E. coli AIF file run twice: all stages fired,
per-run quant + LFQ + manifest produced, the 0-transfer MBR fallback exercised.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ts (opt-in)

Add an opt-in restriction of the search space to charge states a peptide can
physically carry: one proton on the N-terminus plus one per basic residue
(Arg, His, Lys). A precursor is allowed only at charge <= 1 + (#R + #H + #K); a
b/y fragment only at charge <= 1 + (basic residues within that fragment) and never
above the precursor charge. Not charge prediction: a hard chemical bound on the
enumerated charges.

- mumdia-core mass.rs: ParsedPeptidoform::basic_residue_count() and
  fragment_basic_sites(ion, ordinal) (b-ion spans the first k residues, y-ion the
  last k), with unit tests.
- config.rs: peptidoforms.charge_by_basic_residues and
  predict_frag.charge_by_basic_residues, both default false. When set, the cap
  defines the precursor range (z = 1..=1+#basic, ignoring charge_min/charge_max)
  and supersedes charge2_from_precursor_charge for fragments.
- peptidoforms.rs / predict_frag.rs: native enumeration honors the cap.
- import_diann_lib.py --charge-by-basic-residues: same rule applied to an imported
  DIA-NN library, filtering precursor and fragment rows before candidate_id
  assignment so n_fragments stays exact.
- docs/02 config tables updated.

Default off and benchmark-gated: it changes the search / training / FDR
population. Decoys inherit target composition (reverse/shift preserve R/H/K), so
the cap is symmetric across labels and target-decoy exchangeability is preserved.

Measured on the E. coli AIF file (native_tda, DeepLC fine-tune, imported library):
the cap drops 7.87% of library precursors (almost all charge-3 with a single basic
residue) and costs -5.1% confident peptides at 1% FDR (10,730 -> 10,184) with the
decoy fraction unchanged (0.98% -> 0.97%). The lost identifications are real
charge-3 single-basic peptides, so the strict cap is a net negative on this
acquisition and must not be enabled by default. It remains available for
compute-bound searches with many variable modifications, where pruning the
impossible-charge tail can be worth the sensitivity cost.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…re values

Several rescorer features divided by a value guarded only with EPS=1e-9, which is
~14 orders of magnitude below the intensity/RT scale. When the denominator is a
genuinely-absent isotope, a flat baseline, or a zero-width reference peak, the
result was ~1e15 instead of "no signal". Measured on the E. coli AIF features:
iso_ratio_1_0/iso_ratio_2_0 and the derived iso_*_ratio_dev / iso_minus_one_fraction
reached a p99 of ~1e15 (i.e. ~1% of every PSM's vector carried the garbage), and
ms1_ms2_fwhm_ratio / apex_to_boundary_ratio / profile_peak_snr / peak_shoulder_score
had ~1e10-1e16 tails.

Replace the `x/(y+EPS)` pattern with a guarded, clamped `ratio()` helper (0 when the
denominator is not meaningfully positive; clamped to a physical cap otherwise), and
bound the four other sites the same way. Isotope ratios cap at 10 (averagine
T1/T0, T2/T0 stay < 5 up to ~5 kDa), width/SNR/boundary ratios at their own physical
caps. Post-fix the same columns have sane maxima (10, 7, ...) and p99 ~3-6.

Also refresh the ms1.rs doc comments: the MS1 isotope XIC co-elution features are
populated whenever extract runs with MS1 data and window-grid chromatograms
(default on), not "unpopulated / 0.0 until persisted" as the stale comments claimed.

Impact: a correctness fix, ID-neutral on this benchmark (native_tda 10,730 -> 10,730
peptides at 1% FDR; nn_torch 10,778 -> 10,793, within replicate noise; decoy fraction
unchanged at ~0.97-0.98%). Both rescorers already absorbed the outliers via robust
feature scaling, but the 1e15 values are never correct and would mislead any
less-robust consumer (mokapot, feature-importance analysis, a future model).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…AlphaDIA #12)

Add the AlphaDIA-plan item 12 apex-isotope precursor feature: the Pearson
correlation of the observed apex isotope heights [i0,i1,i2] against the
Poisson-averagine model. Gated by `features.ms1_precursor_features` (default
false); when off it returns 0.0 so the vector is unchanged in effect.

Benchmark (E. coli AIF, native_tda deterministic, reusing one extraction): off
10,740 vs on 10,730 confident peptides at 1% FDR, decoy fraction 0.98% either way.
The feature is a net -10 (flat/marginally negative) because it overlaps the existing
`ms1_isotope_cosine_apex` (cosine vs Pearson of the same 3-vector carry nearly the
same information). It therefore stays default-off, consistent with the other
benchmark-negative opt-in features. The item's second half (ms1_prec_mass_dev_ppm,
a precursor mass-accuracy feature) is deliberately not implemented: the precursor
mass axis was already shown flat on this instrument (0.79 ppm systematic offset,
tolerance already at its optimum), so a mass-deviation feature is certain to be flat
here and does not justify the observed-centroid extraction plumbing.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…nfig)

First increment of AlphaDIA plan #7 (top-K peak promotion). The measured lever: the
selected chromatographic apex is correct only ~48-52% of the time while the correct
peak is in the top five ~86-88%, so promoting alternate peaks through features/rescore
lets the scorer pick the correct-but-not-apex peak. Unlike every other AIF lever tried
this session (all flat), this gives the rescorer a candidate it otherwise never sees.

This increment lays the contract without any behaviour change:
- config: extract.promote_top_peaks (default 1), alt_peak_min_area_frac (0.10),
  alt_peak_min_separation_s (5.0). At the default 1 only the selected apex is emitted.
- extract: every emitted PSM carries a peak_rank (CandOut.peak_rank, 0 = selected apex)
  and psms_extracted gains a peak_rank column. PSMS_EXTRACTED schema 1 -> 2.

Verified byte-identical on the E. coli AIF chain (native_tda): peak_rank present and
all-zero over 47,484 rows, 10,740 confident peptides at 1% FDR (unchanged), decoy
0.98%. Subsequent increments thread peak_rank through features/compete/rescore (with a
per-candidate q collapse so the decoy null is not K-inflated) and then emit ranks >= 1;
that emission is behaviour-changing and stays entrapment-gated. Plan in topk_plan.md.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…andidate collapse

Second increment of top-K peak promotion (AlphaDIA #7). Carries the peak_rank
contract through the rest of the pipeline and adds the FDR-critical collapse, all
byte-identical while promote_top_peaks = 1.

- features: peak_rank is passed through as a bookkeeping column (never a predictive
  feature); missing in pre-v2 inputs -> 0.
- compete: peak_rank joins the competition group key (4-tuple), so peaks of one
  candidate compete only within their own rank and a weaker peak cannot eliminate a
  sibling's better peak on prelim score before rescore chooses. At rank 0 everywhere
  the grouping is unchanged. Carries peak_rank forward. PSMS_COMPETED 2 -> 3.
- rescore: after scoring and before every q context, collapse to the best-scoring
  peak per (source, candidate_id) so the rescorer (not the up-front apex) selects the
  peak and the terminal q-null is exactly one row per candidate; promoting K peaks
  therefore does not K-inflate the decoy null. Decoys collapse by the identical rule,
  preserving target/decoy exchangeability. Guarded (best.len() == n at K=1 -> no-op).
  Emits selected_peak_rank. PSMS_SCORED 3 -> 4.

Verified byte-identical on the E. coli AIF chain (native_tda): peak_rank present at
features/compete/scored, selected_peak_rank all-zero, 10,740 confident peptides at 1%
FDR (unchanged), decoy 0.98%. Increment 3 will emit ranks >= 1 (extract closure ->
Vec<CandOut>), which is the behaviour change and stays entrapment-gated.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Third increment of top-K peak promotion (AlphaDIA #7): the behaviour change. The
per-candidate extraction closure now returns Vec<CandOut> instead of one record. At
promote_top_peaks = 1 it returns exactly the rank-0 record (byte-identical). At K > 1
it additionally enumerates alternate chromatographic peaks and emits up to K-1 extra
rows.

- The MS1 apex-isotope computation is factored into an ms1_at(rt) closure so rank 0
  and each promoted peak compute MS1 evidence at their own apex.
- Alternate peaks come from enumerate_peaks over the distinct-fragment count profile
  (the same interference-resistant breadth profile as the diagnostic sidecar). The
  envelope holding the selected apex is excluded, and each alternate must clear the
  area-fraction gate (alt_peak_min_area_frac), the apex-RT separation gate
  (alt_peak_min_separation_s), and the matched-fragment floor.
- Each alternate re-slices only its own apex-dependent scalars (apex_rt, apex
  intensity, matched count, MS1); it carries an empty chromatogram set and shares the
  candidate's chromatograms, which are emitted once on rank 0 and looked up by
  candidate_id downstream, so the largest artifact does not multiply with K.

Benchmark on the E. coli AIF file (native_tda, deterministic, one extraction per K):
K=1 is byte-identical (10,740 peptides at 1% FDR, 47,484 rank-0 rows). K=2 emits
26,170 alternate rows, K=3 emits 47,706; the rank-0 row count is unchanged at every K.
The decoy fraction stays 0.97-0.98% at every K, so the per-candidate best-peak
collapse in rescore keeps the null un-inflated. The rescorer selects an alternate peak
for ~2,800 (K=2) / ~3,600 (K=3) confident targets. Confident precursor count rises
(10,615 -> 10,942 -> 11,001) while the base-peptide count is flat-to-slightly-lower
(10,740 -> 10,706 -> 10,590). Part of the precursor rise is compete letting different
charge/mod siblings win at different ranks, so the gain is not yet trustworthy: this
stays default-off pending entrapment validation (increment 4) on >= 2 acquisition
contexts, which is the FDR-validity gate the decoy-fraction check alone cannot provide.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…(1/8)

Begin the modular, interswappable fragment-competition component from
fragment_competition_strategies.md. Competition is now a swappable part: select a
method via config, compose weight cues, or turn it off entirely (the default).

- config: new PeakClaim::CoelutionMultiCue method + a composable ClaimCues struct
  (extract.claim_cues), each cue label-blind and defaulting off (weight 1.0).
- extract: the two-pass arbitration's scalar elution-profile-height `ph` is replaced
  by a modular per-claimant `weights` vector = ph x product of enabled cues. Every
  non-MultiCue method uses cue = 1.0, so `weights` reduces EXACTLY to `ph` and the
  arbitration is bit-identical to before (CoelutionWinner/Proportional/WinnerMargin
  and the default None are unchanged).
- First cue S3 `mz_close`: exp(-(ppm_err/sigma)^2) on the observed peak's signed ppm
  offset from a claimant's predicted fragment m/z. A shared peak sits at the true
  owner's m/z, so the sub-tolerance offset (which engines use only as a binary gate)
  becomes a novel apportionment weight - high value where ~98% of DIA fragment m/z
  collide within tolerance.
- CoelutionMultiCue ships NON-DESTRUCTIVELY: its cue-weighted split flows into the
  contested/apportioned features (visible under emit_contested_features) without
  altering extracted intensities. A destructive reassign variant is separately
  benchmark+entrapment-gated per the report.

Verified on the E. coli AIF chain (native_tda): default (peak_claim=None) byte-identical
at 10,740 peptides; CoelutionMultiCue + mz_close non-destructive gives 10,733 (feature-
only, within noise) with the decoy fraction unchanged at 0.98% and peak_apportioned_frac
now cue-weighted (100% populated). native_tda barely weights that feature; the nn_torch +
entrapment benchmark (increment 8) is the real test. Plan for the remaining 7 strategies
in frag_competition_impl_plan.md.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Second composable claim-weight cue for the modular fragment-competition framework.

- config: ClaimCues.rt_prior + rt_prior_tau_s (default 30 s).
- extract: thread rt_cal (calibrated predicted RT) into extract_twopass_windows; the
  rt_prior cue multiplies each claimant's weight by exp(-(rt - rt_pred)^2 / 2 tau^2),
  down-weighting a co-isolated interferent whose DeepLC-predicted RT is far from the
  current scan. Composes multiplicatively with mz_close. No-op where rt_pred is unset;
  label-blind; default off.

Verified on the E. coli AIF chain (native_tda, non-destructive CoelutionMultiCue +
emit_contested): mz_close+rt_prior gives 10,734 peptides at 1% FDR (baseline 10,740,
feature-only/neutral), decoy 0.98%, peak_apportioned_frac cue-weighted. Default
(peak_claim=None) remains byte-identical.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Third composable claim-weight cue: cross-dimension MS1 precursor support (S4).

- config: ClaimCues.ms1_support.
- extract: thread ms1_scans + ms1_rts into extract_twopass_windows; the cue looks up
  the claimant's OWN precursor isotope envelope (mono + a plausible +1/mono ratio) at
  the nearest MS1 scan and down-weights a claimant with no real precursor signal
  (mono absent -> x0.5; present but implausible ratio -> x0.75). A shift/reverse decoy
  has a defined precursor m/z but no co-eluting MS1 precursor, so this starves its MS2
  claim through an orthogonal dimension it cannot fake. Composes multiplicatively with
  mz_close and rt_prior; no-op when no MS1 is provided; label-blind; default off.

Verified on the E. coli AIF chain (native_tda, non-destructive CoelutionMultiCue +
emit_contested): mz_close+rt_prior+ms1_support gives 10,727 peptides at 1% FDR (baseline
10,740, feature-only), decoy 0.98%, peak_apportioned_frac reflecting all three cues.
Default (peak_claim=None) byte-identical.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…on 4/8)

S2: a deterministic fixed-point that fixes a concrete flaw in the co-elution profile.
The profile is built from FULL peak intensities, so a borrowing candidate's per-scan
profile height is inflated by the very peaks it borrows, which then biases every
profile-height-weighted decision (and the apportioned/contested features).

- config: ClaimCues.apportion_em_iters (default 0 = off).
- extract: under CoelutionMultiCue with N > 0, run N refinement passes between the
  profile build and the final arbitration. Each pass re-seeds every candidate's profile
  from its cue-weighted APPORTIONED share (not the full peak). An uncontested
  single-claimant peak contributes full intensity every iteration, an immovable anchor,
  so candidates with real unique evidence converge high and pure borrowers collapse
  toward zero. Deterministic (fixed N, ordered f32 reductions, sorted claimant order);
  non-destructive; label-blind; no-op at N=0.

Verified on the E. coli AIF chain (native_tda, non-destructive, all cues + EM=3): 10,731
peptides at 1% FDR (baseline 10,740, feature-only), decoy 0.98%, peak_apportioned_frac now
borrow-corrected. Default (peak_claim=None) byte-identical.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…g (T2 foundation)

Foundation for the spectrum-centric deconvolution strategies (D2 NNLS, D3 elastic-net,
S1 ridge) of the fragment-competition report. Solves min_{x>=0} ||A x - b||^2 where A is
the co-isolated candidate x fragment-channel design matrix and b the observed spectrum;
x_p is candidate p's interference-corrected abundance (the CHIMERYS coefficient idea,
reimplemented clean-room).

Lawson-Hanson active set with the passive-set least squares solved by RIDGE-regularized
normal equations (A_P^T A_P + lambda I) z = A_P^T b via a hand-rolled Cholesky. Ridge
keeps the Gram strictly PD, so the solve is unique and bit-stable even under the ~98%
column collinearity of wide-window DIA - the deliberate divergence from CHIMERYS's
LASSO+AICc. All reductions run in fixed index order; no BLAS/LAPACK (Windows/OneDrive
build bars C deps); no randomness. Deterministic by construction.

Tests: identity recovery, negative-clamp-to-zero, two-column apportionment, collinear
columns bit-stable across runs under ridge, empty problem. The demix wiring (assemble the
per-window apex design matrix, emit demix features, plumb to rescore) follows.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…tion 5/8, D2)

Wire the deterministic NNLS solver into a non-destructive demix feature family - the
CHIMERYS-style spectrum-centric competition, made chromatographic and clean-room.

- extract: demix_at_apex assembles the co-isolated candidate x fragment design matrix at
  each accepted candidate's apex scan (rows = observed peaks, columns = candidates that
  claim a peak in sorted cid order, D[peak,cand] = predicted intensity), solves
  min_{beta>=0}||D beta - y||^2 (crate::solve::nnls, ridge), and emits three features:
  deconv_explained_frac (joint residual-explained), deconv_active (this candidate survived
  the active set), deconv_share (its fraction of the total demixed abundance). Gated by
  extract.emit_demix_features (+ demix_lambda, demix_max_candidates); intensities untouched;
  a per-run rt_bits->scan-index map is built only when demixing is on. Deterministic:
  peaks in scan order, columns sorted by cid, ordered reductions.
- features: new `demix` family (3 columns) reading the values from psms via Evidence; 0
  when demixing is off, so the vector is unchanged by default.

Verified on the E. coli AIF chain (native_tda, non-destructive): demix-on gives 10,725
peptides at 1% FDR (baseline 10,740, feature-only), decoy 0.98%. The features DISCRIMINATE
clearly - deconv_share target mean 0.122 vs decoy 0.040 (a decoy borrower is explained
away by better-supported co-eluters), deconv_active 0.75 vs 0.60. native_tda barely weights
new features, so IDs are flat; the nn_torch + entrapment benchmark is the real test.
Default (emit_demix_features=false) byte-identical.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…(competition 6/8, D3)

D3's key novel contribution: an explicit identifiability signal that no incumbent engine
emits. The ridge in the NNLS solver already provides D3's L2 collinearity handling (mass
shared smoothly among collinear candidates rather than coin-flipped); this exports the
"how identifiable was this?" measure that ridge alone cannot report.

- extract: demix_at_apex additionally computes, for the target candidate's design column,
  the maximum cosine similarity with any co-isolated rival column (O(m n), reusing the
  assembled matrix). Near 1 = near-degenerate with a rival, so the demix split is arbitrary
  and the coefficient should be distrusted. Emitted as deconv_max_collinearity under the
  same emit_demix_features gate.
- features: the demix family gains the 4th column, read via Evidence, clamped [0,1], 0 when
  demixing is off.

Verified on the E. coli AIF chain: deconv_max_collinearity is populated (60% nonzero,
bounded 0-1, target 0.33 / decoy 0.29). Default (emit_demix_features=false) byte-identical.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…petition 7/8, D1)

D1: a distinct deconvolution mechanism from the NNLS solve. Estimate each co-eluter's
abundance from the channels it ALONE claims (its unique ions, a_p = median y/D over
unique rows), subtract those contributions from the candidate's channels, and measure the
fraction of the candidate's observed apex intensity that survives. This sidesteps the
shared-peak circularity by seeding abundances from unique ions only, and distinguishes a
real second peptide (keeps its signal) from a pure borrower (subtracted toward zero).
Computed inside demix_at_apex, reusing the assembled co-isolated design matrix (no extra
probing); deterministic (unique-row medians over sorted ratios, ordered reductions).

- extract: demix_at_apex returns shadow_kept in addition to the demix/identifiability
  outputs; emitted as shadow_kept_frac under emit_demix_features.
- features: the demix family gains the column (read via Evidence, clamped [0,1]).

Verified on the E. coli AIF chain: shadow_kept_frac discriminates clearly (target mean
0.658 vs decoy 0.486). All five demix features now separate targets from decoys
(deconv_share 0.12/0.04, deconv_active 0.75/0.60, deconv_explained 0.44/0.34, shadow_kept
0.66/0.49, deconv_max_collinearity 0.33/0.29). Default (emit_demix_features=false)
byte-identical; native_tda IDs neutral, nn_torch + entrapment is the real test.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ion, destructive A)

The impactful half of the modular fragment-competition component: competition that
rewrites the extracted evidence, not just emits features. When claim_cues.reassign is set,
CoelutionMultiCue redistributes each shared peak winner-take-all on the composite cue
weight (mz_close x rt_prior x ms1_support), so the competed intensities feed EVERY
downstream feature (co-elution, spectral, mass accuracy, ...) rather than only the
apportioned/contested columns.

- config: ClaimCues.reassign (default false).
- extract: MultiCue joins the reassign path when the flag is set; its reassign arm is
  winner-take-all on the composite weight (CoelutionWinner is the same rule on the plain
  profile height).

Verified on the E. coli AIF chain (native_tda): destructive MultiCue gives 10,326 peptides
at 1% FDR vs 10,740 baseline (-3.9%), decoy 0.98% (FDR still valid). Winner-take-all is
aggressive on chimeric wide-window DIA - it strips intensity from real co-eluting peptides,
exactly the failure mode the report predicts for hard redistribution on wide windows. The
smoother demix beta-apportionment (splits rather than strips) is the next, more principled
destructive mode. Default off; entrapment-gated per CLAUDE.md.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…x (destructive B)

The smooth, principled destructive mode: split each shared peak by the joint NNLS
solution rather than stripping winner-take-all. New PeakClaim::CoelutionDemix (two-pass,
always redistributes). At each scan, assemble the co-isolated candidate x fragment design
matrix, solve min_{beta>=0}||D beta - y||^2 (crate::solve::nnls, ridge), and give each
claimant of a peak the share beta_c * D[peak,c] / sum_p beta_p * D[peak,p]; the competed
intensities then feed every downstream feature. Contested features (won/lost/apportioned)
are accrued from the same split. Deterministic (sorted candidate columns, ridge NNLS,
ordered reductions). Default off.

KNOWN COST: this solves one NNLS PER SCAN. On the full E. coli AIF file (465,806 scans) the
extract does not complete in a practical time - the report's predicted "heaviest tier;
bound T; solve only at apexes" outcome. The mechanism is correct (shares the unit-tested
solver and the same assembly that produced the strongly target/decoy-discriminating demix
FEATURES), but as a full-file per-scan solve it is impractical for wide-window DIA. The
required optimization is to restrict the solve to accepted-candidate apex neighborhoods (or
per-window elution regions) and reuse across nearby scans, keeping the redistribution while
bounding the solve count. Shipped as an opt-in foundation with the cost caveat; the
apex-restricted variant is the follow-up. The gentler alternative already validated is the
demix FEATURE path (emit_demix_features, non-destructive) which feeds the rescorer the
interference-corrected coefficients without the per-scan redistribution cost.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ubtraction mode

Two follow-ups to the destructive fragment-competition work.

1. demix_scan_stride: the per-scan NNLS of CoelutionDemix is impractical over a full
   wide-window run (~465k scans). Solve the NNLS every Nth scan and reuse the candidate
   abundances to apportion the intervening scans, forcing a re-solve whenever a new
   candidate enters the co-isolated set (so accuracy holds where the population changes).
   Each scan's own peaks + predicted intensities still drive the split. Default 1 (solve
   every scan, unchanged). Per-window state, so par_iter stays deterministic.

2. PeakClaim::CoelutionShadow: the gentle, solver-free destructive mode. Each candidate's
   abundance is estimated from the channels it ALONE claims (a_p = median observed/predicted
   over its unique ions); every candidate then keeps max(0, y - sum_{p!=c} a_p*D) at each of
   its channels. Several real co-eluters can both retain signal at a shared peak (unlike
   winner-take-all), and no solve is needed. Deterministic.

Benchmark on the E. coli AIF file (native_tda, all FDR-valid at decoy ~0.98%):
  baseline                    10,740
  winner-take-all (destr A)   10,326  (-3.9%, aggressive - strips real co-eluting signal)
  shadow (subtract)           10,660  (-0.7%, extract 12s)
  demix beta-split (stride=8) 10,763  (+0.2%, extract 613s)
The smooth NNLS beta-apportionment is the only destructive mode that does not damage IDs on
chimeric wide-window DIA - it reshapes every downstream feature yet stays neutral/safe. On
AIF all sit near the data/library ceiling; the payoff is narrow-window acquisitions where
interference (hence competition) actually bites. Every destructive mode is default off and
entrapment-gated per CLAUDE.md.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Adds Farnesyl (204.187801), GeranylGeranyl (272.250401) and Hydroxyfarnesyl
(220.182716) to the `unimod_mass` table, and teaches the DIA-NN importer to map
and KEEP those UniMod ids instead of dropping their precursors.

The importer previously mapped only UniMod 4/35 and discarded every precursor
carrying anything else, so prenylated peptides could never enter a search at all.
Replacement is substring-exact including the closing paren, so "(UniMod:4)" does
not match inside "(UniMod:44)".

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Extract is allocation-heavy inside rayon workers and measured only ~1.07x parallel
scaling under the Windows system allocator's shared heap lock. Swapping in mimalloc
(per-thread arenas) plus thin LTO and codegen-units=1 gives ~1.98x on extract with
BIT-IDENTICAL output (2,244,996 accepted / 33,553,770 chromatograms; 21/21 columns
compared equal).

Measured with 3 ALTERNATING passes, fresh-binary-first so the old binary got the
warm-cache advantage: ratios 2.02 / 1.92 / 1.99 (means 107.8s vs 213.5s).

Benchmarking note: the SAME binary runs 137s cold vs 77s warm on a 3 GB library, so
never trust absolute numbers from two arms run back-to-back -- only same-pass ratios
(both arms shift together and the ratio is invariant).

Deviation from the repo's pure-Rust dependency preference: mimalloc compiles C. It is
kept because the gain is large, reproducible and behaviour-neutral; see
optimization_audit.md.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
table::nrows() reads the row count from the parquet footer. Five run.rs sites were
calling Table::read(path)?.nrows, which fully decodes the file -- including the
fragment library, hundreds of millions of rows -- purely to learn its length.

index.rs: the typed column getters return owned Vecs, so the decoded Arrow batches are
dead weight for the rest of Library::load; drop both tables right after extraction
(the fragment table alone is ~23 GB of batches at 657M rows). Fragment names are now
MOVED out of the source vector instead of cloned, so only one copy of every name is
live at peak, and the emptied source plus the per-candidate grouping are released
before the index build allocates.

Extract output verified bit-identical after these changes (21/21 columns).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
… its write

write_pin now takes one slice reference per feature column, resolved once, instead of a
transposed Vec<Vec<f64>>. At 1.5M rows x 387 features that transpose allocated ~4.6 GB
and did ~580M string-keyed HashMap lookups. Byte output is unchanged.

Also adds `features.emit_pin` (default TRUE, so no artifact silently disappears). No
MuMDIA stage reads run.pin -- rescore builds its own PIN for the sidecars -- so setting
it false skips a ~5.4 GB text write.

FeatureSchema::read now degrades gracefully: if the .schema.json companion is missing or
corrupt, the feature list is reconstructed from the artifact's own parquet schema by
excluding the known metadata columns, with a warning. A missing companion previously
aborted the run, which happened for real when a competed table was rewritten by an
external tool that did not copy the sidecar.

The metadata exclusion list is verified against a real artifact to reproduce the recorded
feature_columns exactly (387, same order). Two traps it encodes: `charge` IS a feature
(carried as f64) while elution_lo/elution_hi are quant metadata -- misclassifying either
would silently change the rescorer's trained population.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…rker child

Three defects, two of them observed live this session.

1. Fixed sidecar filenames. `rescore.pin` / `rescore_sidecar_out.parquet` /
   `*.feat.mm` were the same paths for every invocation, so two concurrent rescores
   clobbered each other and one died with a bare OSError [Errno 22] naming a file it did
   not own. Paths are now keyed on the output artifact stem plus the PID.

2. Orphaned worker poisoned later runs. Killing the mumdia parent left the Python worker
   alive holding its multi-GB feature memmap, so the file could not be deleted and EVERY
   subsequent rescore failed on the stale lock -- with an error naming the path, not the
   PID holding it. The child is now owned by a guard that kills it on drop/unwind.

3. The PIN was accumulated in ONE un-reserved String before being written: a >5 GB
   allocation plus realloc churn held entirely in RAM at ~1M rows x 387 features. It is
   now streamed through a BufWriter. Byte output unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
RobbinBouwmeester and others added 21 commits July 28, 2026 08:02
…early-stop (3.11x)

Production measurement on the full HYE PIN (959,426 rows x 387 features, 3.5 GB):
  before  43.5 min  psm@1%=68,686  peptides@1%=61,262  decoy 0.99%
  after   14.0 min  psm@1%=68,588  peptides@1%=61,172  decoy 0.99%
  => 3.11x, -90 peptides (-0.15%, within NN run-to-run noise), decoy% identical.

The MLP is tiny (387->128->64->1, 58k params) but runs ~160,500 times, so cost is step
count and data movement, not model size -- cutting the model to 15k params changed
runtime by 1.02x, and Adam foreach/fused/SGD are all within noise.

- PREGATHER (the win, 2.72x alone): the training row set is fixed for all EPOCHS, so its
  features are gathered ONCE into a contiguous tensor rather than fancy-indexing the full
  matrix per minibatch. MUMDIA_NN_PREGATHER_GB=8 bounds it and falls back to the old path
  above that, so the streaming backend keeps its low-RAM guarantee.
- MUMDIA_NN_THREADS=16 default. Note it is a REGRESSION on its own (0.94x) and only pays
  off after the pregather; the MLP saturates at 16 threads and 32 is slower than 16.
- Churn-tolerance early stop (MUMDIA_NN_EARLY_STOP_TOL, default 0.01). Exact equality
  never fires: dropout plus retraining flips a few borderline PSMs every iteration
  forever. Worth 1.59x on a 40k pool but INERT at production scale, where churn bottoms at
  1.21% while the positive set is still growing (38,950 -> 47,231) -- correctly declining
  to stop while the model is still learning.
- Resolve the memmap path absolutely (it was relative to the caller's cwd) and log which
  backend was chosen and why, since the auto-threshold is on PIN text size and silently
  switched behaviour at ~1M PSMs.

Rejected after measurement: batch 4096->16384 is 1.43x but costs 14% of peptides.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…face inert config

Speed (all bit-identical: per-row/per-candidate only, no cross-row float reduction, and
rayon's indexed collect/stable sort preserve order):
- predict-frag: native intensity and native iRT assignment now par_iter_mut, and the
  precursor-m/z sort is par_sort_by (stable, so ties keep input order). Three whole-library
  serial loops in a file that already used rayon everywhere else.
- quant: the per-candidate peak_window phase is now par_iter, collected into a Vec and
  folded into the BTreeMap afterwards so the map is built from a deterministic sequence.
- fragment_cardinality: replaced BTreeMap<bin, BTreeSet<candidate>> with a sort of
  (bin, candidate) pairs plus a single run-scan. Same values, same determinism, flat memory
  instead of a tree node and a set per bin over ~86M fragment rows.
- peptidoforms: dropped the `seen_forms` HashSet that hashed and cloned every emitted
  ProForma string. It existed to guard one degenerate case -- a config listing the same
  variable mod twice on a residue -- which is now handled by deduplicating the handful of
  alternatives per site instead. Every combination is unique by construction once the
  alternatives are.

Correctness/provenance:
- run: fold the conversion caps (max_spectra, top_peaks_*) into the convert artifacts'
  config_hash exactly as the standalone `convert` subcommand does. Two runs with different
  caps previously recorded the same hash, disagreeing with the standalone entry point.
- compete: mode=unique_evidence silently behaved like mode=none whenever only the raw
  matched-count fallback was available, because extract's presence_min_matched (3) already
  exceeds unique_evidence_min_fragments (2), so every candidate passed and nothing was
  removed. It now warns and names the fix.
- config: warn when mbr.rt_window_s / mbr.decoy_transfer / mbr.requant_all are set, since
  no stage reads them yet, and when a non-none mbr.strategy is chosen, since only
  none-vs-not-none is ever tested. They are documented as functional but are inert; warning
  beats rejecting because deny_unknown_fields means deleting the fields breaks old configs.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ouping; window-relative seed scratch

Library load and seed scoring were paying for structures the default backend never reads.
Extract on the 13.6M-precursor focused library: 137s -> 84s, output BIT-IDENTICAL
(2,244,996 rows, 21/21 columns compared equal).

- The bucketed `page_search` index is now built only when the `bucketed` matcher is
  actually selected. On the default `fragindex` path it was a full (f32,u32,f32) copy of
  every library fragment, a global sort of it, and three more full arrays -- all dead,
  since `page_search` is never called. It is still built for `matcher = bucketed`, which
  extract and search-seed do use as a fallback, so this is gated on the matcher rather
  than removed.
- Fragment grouping by candidate_id uses a counting sort into two flat arrays instead of
  `Vec<Vec<usize>>`, which performed one heap allocation per candidate (54.8M on the
  profiled library). Scattering in ascending row order preserves stored order, so the
  resulting per-candidate layout is unchanged.
- `SeedScratch` is indexed window-relative. Sized by library candidate count it cost
  16 B x n_cand PER rayon worker (877 MB/worker at 54.8M candidates) while a worker only
  ever touches one isolation window. The caller now sizes it to the widest window, and the
  arrays grow on demand so an underestimate cannot panic.

Still outstanding from the audit: interning `Library::frag_name` to `Vec<u16>` + a
dictionary. That is the largest remaining RAM item (a `Vec<String>` is ~24 B of struct per
fragment regardless of content, ~16 GB per copy at 657M fragments), but it changes a public
field that extract reads, so it wants its own verified change.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…gment)

`Library::frag_name` was a `Vec<String>`, which costs ~24 bytes of String struct per
fragment BEFORE any text -- about 16 GB per copy at 657M library fragments -- to store
names drawn from a tiny repeating vocabulary (b1, y7, y12^2, ...). It is now a
`Vec<u16>` of dictionary indices plus a `Vec<String>` dictionary, resolved via the new
`Library::frag_name_str`.

Only one consumer ever read the names (the chromatogram emission in extract); every other
`cand_frags` caller discarded them, so the API change is narrow. `cand_frags` now returns
`&[u16]` for its third element and the String is materialised only for the emitted row.

Fails loudly if a library ever contains more than u16::MAX distinct fragment names rather
than silently truncating the id.

Verified on the 13.6M-precursor focused library: extract 84s -> 71.5s, and BOTH outputs
bit-identical -- psms 2,244,996 rows x 21 columns and the 33,553,770-row chromatogram
table x 7 columns (which carries the frag_name strings this change touches).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…t builders

Table::read_cols(path, &["a","b"]) builds a parquet ProjectionMask so unselected columns
are never fetched or decoded. Wired into the two chromatogram readers: quant needs only 4
of the artifact's 7 columns (it was decoding and holding frag_mz, frag_obs_mz and
predicted_intensity for nothing) and features names all 7 explicitly, so a future added
column is not silently decoded on that hot path. The chromatogram table is the largest
artifact in a run -- tens of millions of rows with two big list columns.

A projection-equivalence test caught a real bug in the first version: the schema was taken
from the builder, which reports the FULL file schema, so under a projection it disagreed
with the batches (which carry only the selected columns) and the typed getters resolved a
name to the wrong column index -- i.e. silently reading the wrong data. The schema now
comes from the reader. The test asserts read_cols matches read value-for-value across
u32/f64/str/list columns and that an unprojected column is absent rather than zero-filled.

List builders (ListF32, ListF64, LargeListF32 -- the last carries the chromatograms) now
reserve the exact total value count instead of reallocating as they grow, and CONSUME the
source rows so each inner Vec is freed as it is copied rather than keeping the whole
Vec<Vec<f32>> alive beside the finished array.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…hains

`run-experiment` executed its N per-run chains strictly sequentially even
though the runs are fully independent: `process_run` derives every path from
its own per-run output directory and shares nothing but the read-only library.

Add `experiment.parallel_runs` (default 1) to run chains concurrently. The
default reproduces the previous sequential loop exactly -- it is the same code
path, so opting out is not merely equivalent but identical. Concurrency is
bounded rather than a flat `par_iter` because each concurrent run holds its own
extraction working set (tens of GB on a large library), so the practical
ceiling is memory, not cores; an unbounded fan-out over an 83-file experiment
would exhaust RAM long before it saturated the CPU.

Output ordering is preserved: the downstream combined rescore keys rows by
`source` index, chunks are processed in index order, and rayon's indexed
`collect` preserves within-chunk order, so completion order never reaches the
output.

Verified on a two-run experiment (two disjoint mid-gradient slices of
LFQ_Orbitrap_AIF_Ecoli_01, DIA-NN library, native_tda rescoring and DeepLC
fine-tuning disabled so the reference is deterministic):

  all 30 parquet/tsv outputs byte-identical between parallel_runs=1 and 2
  remaining diffs are elapsed_ms, output paths, and config_hash only
  per-run section 19.7s -> 11.4s (1.73x on 2 runs)

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`quant` reuses the identification apex carried in `psms_scored`, and falls back to
detecting an apex from the chromatogram when the column is absent (pre-v3 scored
artifact) or a row's apex is null / non-finite.

That fallback is not equivalent to the reused apex. Re-detection reproduces the
identification's apex only about half the time (CLAUDE.md: "correct/strongest only
about 48-52% of the time"), so a quantity integrated around a re-detected apex can
belong to a different peak than the one that was identified. The fallback was
silent, ungated, and its comment named schema v3 while the current scored schema
is v4.

Warn when the column is missing, and record `apex_rt_column_present` plus
`candidates_with_scored_apex` in the quant artifact report so a reader can tell
which apex source a quantity actually used. Behaviour of the numeric path is
unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
… hash

`run_experiment` passed the bare config hash to `convert`, while both the
standalone `convert` subcommand and single-run `run` fold `max_spectra` and
`top_peaks_ms2` into it. Those caps change the spectra output but are not part of
the config, so two experiments converted with different caps recorded identical
convert provenance -- and disagreed with the other two entry points for the same
inputs.

This is one of the divergences between `process_run` and the `run.rs` chain it was
copied from; the copies remain to be deduplicated.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…uffers

The Lawson-Hanson passive solve rebuilt `A_P^T A_P` and `A_P^T b` from scratch on
every active-set move, an O(m k^2) pass over all m rows per move, and allocated a
fresh passive-column list, solution vector, Gram, right-hand side and Cholesky
factor each time.

Compute `A^T A` and `A^T b` once per solve and gather the passive sub-Gram out of
them, and hold every working buffer in an `NnlsScratch`. The precompute is
O(m n^2), i.e. at most the cost of the single move where all n columns are
passive, and callers bound n by `extract.demix_max_candidates` (64).

This is exactly bit-identical, not merely algebraically equal: each Gram entry is
still summed over rows in ascending order, so `g[r][c]` gathers the same terms in
the same order the per-move accumulation produced.

The gradient `w = A^T (b - A x)` is deliberately NOT rewritten as
`A^T b - (A^T A) x`. Those are algebraically equal but differ in the last bits,
and `w` decides which column enters the passive set -- a one-ULP difference near
the activation tolerance, or between two nearly-equal gradients, would pick a
different column and change fragment apportionment, hence extracted intensities
and identifications. It keeps the explicit residual form.

Verified by a new test that keeps the previous implementation verbatim as a
reference and asserts bit-equal output (`to_bits`) over 324 problems: 9 shapes x
12 sparsity/collinearity patterns x 3 ridge values, including duplicated columns,
all-zero columns, negative `b`, and m < n.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…r PSM

Interference features 15 (profile_corr_full_vs_peak_delta), 16
(frac_frag_ref_corr_below_0_5), 17 (explained_apex_intensity_frac) and 18
(apex_purity) each evaluated `pearson(&e.traces[f], r)` in their own loop over the
same (trace, reference-profile) pair, so the identical correlation was computed
four times per fragment per PSM. Three independent audit agents flagged it.

Hoist it into one `corr_ref` vector built once per PSM. This strictly removes work:
feature 18 already sweeps all k fragments, so all k correlations were computed
regardless, and features 15-17 recomputed the `matched` subset three more times.
Net saving is 3 * |matched| correlation evaluations per PSM.

Exactly value-preserving: `pearson` is pure, and the `0.0` filler for a fragment
without a trace is unreachable -- features 15-17 iterate `matched`, whose
construction already requires `f < e.traces.len()`, and feature 18 keeps its own
bounds check.

Verified: features.parquet byte-identical (blake2b 5224424578dda6e1) over 30,031
rows x 387 Extended features.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ot per candidate

`demix_at_apex` assembled the co-isolated design matrix from every peak of a
candidate's apex scan and ran the NNLS, but the candidate id only selected a column
and gated an early return: the matrix, the observed vector, the solution, the
explained fraction and the unique-ion abundances are all functions of the SCAN.
Every candidate apexing in the same scan therefore re-probed all of that scan's
peaks and re-solved the same problem. Two audit agents flagged it independently.

Split it into `demix_solve_scan` (per scan) and `demix_features_for` (per
candidate), and run it as a second pass over the collected rows grouped by resolved
apex-scan index. The group key is the scan, not the RT, because two candidates
sharing an apex RT in different isolation windows resolve to different scans. Rows
are patched by candidate id, so the result order is untouched.

Verified byte-identical (psms_extracted and chromatograms, blake2b
7bf04c61b856be83 / e6dfae90e97f7995) on a competition-config extract
(peak_claim=coelution_multi_cue, emit_demix_features=true,
emit_contested_features=true), where it cut the solve count from 30,031 to 5,527.

Scope of the win, stated honestly: the reduction is #candidates / #distinct apex
scans, which on that small input is only 5.4x, and demix there is ~0.5% of the
stage (7,593 ms vs 7,553 ms with the features off), so no wall-clock change is
visible at that size. The ratio grows with candidate density per scan; a
production-scale attribution run is in progress.

Also carries the `probe_peak` per-window narrowing cache (`WindowNarrow`) and the
`Prober` bundle it needed. The narrowing measured NEUTRAL on the default
accumulate path (before 78.0/77.1/77.7 s vs after 79.1/76.9/73.0 s, three
alternating passes on a 13.6M-candidate extract) -- the binary searches walk the
same `post_cand` cache lines the emit loop already reads, so they were close to
free. It is retained only because the two-pass path probes each peak
1+em_iters+1 times; if that measurement is also neutral it should be reverted.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
… items

An independent verification pass over the optimization audit found four items
credited as done that were only half done. Each sibling path is now fixed:

quant Phase 2 was still serial. Only the peak-window phase had been
parallelized; the integration loop that walks every chromatogram sample of every
fragment trace was untouched, and it is per-candidate independent by the same
argument. Parallel across candidates with an indexed `collect`, then merged in
`cand_rows` (BTreeMap) order, so the maps fill exactly as the serial loop filled
them. Verified byte-identical: peptide, protein-group and fragment quant outputs.

predict-frag was parallelized in the Native arms only. The MS2PIP arm calls the
native predictor per row as the charge-2 fallback and then does two per-charge-group
max-normalizations -- real CPU work, not sidecar wait -- and it is the arm the
shipped `docker/config.dia.json` selects (`predictor: ms2pip`), so neither
parallelized branch ever ran in that image. Now `par_iter_mut().enumerate()`;
each row writes only its own `frag_int`.

Wide-table reads that decoded ~390 columns for a handful. `audit` read the
competed, extracted and scored artifacts unprojected for one to three columns.
`FeatureSchema::read`'s reconstruction fallback (added earlier this session)
decoded every column of the artifact and then used only the column NAMES -- add
`table::column_names`, the footer-only counterpart to `table::nrows`, and use it.
The fragment library read now projects too: `ion_type`, `ordinal`, `frag_charge`
and `cardinality` are never read by `Library`, but were decoded and copied at
hundreds of millions of rows, one of them a string column.

`run`'s manifest still stamped the bare config hash on the four convert
artifacts, at literally the line the audit item cited. The cap-folded
`convert_hash` reached the artifact reports but not the manifest, so two runs
differing only in `--top-peaks-ms2` recorded identical provenance for their
spectra there, disagreeing with the report written beside the same file.

Extract verified byte-identical after the library projection (blake2b
7bf04c61b856be83 / e6dfae90e97f7995).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`Fragment` carried a `name: String` that is a pure function of the `ion_type`,
`ordinal` and `charge` it already holds, and it was eagerly `format!`ed for every
theoretical ion of every peptidoform. Top-N truncation then discards the large
majority of those ions, so most of the allocation existed only to be dropped --
the single most cross-corroborated finding in the audit (three independent agents:
index, mass, predict_frag).

Replace the field with a `name()` method. The only consumer that needs the text is
the library writer, which now materialises it per row it actually keeps.

Verified byte-identical: fragment library precursors and fragments (blake2b
4ce17f467df034a2 / e5295ffaf99bb905) over 1,792 candidates / 10,752 fragments.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…p by value

The audit asked for the misleading config to be wired or deleted. Deleting is not
available: `deny_unknown_fields` means removing a field breaks every config that
sets it, and these are the recorded design values for a planned tier. So state the
truth in the schema itself, where a config author reads it.

`MbrConfig::rt_window_s`, `decoy_transfer` and `requant_all` are now documented as
NOT YET WIRED, with the reason the default has the value it does. `MbrStrategy`'s
doc says plainly that only `None` vs not-`None` is implemented and that selecting
`RtTransfer` or `Full` today behaves exactly like `EmpiricalLibrary`, instead of
describing a three-stage ladder that does not exist. `validate()` already warns in
both cases.

The `mode=unique_evidence` no-op warning was keyed on which column the evidence
estimate came from, and so never fired in production: `peak_contested_frac` is part
of the Extended feature set unconditionally, so the estimate always reported itself
as "contested-discounted" -- even when that column is all zeros because the
competition features are off, which discounts nothing and leaves the raw matched
count, exactly the no-op the warning exists to announce. Test the VALUES instead:
warn when no candidate's evidence falls below the threshold, so nothing can be
removed. That also covers cases the name-based test never could.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…the real cost

The fine-tune sidecar is now the largest single cost in the pipeline (~40 min per
file, against 212 s for the whole extract stage), and almost none of it is the
fine-tuning: training runs on a few thousand confident seed peptides, while
predicting a new iRT for every peptidoform in the library dominates.

Three changes, none altering a predicted value:

`is_std` ran `re.sub` on every library row to strip modifications before checking
the residue alphabet. It now only pays the regex when the peptidoform actually
contains a bracketed modification, which most rows do not. Verified equivalent to
the original predicate on 40,021 peptidoform strings including decoy prefixes,
multiple and adjacent modifications, terminal mods, empty strings, non-standard
residues (U/X/B/Z), underscores and lowercase.

`--predict-threads` (default 0 = reuse `--threads`, i.e. unchanged) separates the
thread count for the prediction phase from the training phase. The documented crash
this file guards against was OpenMP oversubscription during fine-tuning's sustained
BACKWARD pass; prediction is forward-only and is the phase that actually dominates,
so it can usually take more threads. Left at the old value by default because that
claim is untested here -- it needs a measured trial.

The prediction loop now logs per-chunk seconds, throughput and an ETA, plus a phase
total. The audit's measurement plan asks for the fine-tune to be attributable and it
previously printed only a bare progress counter, so there was no way to tell whether
a thread change helped.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
… per run

`run_experiment` fine-tuned DeepLC inside `process_run`, so with
`rt_im_train.finetune_deeplc` set it paid one full fine-tune PER RUN. That is the
most expensive step in the pipeline by a wide margin -- roughly 40 min per run
against 212 s for the whole extract stage -- so an 80-run experiment spent about
53 hours fine-tuning and 4.7 hours searching.

Add `experiment.finetune_scope`, defaulting to `first_run_only`: the first run is
processed alone, its fine-tuned precursor library is handed to every other run, and
each run still fits its OWN retention-time calibration against it. That per-run fit
is LOESS by default (`CalibrationMethod::Loess`) and already existed -- it is what
absorbs run-to-run chromatographic drift, which is why per-run weight adaptation
buys little once every run is calibrated separately. `per_run` restores the old
behaviour.

Cost on 80 runs: one fine-tune instead of eighty. Combined with
`experiment.parallel_runs`, the remaining per-run chains still fan out; only the
first run is forced to be sequential, because the others depend on its output.

Verified end-to-end on a two-run experiment with the fine-tune enabled: run 1 runs
the sidecar and reports `status="loess"`, run 2 logs no sidecar at all, reuses the
shared library and fits its own LOESS. If the first run somehow produces no library
the others fall back to fine-tuning themselves, so the path is never silently wrong.

Caveat recorded for whoever validates this on real data: the two runs in that test
are disjoint mid-gradient slices of one file, i.e. the WORST case for reuse, since
run 2's peptides elute outside run 1's training range. Its calibrated RT window
widened accordingly (58 s vs 167 s). On a batch of real runs sharing one method the
drift is far smaller, but the quality side of this tradeoff should be confirmed on
two or three real files before committing an 80-run batch to it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`FinetuneScope::FirstRunOnly` claimed the per-run LOESS calibration is "what
absorbs run-to-run chromatographic drift", implying reuse is close to free. It is
not, and the claim was reasoning rather than measurement. Six-run ProteoBench HYE
AIF, 2026-07-28:

  A_01 (owns the fine-tune)  median |RT residual| 15.2 s, window 145 s
  A_02 A_03 B_01 B_02 B_03   20.5  20.3  20.9  24.9  25.4 s, windows 179-227 s

+7.2 s and +47% on average, and MONOTONIC in acquisition order -- real drift that
one fine-tune cannot track. Per-run LOESS corrects the slope (0.96-0.99) but not the
scatter. The wider windows also cost compute: extract roughly doubled (126 -> 203-242
s) and features up to tripled (116 -> 215-388 s), clawing back part of the saving,
so the honest figure for an 80-run batch is ~48 h saved against ~6.5 h added, not the
clean ~53 h the original message implied.

Left as the default, because one 36.5 min fine-tune instead of N still dominates that
arithmetic, but the doc now says plainly that `PerRun` is preferable when the hours
are affordable, and that periodic re-fine-tuning is the better design for very large
batches, since the drift grows with batch length.

Identification impact is not yet known; the experiment-wide rescore is still running.
No behaviour change.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The experiment-wide rescore was 88.5% of a 6-run HYE experiment (671.6 of 759 min).
The cause was not the model and not the training volume:

  * the MLP is 387-128-64-1, 57,985 parameters; a step needs ~29 us of arithmetic,
    and shrinking it measured 1.02x
  * `run_pin_sidecar` serialised the feature matrix as a 30.18 GB tab-separated
    PIN. That exceeded the worker's 4 GB streaming threshold, so it fell back to a
    disk-backed memmap and then re-read 12.77 GB on EVERY self-training iteration
    -- roughly 385 GB of reads across 30 iterations. Measured 4.55 ms/PSM
    streaming against 3.63 ms/PSM in memory.

Add `rescore.handoff = tsv | parquet`. Under `parquet` the feature table is written
as Parquet with f32 features, streamed a batch at a time through a new
`mumdia_io::table::BatchWriter`, so nothing beyond one block is materialised --
`feats` is already resident as Vec<Vec<f64>> (~27 GB on an experiment-wide pool).

f32, not f64, deliberately: the TSV wrote `{:.6}` and the worker casts to float32
regardless, so f32 is what is actually consumed. It halves the file rather than
silently raising precision above the reference the result was validated against.

Parquet is nn_torch only. `mokapot_worker.py` goes through `mokapot.read_pin()`,
which requires the tab-separated form, so a mokapot run warns and writes the PIN
instead of failing a configured job.

The worker reads either format, chosen by extension, with its own low-memory path:
feature columns are read batch-wise straight into one preallocated float32 array and
standardised in place. Routing Parquet through `read_table().to_pandas()` first held
~90 of 95.7 GB and had to be killed; batched it peaks near 34 GB. It standardises by
mean/std, matching the streaming backend the reference run used, so the comparison
is not confounded by a different standardisation.

Measured on the same 8,858,206-PSM pool, scored by an identical metric:

  original (CPU, TSV, streaming, cold, 10 iter)  671.6 min  74,430 pep  0.988% decoy
  parquet + in-memory + GPU + warm, 10 iter       12   min  74,048 pep  0.988%
  parquet + in-memory + GPU + warm, 20 iter       26   min  74,339 pep  0.988%
  parquet + in-memory + GPU + warm,  5 iter        8   min  72,465 pep  0.988%

56x at 10 iterations. FDR is unaffected -- 0.988% in every configuration. The
init-feature scan reproduces EXACTLY (per-fold counts 8103 / 7918 / 8179, identical
to the original), which also demonstrates the Parquet path yields the same feature
values as the TSV path.

Default stays `tsv`: the numbers above want confirming on a second dataset before a
default changes, and `parquet` currently benefits only the nn path.

Also in the worker, all default-off: `MUMDIA_NN_DEVICE` (auto|cuda|cpu; errors
rather than silently using CPU on a CPU-only torch build -- which is why an RTX 4090
sat idle through the 11-hour run), `WARM_START`/`WARM_EPOCHS` (carry weights and Adam
state forward; 1.65x on training for -0.2 to -0.5% peptides), `TRAIN_SUB`,
`NEG_RATIO`, `FEATURES` (applied before the file is read, so dropped columns are
never parsed), `INIT_TOPK` (7.2x on the init scan but NOT exact -- tie ordering at the
window edge shifted 30 of 774 counts by 1-2), and per-phase wall-clock timers, which
are what located the bottleneck in the first place.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A portability audit of the deployment target (Linux, no GPU, 83 mzML files) found three
defects that could each waste the whole batch, plus operator traps worth writing down.

Fixed:

BatchNorm1d raises "Expected more than 1 value per channel when training" on a
trailing minibatch of exactly one row, and nothing guarded it. The training-set size
changes with the selected positive set every iteration, so this was a latent crash that
could surface at an arbitrary point -- including the final experiment-wide rescore,
after every run's compute is already spent. Reproduced the ValueError directly. The
single row is now dropped for that epoch; the permutation reshuffles next epoch, so no
row is systematically excluded.

MUMDIA_NN_THREADS was silently overridden by any site-wide OMP_NUM_THREADS.
OMP_NUM_THREADS=1 is a common cluster module default, which would have pinned an 83-run
rescore to one torch thread with no way to override it from MuMDIA's own knob.
Precedence is now explicit-MuMDIA-knob-wins, OMP honoured only when the MuMDIA one is
unset, and the worker prints which source it used.

`run_experiment` preflight checked only that interpreter fields were SET, not that they
exist, while the rescore runs after all per-run chains -- so a mistyped `rescore.python`
discarded days of compute at the last stage. It now fails fast and points at `mumdia
doctor`. The sidecar spawn also reported a bare "No such file or directory" naming
nothing; it now names the interpreter and script.

Added `config.strasbourg-linux.json` and `docs/21_linux_deployment.md` covering the
things the existing docs get wrong for this target (`19_getting_started.md` is entirely
Windows, and `17_troubleshooting.md` recommends a native Windows path):

  * a Windows-style `c:/...` path is RELATIVE on Linux, so script resolution builds
    nonsense paths before failing with a message naming a Windows directory
  * `.cargo/config.toml` is gitignored but travels with rsync, and its Windows
    `target-dir` makes cargo build into `rust/mumdia/C:/Users/...` where nothing looks
  * `MUMDIA_NN_FOLDS/ITERS/TRAIN_FDR` in the shell are silently overridden by the Rust
    caller -- the config fields are the real control
  * `DEEPLC_FT_THREADS` is the ONLY handle on the DeepLC phase (the caller never passes
    `--threads`/`--predict-threads`) and defaults to 8 regardless of core count
  * `--out-dir` must be local disk: the streaming memmap lands in `sidecar_work/`
  * `sidecar_work/` is never cleaned up
  * memory as ~28 B/fragment plus hit accumulation, with `parallel_runs` guidance by RAM,
    noting the fragment index is SHARED across concurrent runs rather than duplicated

The config sets `handoff: parquet` (671.6 min -> 12 min on the 6-run HYE pool, decoy
fraction unchanged at 0.988%), `finetune_scope: first_run_only` with per-run LOESS as
requested, and `parallel_runs: 1` as a deliberately conservative start.

Flagged in the doc rather than assumed: the Parquet handoff and BatchWriter have only
ever run on Windows, so section 5 is a two-file smoke test that exercises them before
the batch, with `handoff: tsv` as the one-field fallback.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Added when config validation started emitting warnings for inert MBR fields. The
Linux deployment runbook builds with `--locked`, so the lockfile has to match the
manifest or the build fails on a clean checkout.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@RobbinBouwmeester
RobbinBouwmeester merged commit c236a87 into main Jul 30, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant