diff --git a/rust/apex_sim/examples/bench.rs b/rust/apex_sim/examples/bench.rs index 284f6cb6..cab39c8e 100644 --- a/rust/apex_sim/examples/bench.rs +++ b/rust/apex_sim/examples/bench.rs @@ -1,42 +1,50 @@ -//! Apex-finder bench: runs the whole catalog of canonical scenarios -//! (`apex_sim::bench::canonical_suite`) at once and prints per-scenario detail -//! plus a comparison summary table. +//! Apex-finder + score benches. //! -//! Each scenario runs `n_runs` seeds of one fixed config; "sensitivity" is the -//! fraction of runs whose pass-2 apex lands within `tol` cycles of truth. +//! Runs three suites: +//! * canonical (historical 245-cyc scenarios; apex-recovery only), +//! * broad (production Phase-1 window ~1695 cyc; apex-recovery), and +//! * narrow (production Phase-3 window ~150 cyc; apex-recovery + score AUC). +//! "recovery" = fraction of runs whose pass-2 apex lands within `tol` cycles of +//! the (jittered) truth. "AUC" = ROC-AUC separating the Pass-2 score of a real +//! peptide from a pure-noise twin (signal-vs-noise discrimination). //! //! Run: //! cargo run -p apex_sim --release --example bench -//! cargo run -p apex_sim --release --example bench -- 500 2 -//! (positional: ) +//! cargo run -p apex_sim --release --example bench -- 2000 2 500 +//! (positional: ) use apex_sim::bench::{ + SensitivityReport, + broad_suite, canonical_suite, + narrow_suite, + run_discrimination, run_sensitivity, }; +use apex_sim::sim::SimParams; -fn main() { - let args: Vec = std::env::args().collect(); - let n_runs: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(2000); - let tol: i64 = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(2); - - let suite = canonical_suite(); - let mut rows = Vec::with_capacity(suite.len()); - for (name, cfg) in &suite { - let report = run_sensitivity(cfg, n_runs, tol); - report.print(name); - println!(); - rows.push((*name, report)); - } +/// Run a suite and print a comparison summary table under `=== {title} ===`. +/// The canonical suite passes `title = "summary"` so `bench_out/score.sh` +/// (which keys on `=== summary`) keeps working. +fn run_and_print_recovery( + suite: &[(&'static str, SimParams)], + n_runs: usize, + tol: i64, + title: &str, +) { + let rows: Vec<(&str, SensitivityReport)> = suite + .iter() + .map(|(name, cfg)| (*name, run_sensitivity(cfg, n_runs, tol))) + .collect(); - println!("=== summary (n={n_runs}, tol=±{tol} cycles) ==="); + println!("=== {title} (n={n_runs}, tol=±{tol} cycles) ==="); println!( - " {:<26} {:>7} {:>7} {:>8} {:>9}", + " {:<28} {:>7} {:>7} {:>8} {:>9}", "scenario", "pass2%", "pass1%", "medErr", "us/run" ); for (name, r) in &rows { println!( - " {:<26} {:>7.1} {:>7.1} {:>8} {:>9.2}", + " {:<28} {:>7.1} {:>7.1} {:>8} {:>9.2}", name, r.pass2_pct(), r.pass1_pct(), @@ -44,4 +52,36 @@ fn main() { r.score_us_per_run(), ); } + println!(); +} + +fn main() { + let args: Vec = std::env::args().collect(); + let n_runs: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(2000); + let tol: i64 = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(2); + // Broad sims are ~42x heavier per run; default to fewer. + let broad_runs: usize = args.get(3).and_then(|s| s.parse().ok()).unwrap_or(500); + + // Canonical (historical) apex-recovery suite. Header stays `=== summary`. + run_and_print_recovery(&canonical_suite(), n_runs, tol, "summary"); + + // Broad apex-finding (production Phase-1 window ~1695 cyc). + run_and_print_recovery(&broad_suite(), broad_runs, tol, "broad apex-finding"); + + // Narrow scoring (production Phase-3 window ~150 cyc): apex recovery. + run_and_print_recovery(&narrow_suite(), n_runs, tol, "narrow recovery"); + + // Narrow score discrimination: signal-present vs pure-noise (ROC-AUC). + println!("=== narrow score discrimination (AUC, n_pairs={n_runs}) ==="); + println!( + " {:<28} {:>7} {:>12} {:>12}", + "scenario", "AUC", "med+signal", "med-noise" + ); + for (name, cfg) in &narrow_suite() { + let d = run_discrimination(cfg, n_runs); + println!( + " {:<28} {:>7.3} {:>12.3e} {:>12.3e}", + name, d.auc, d.median_present, d.median_absent + ); + } } diff --git a/rust/apex_sim/src/bench.rs b/rust/apex_sim/src/bench.rs index 8237e424..8b1dda64 100644 --- a/rust/apex_sim/src/bench.rs +++ b/rust/apex_sim/src/bench.rs @@ -94,6 +94,13 @@ pub fn canonical_suite() -> Vec<(&'static str, SimParams)> { absent.random_peaks.count = 100; absent.real_fragments[0].obs_scale = 0.0; + // Absent precursor at stress level: MS1 precursor undetected (intensity 0), + // fragments intact. Guards against over-reliance on precursor signal. + let mut absent_prec = base(); + absent_prec.noise_floor = 0.5; + absent_prec.random_peaks.count = 100; + absent_prec.precursor_intensity = 0.0; + vec![ ("clean", clean), ("moderate_noise", moderate), @@ -101,6 +108,101 @@ pub fn canonical_suite() -> Vec<(&'static str, SimParams)> { ("heavy_interference", heavy), ("mismatched_library", mismatched), ("absent_top_fragment", absent), + ("absent_precursor", absent_prec), + ] +} + +/// Production-realistic BROAD apex-finding suite (Phase-1 window ~1695 cycles, +/// RT-unrestricted). Interferents are DENSITY-scaled (fixed per-cycle rate) and +/// the apex is jittered per seed with a sub-cycle offset, so widening the window +/// does not make the task artificially easy and the finder faces a moving, +/// off-grid target. Each scenario is its own line so it can be commented out. +pub fn broad_suite() -> Vec<(&'static str, SimParams)> { + let base = || { + let mut p = SimParams::default(); + p.n_cycles = 1695; + p.apex_cycle = 847.0; + p.apex_jitter = Some(800.0); + p.width_sigma = 1.0; + p.seed = 680; + p.random_peaks.count = 0; + p.random_peaks.density_per_cycle = 0.41; // matches canonical 100/245 + p + }; + + let mut clean = base(); + clean.noise_floor = 0.02; + clean.random_peaks.enabled = false; + + let mut moderate = base(); + moderate.noise_floor = 0.2; + + let mut stress = base(); + stress.noise_floor = 0.5; + + let mut hard = base(); + hard.noise_floor = 0.5; + hard.random_peaks.hardness = 3.0; // 3x interferent density + + let mut mismatched = base(); + mismatched.noise_floor = 0.5; + for f in &mut mismatched.real_fragments { + f.obs_scale = f.theo_intensity; + f.theo_intensity = 1.0; + } + + vec![ + ("broad_clean", clean), + ("broad_moderate_noise", moderate), + ("broad_high_noise+interf", stress), + ("broad_hard_3x_density", hard), + ("broad_mismatched_library", mismatched), + ] +} + +/// Production-realistic NARROW scoring suite (Phase-3 calibrated window ~150 +/// cycles). Same density-scaled interferents + jittered sub-cycle apex; used +/// for BOTH apex recovery and score discrimination (`run_discrimination`). +pub fn narrow_suite() -> Vec<(&'static str, SimParams)> { + let base = || { + let mut p = SimParams::default(); + p.n_cycles = 150; + p.apex_cycle = 75.0; + p.apex_jitter = Some(60.0); + p.width_sigma = 1.0; + p.seed = 680; + p.random_peaks.count = 0; + p.random_peaks.density_per_cycle = 0.41; + p + }; + + let mut clean = base(); + clean.noise_floor = 0.02; + clean.random_peaks.enabled = false; + + let mut moderate = base(); + moderate.noise_floor = 0.2; + + let mut stress = base(); + stress.noise_floor = 0.5; + + let mut hard = base(); + hard.noise_floor = 0.5; + hard.random_peaks.hardness = 3.0; + + let mut mismatched = base(); + mismatched.noise_floor = 0.5; + for f in &mut mismatched.real_fragments { + f.obs_scale = f.theo_intensity; + f.theo_intensity = 1.0; + } + + vec![ + ("narrow_clean", clean), + ("narrow_moderate_noise", moderate), + ("narrow_high_noise+interf", stress), + ("narrow_hard_3x_density", hard), + ("narrow_mismatched_library", mismatched), ] } @@ -108,7 +210,9 @@ pub fn canonical_suite() -> Vec<(&'static str, SimParams)> { /// collect apex-recovery + timing stats. A hit = pass-2 apex within `tol` /// cycles of the configured `apex_cycle`. pub fn run_sensitivity(base: &SimParams, n_runs: usize, tol: i64) -> SensitivityReport { - let true_apex = base.apex_cycle.round() as i64; + // Echoed in the report header; each run is scored against its OWN realized + // apex (jitter varies it per seed), not this nominal value. + let nominal_apex = base.apex_cycle.round() as i64; // rt mapping is seed-independent, so one mapper serves every run. let map = base.rt_mapper(); @@ -131,6 +235,7 @@ pub fn run_sensitivity(base: &SimParams, n_runs: usize, tol: i64) -> Sensitivity let mut scorer = TraceScorer::new(base.n_cycles, base.real_fragments.len().max(1)); let t_score = Instant::now(); for data in &sims { + let true_apex = data.realized_apex_cycle.round() as i64; match scorer::run_with(&mut scorer, &data.extraction, &map) { Ok((pass1, pass2)) => { if (pass1.apex_cycle as i64 - true_apex).abs() <= tol { @@ -151,7 +256,7 @@ pub fn run_sensitivity(base: &SimParams, n_runs: usize, tol: i64) -> Sensitivity SensitivityReport { n: n_runs, tol, - true_apex, + true_apex: nominal_apex, build_ms, score_ms, errors, @@ -246,3 +351,118 @@ impl SensitivityReport { println!(); } } + +/// Present-vs-absent score-discrimination result for one scenario. +pub struct DiscriminationReport { + pub n_pairs: usize, + /// ROC-AUC = P(score_present > score_absent), 0.5 tie credit. 1.0 = perfect + /// signal/noise separation, 0.5 = useless. + pub auc: f64, + pub median_present: f32, + pub median_absent: f32, +} + +/// ROC-AUC of `present` vs `absent` score populations via the average-rank +/// Mann-Whitney statistic: `P(present > absent)` with 0.5 credit for ties. +pub fn roc_auc(present: &[f32], absent: &[f32]) -> f64 { + let (np, na) = (present.len(), absent.len()); + if np == 0 || na == 0 { + return f64::NAN; + } + // Combined values tagged present(true)/absent(false), sorted ascending. + let mut all: Vec<(f32, bool)> = present + .iter() + .map(|&v| (v, true)) + .chain(absent.iter().map(|&v| (v, false))) + .collect(); + all.sort_by(|a, b| a.0.total_cmp(&b.0)); + // Sum of average ranks (1-based) assigned to present values. + let mut rank_sum_present = 0.0f64; + let mut i = 0usize; + while i < all.len() { + let mut j = i + 1; + while j < all.len() && all[j].0 == all[i].0 { + j += 1; + } + // Mean rank of the tie block covering positions (i+1..=j). + let avg_rank = ((i + 1 + j) as f64) / 2.0; + for entry in &all[i..j] { + if entry.1 { + rank_sum_present += avg_rank; + } + } + i = j; + } + let u = rank_sum_present - (np * (np + 1)) as f64 / 2.0; + u / (np as f64 * na as f64) +} + +/// Score `n_pairs` matched present/absent realizations and report how well the +/// production Pass-2 score separates real signal from pure noise. Each pair +/// shares a seed, so the "absent" twin (all real fragments `obs_scale = 0`) has +/// IDENTICAL noise + interferents — only the real peak differs. +pub fn run_discrimination(base: &SimParams, n_pairs: usize) -> DiscriminationReport { + let map = base.rt_mapper(); + let mut scorer = TraceScorer::new(base.n_cycles, base.real_fragments.len().max(1)); + let mut present: Vec = Vec::with_capacity(n_pairs); + let mut absent: Vec = Vec::with_capacity(n_pairs); + for i in 0..n_pairs { + let mut pp = base.clone(); + pp.seed = base.seed.wrapping_add(i as u64); + // Pure-noise twin: NO real signal at all (fragments AND precursor + // zeroed), identical seed => identical noise + interferents. Zeroing + // only fragments would leave the precursor isotope peaks in place and + // the score would still see real signal. + let mut ap = pp.clone(); + for f in &mut ap.real_fragments { + f.obs_scale = 0.0; + } + ap.precursor_intensity = 0.0; + if let Ok((_, s)) = scorer::run_with(&mut scorer, &sim::build(&pp).extraction, &map) { + present.push(s.score); + } + if let Ok((_, s)) = scorer::run_with(&mut scorer, &sim::build(&ap).extraction, &map) { + absent.push(s.score); + } + } + let median = |mut v: Vec| { + v.sort_by(f32::total_cmp); + v.get(v.len() / 2).copied().unwrap_or(f32::NAN) + }; + let auc = roc_auc(&present, &absent); + DiscriminationReport { + n_pairs, + auc, + median_present: median(present), + median_absent: median(absent), + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn roc_auc_basic() { + assert!((roc_auc(&[3.0, 4.0, 5.0], &[0.0, 1.0, 2.0]) - 1.0).abs() < 1e-9); + assert!((roc_auc(&[0.0, 1.0], &[2.0, 3.0]) - 0.0).abs() < 1e-9); + assert!((roc_auc(&[1.0, 1.0], &[1.0, 1.0]) - 0.5).abs() < 1e-9); + } + + #[test] + fn discrimination_high_when_signal_present() { + // Genuinely clean: low noise, no injected interferents. A real peptide + // vs pure noise must separate near-perfectly. + let mut base = SimParams::default(); + base.n_cycles = 150; + base.noise_floor = 0.05; + base.random_peaks.enabled = false; + let rep = run_discrimination(&base, 200); + assert!( + rep.auc > 0.9, + "clean signal should separate: auc={}", + rep.auc + ); + assert!(rep.median_present > rep.median_absent); + } +} diff --git a/rust/apex_sim/src/sim.rs b/rust/apex_sim/src/sim.rs index 94eac3ab..324c71e1 100644 --- a/rust/apex_sim/src/sim.rs +++ b/rust/apex_sim/src/sim.rs @@ -67,6 +67,14 @@ pub struct RandomPeaks { pub height_frac_max: f32, /// If true, precursor rows are also eligible targets. pub hit_precursors: bool, + /// Interferents per cycle; when > 0 this overrides `count` as + /// `round(density_per_cycle * n_cycles * hardness)`, so interferent DENSITY + /// (not absolute count) stays fixed as the window widens — a fixed count in + /// a wide window is an artificially easy task. + pub density_per_cycle: f32, + /// Multiplier (>=1) to deliberately over-load interferents above realistic + /// density for stress-testing. + pub hardness: f32, } impl Default for RandomPeaks { @@ -77,16 +85,33 @@ impl Default for RandomPeaks { height_frac_min: 0.1, height_frac_max: 10.0, hit_precursors: true, + density_per_cycle: 0.0, + hardness: 1.0, } } } +/// Effective interferent count for a window of `n_cycles`: density-scaled when +/// `density_per_cycle > 0`, else the absolute `count`. +pub fn effective_random_count(rp: &RandomPeaks, n_cycles: usize) -> usize { + if rp.density_per_cycle > 0.0 { + (rp.density_per_cycle * n_cycles as f32 * rp.hardness.max(1.0)).round() as usize + } else { + rp.count + } +} + /// Full set of simulation knobs. Deterministic given `seed`. #[derive(Debug, Clone)] pub struct SimParams { pub n_cycles: usize, /// Shared apex position (cycle index) for all real fragments + precursor. pub apex_cycle: f32, + /// When `Some(range)`, each run draws the realized apex from the seeded rng + /// as `apex_cycle + U(-range, range) + U(-0.5, 0.5)` (position varies per + /// seed AND lands off-grid). `None` pins it exactly to `apex_cycle` (the + /// historical, unrealistically-easy on-grid, fixed-position behavior). + pub apex_jitter: Option, /// Elution peak width (gaussian sigma, in cycles). pub width_sigma: f32, /// Global height scale applied on top of each fragment's rel_intensity. @@ -134,6 +159,7 @@ impl Default for SimParams { Self { n_cycles: 60, apex_cycle: 30.0, + apex_jitter: None, width_sigma: 1.0, height: 1000.0, real_fragments, @@ -172,6 +198,9 @@ pub struct SimData { pub extraction: Extraction, pub fragment_rows: Vec, pub precursor_rows: Vec, + /// The apex cycle actually used to generate this realization (equals + /// `apex_cycle` when `apex_jitter` is `None`, else the jittered value). + pub realized_apex_cycle: f32, } /// Gaussian elution value at `cycle` for a peak centered at `center`. @@ -186,6 +215,19 @@ pub fn build(params: &SimParams) -> SimData { let n = params.n_cycles; let dummy_mz = 500.0_f64; + // Realized apex: jittered per-seed + sub-cycle when enabled, else pinned. + // Drawn BEFORE any signal/noise sampling; the `None` branch consumes no rng + // so historical scenarios stay byte-identical. + let realized_apex = match params.apex_jitter { + Some(range) => { + let lo = 3.0f32; + let hi = (n as f32 - 4.0).max(lo); + (params.apex_cycle + rng.random_range(-range..=range) + rng.random_range(-0.5..=0.5)) + .clamp(lo, hi) + } + None => params.apex_cycle, + }; + // --- Expected intensities: THEORETICAL (library) ratios. --- // These drive cosine/scribe. Observed peaks below may deviate (obs_scale). let expected = ExpectedIntensities::try_from_pairs( @@ -211,7 +253,7 @@ pub fn build(params: &SimParams) -> SimData { let noise = params.noise_floor * params.height * f.noise_mult; let intensities = (0..n) .map(|c| { - let signal = gaussian(c as f32, params.apex_cycle, params.width_sigma, peak); + let signal = gaussian(c as f32, realized_apex, params.width_sigma, peak); sample_cell(&mut rng, signal, noise) }) .collect(); @@ -231,7 +273,7 @@ pub fn build(params: &SimParams) -> SimData { let noise = params.noise_floor * params.height; let intensities = (0..n) .map(|c| { - let signal = gaussian(c as f32, params.apex_cycle, params.width_sigma, peak); + let signal = gaussian(c as f32, realized_apex, params.width_sigma, peak); sample_cell(&mut rng, signal, noise) }) .collect(); @@ -271,7 +313,7 @@ pub fn build(params: &SimParams) -> SimData { let chromatograms = ChromatogramCollector:: { id: 0, mobility_ook0: 1.0, - rt_seconds: (map((params.apex_cycle as usize).min(n - 1)) as f32) / 1000.0, + rt_seconds: (map((realized_apex as usize).min(n - 1)) as f32) / 1000.0, precursor_mono_mz: dummy_mz, precursor_charge: 2, precursor_mz_limits: (dummy_mz - 1.0, dummy_mz + 1.0), @@ -310,6 +352,7 @@ pub fn build(params: &SimParams) -> SimData { extraction, fragment_rows: frag_rows, precursor_rows: prec_rows, + realized_apex_cycle: realized_apex, } } @@ -323,10 +366,11 @@ fn inject_random_peaks( params: &SimParams, ) { let rp = ¶ms.random_peaks; - if !rp.enabled || rp.count == 0 { + let n = params.n_cycles; + let count = effective_random_count(rp, n); + if !rp.enabled || count == 0 { return; } - let n = params.n_cycles; let n_frag = frag_rows.len(); let n_targets = n_frag + if rp.hit_precursors { @@ -338,7 +382,7 @@ fn inject_random_peaks( return; } - for _ in 0..rp.count { + for _ in 0..count { let target = rng.random_range(0..n_targets); let row = if target < n_frag { &mut frag_rows[target] @@ -382,6 +426,44 @@ fn fill_array( mod tests { use super::*; + #[test] + fn density_scales_interferent_count_with_window() { + let mut p = SimParams::default(); + p.random_peaks.count = 0; + p.random_peaks.density_per_cycle = 2.0; + p.random_peaks.hardness = 1.0; + assert_eq!(effective_random_count(&p.random_peaks, 100), 200); + assert_eq!(effective_random_count(&p.random_peaks, 1000), 2000); + p.random_peaks.hardness = 3.0; + assert_eq!(effective_random_count(&p.random_peaks, 100), 600); + // count path preserved when density is 0 + let mut q = RandomPeaks::default(); + q.density_per_cycle = 0.0; + q.count = 77; + assert_eq!(effective_random_count(&q, 5000), 77); + } + + #[test] + fn apex_jitter_is_deterministic_and_subcycle() { + let mut p = SimParams::default(); + p.n_cycles = 200; + p.apex_cycle = 100.0; + p.apex_jitter = Some(20.0); + let a = build(&p).realized_apex_cycle; + let b = build(&p).realized_apex_cycle; // same seed => identical + assert_eq!(a, b); + assert!((a - 100.0).abs() <= 20.5 + 1e-3); + // vary seed => different position + let mut q = p.clone(); + q.seed = p.seed.wrapping_add(1); + let c = build(&q).realized_apex_cycle; + assert!((a - c).abs() > 1e-6); + // no jitter => exact + let mut r = p.clone(); + r.apex_jitter = None; + assert_eq!(build(&r).realized_apex_cycle, 100.0); + } + #[test] fn clean_peak_apex_lands_at_configured_cycle() { let mut p = SimParams::default(); diff --git a/rust/timsseek/src/scoring/apex/mod.rs b/rust/timsseek/src/scoring/apex/mod.rs new file mode 100644 index 00000000..10ae4f10 --- /dev/null +++ b/rust/timsseek/src/scoring/apex/mod.rs @@ -0,0 +1,32 @@ +//! Tunable knobs for the composite apex profile. +//! +//! The apex profile (`cos^cos_pow * I^i_exp * (s_ratio + s_norm)`, optionally +//! blurred) is built in [`crate::scoring::apex_finding`]. This module holds the +//! bench-validated shipping values for its knobs. + +/// Tunable apex-profile knobs. `Default` holds the bench-validated shipping +/// values. +#[derive(Debug, Clone)] +pub struct ApexConfig { + /// Cosine exponent in the base profile: `cos^cos_pow`. + pub cos_pow: f32, + /// Log-intensity exponent in the base profile: `I^i_exp`. + pub i_exp: f32, + /// Additive scribe floor: `apex = C * (s_ratio + s_norm)`. + pub s_ratio: f32, + /// Gaussian-blur passes applied to the base profile (0 = none). + pub blur_passes: usize, +} + +impl Default for ApexConfig { + fn default() -> Self { + // Optuna-tuned on the apex_sim canonical+broad+narrow recovery suites + // (summed pass2%), rounded to robust non-degenerate values. + Self { + cos_pow: 0.25, + i_exp: 1.0, + s_ratio: 0.25, + blur_passes: 2, + } + } +} diff --git a/rust/timsseek/src/scoring/apex_finding.rs b/rust/timsseek/src/scoring/apex_finding.rs index 4706d954..b5dd2898 100644 --- a/rust/timsseek/src/scoring/apex_finding.rs +++ b/rust/timsseek/src/scoring/apex_finding.rs @@ -1,8 +1,7 @@ //! Logic for finding the apex of a peptide elution profile. //! -//! This module replaces the old `LocalizationBuffer` and `calculate_scores.rs`. -//! It implements an efficient, accumulator-based scoring engine that avoids -//! intermediate data transposition. +//! An efficient, accumulator-based scoring engine that avoids intermediate +//! data transposition. //! //! # Usage //! @@ -27,6 +26,7 @@ use std::fmt::Display; +use super::apex::ApexConfig; use super::{ NUM_MS1_IONS, NUM_MS2_IONS, @@ -42,7 +42,6 @@ use crate::scoring::scores::apex_features::{ compute_apex_features, compute_split_product, compute_weighted_score, - find_joint_apex, }; use crate::utils::top_n_array::TopNArray; use serde::Serialize; @@ -288,6 +287,7 @@ pub struct TraceScorer { coel_scratch: crate::scoring::scores::apex_features::CoelutionScratch, cosine_profile: Vec, scribe_profile: Vec, + cfg: ApexConfig, } #[derive(Debug)] @@ -343,6 +343,7 @@ impl TraceScorer { ), cosine_profile: Vec::with_capacity(capacity), scribe_profile: Vec::with_capacity(capacity), + cfg: ApexConfig::default(), } } @@ -465,13 +466,10 @@ impl TraceScorer { let delta_next = cycle_val - next_val; let delta_second_next = cycle_val - second_next_val; - // Joint apex: find precursor-fragment agreement, but use clicked cycle if far from it - let joint_apex = find_joint_apex(&self.cosine_profile, &self.traces.ms1_precursor_trace); - let effective_apex = if (joint_apex as i64 - cycle as i64).unsigned_abs() as usize <= 3 { - joint_apex - } else { - cycle - }; + // Unified apex source: Pass 2 scores at Pass 1's weighted-apex_profile + // location. `cycle` is `suggested.apex_cycle` (the weighted argmax) — the + // single, best-validated apex finder. No re-location on a weaker profile. + let effective_apex = cycle; // 11 features at effective apex let n_cycles = self.cosine_profile.len(); @@ -529,7 +527,7 @@ impl TraceScorer { }) } - /// Convenience: compute_traces + suggest_apex. Migration aid. + /// Phase-1 entry: compute_traces + suggest_apex (apex location only). #[cfg_attr( feature = "instrumentation", tracing::instrument(skip(self, scoring_ctx, rt_mapper), level = "trace") @@ -544,7 +542,7 @@ impl TraceScorer { self.suggest_apex(rt_mapper, cycle_offset) } - /// Convenience: compute_traces + suggest_apex + score_at. Migration aid. + /// Phase-3 entry: compute_traces + suggest_apex + score_at (full score). #[cfg_attr( feature = "instrumentation", tracing::instrument(skip(self, scoring_ctx, rt_mapper), level = "trace") @@ -711,8 +709,8 @@ impl TraceScorer { /// Compute the apex profile from cosine and scribe traces. /// - /// apex_profile(t) = C(t) * (0.5 + S_norm(t)) - /// where C(t) = cosine(t)^3 * I(t) + /// apex_profile(t) = C(t) * (s_ratio + S_norm(t)), then optional gaussian blur. + /// where C(t) = cosine(t)^cos_pow * I(t)^i_exp (see ApexConfig) /// S(t) = scribe(t) * I(t) /// S_norm = (S - min(S)) / (max(S) - min(S)) #[cfg_attr( @@ -740,7 +738,7 @@ impl TraceScorer { for i in 0..len { let cos = self.traces.cosine_trace[i]; let intensity = self.traces.ms2_log_intensity[i]; - let c = cos * cos * cos * intensity; // cos^3 * I + let c = cos.powf(self.cfg.cos_pow) * intensity.powf(self.cfg.i_exp); let s = self.traces.ms2_scribe[i] * intensity; let s_norm = if s_range > 0.0 { @@ -749,7 +747,13 @@ impl TraceScorer { 0.5 // Degrade to cosine-only when scribe is constant }; - self.traces.apex_profile.push(c * (0.5 + s_norm)); + self.traces + .apex_profile + .push(c * (self.cfg.s_ratio + s_norm)); + } + + for _ in 0..self.cfg.blur_passes { + gaussblur_in_place(&mut self.traces.apex_profile); } } diff --git a/rust/timsseek/src/scoring/full_results.rs b/rust/timsseek/src/scoring/full_results.rs deleted file mode 100644 index 2f8f4d0c..00000000 --- a/rust/timsseek/src/scoring/full_results.rs +++ /dev/null @@ -1,13 +0,0 @@ -use super::apex_finding::ElutionTraces; -use super::results::ScoredCandidate; -use crate::IonAnnot; -use serde::Serialize; -use timsquery::models::aggregators::ChromatogramCollector; - -#[derive(Debug, Clone, Serialize)] -pub struct ViewerResult { - pub traces: ElutionTraces, - pub longitudinal_apex_profile: Vec, - pub chromatograms: ChromatogramCollector, - pub scored: ScoredCandidate, -} diff --git a/rust/timsseek/src/scoring/mod.rs b/rust/timsseek/src/scoring/mod.rs index 21d8f7ac..095db848 100644 --- a/rust/timsseek/src/scoring/mod.rs +++ b/rust/timsseek/src/scoring/mod.rs @@ -1,7 +1,7 @@ mod accumulator; +pub mod apex; pub mod apex_finding; pub mod extraction; -pub mod full_results; mod maybe_par; pub mod offsets; pub mod parquet_writer; diff --git a/rust/timsseek/src/scoring/pipeline.rs b/rust/timsseek/src/scoring/pipeline.rs index 08b0fe15..d66c45a8 100644 --- a/rust/timsseek/src/scoring/pipeline.rs +++ b/rust/timsseek/src/scoring/pipeline.rs @@ -51,7 +51,6 @@ use super::apex_finding::{ RelativeIntensities, TraceScorer, }; -use super::full_results::ViewerResult; use super::hyperscore::single_lazyscore; use super::offsets::MzMobilityOffsets; use super::results::{ @@ -367,42 +366,6 @@ pub struct Scorer { } impl Scorer { - #[cfg_attr( - feature = "instrumentation", - tracing::instrument(skip_all, level = "trace") - )] - fn build_broad_extraction( - &self, - item: &QueryItemToScore, - ) -> Result< - ( - super::apex_finding::PeptideMetadata, - super::apex_finding::Extraction, - ), - SkipReason, - > { - let extraction = super::extraction::build_extraction( - &item.query, - item.expected_intensity.clone(), - &self.index, - &self.broad_tolerance, - Some(TOP_N_FRAGMENTS), - )?; - - let library_rt = item.query.rt_seconds(); - let metadata = super::apex_finding::PeptideMetadata { - digest: item.digest.clone(), - charge: item.query.precursor_charge(), - library_id: extraction.chromatograms.id as u32, - library_rt, - calibrated_rt_seconds: library_rt, // no calibration in broad path - ref_mobility_ook0: item.query.mobility_ook0(), - ref_precursor_mz: item.query.mono_precursor_mz(), - }; - - Ok((metadata, extraction)) - } - /// Calculates the weighted mean ion mobility across fragments and precursors. fn get_mobility(item: &SpectralCollector) -> f64 { // Calculate weighted mean mobility from fragments @@ -525,50 +488,6 @@ impl Scorer { } impl Scorer { - pub fn score_for_viewer( - &self, - item: QueryItemToScore, - calibration: &CalibrationResult, - ) -> Result { - // One-shot worker for the viewer path (not hot). - let mut worker = - ScoringWorker::new(self.num_cycles(), item.expected_intensity.fragment_len()); - - let (metadata, scoring_ctx) = self.build_broad_extraction(&item).map_err(|_| { - DataProcessingError::ExpectedNonEmptyData { - context: Some("RT out of bounds".into()), - } - })?; - - let apex_score = worker - .scorer - .find_apex(&scoring_ctx, &|idx| self.map_rt_index_to_milis(idx))?; - let spectral_tol = calibration.get_spectral_tolerance(); - let isotope_tol = calibration.get_isotope_tolerance(); - self.execute_secondary_query(&item, &apex_score, &spectral_tol, &isotope_tol, &mut worker); - let inner_collector = worker.inner_collector.as_ref().expect("set by secondary"); - let isotope_collector = worker.isotope_collector.as_ref().expect("set by secondary"); - - let nqueries = scoring_ctx.chromatograms.fragments.num_ions() as u8; - let scored = self.finalize_results( - &metadata, - nqueries, - &apex_score, - inner_collector, - isotope_collector, - )?; - - // Extract chromatograms before it's consumed - let chromatograms = scoring_ctx.chromatograms; - - Ok(ViewerResult { - traces: worker.scorer.traces.clone(), - longitudinal_apex_profile: worker.scorer.traces.apex_profile.clone(), - chromatograms, - scored, - }) - } - /// Build a chromatogram extraction using calibrated RT and per-query tolerance. /// The speclib is NOT mutated — CalibrationResult provides the RT conversion. #[cfg_attr( diff --git a/rust/timsseek/src/scoring/scores/apex_features.rs b/rust/timsseek/src/scoring/scores/apex_features.rs index 65f3ee22..4231d48c 100644 --- a/rust/timsseek/src/scoring/scores/apex_features.rs +++ b/rust/timsseek/src/scoring/scores/apex_features.rs @@ -447,39 +447,13 @@ pub fn compute_split_product( } } -/// Find the joint precursor-fragment apex (METHODS.md Section 3.3). -/// -/// joint(t) = C(t) * (0.5 + P(t) / max(P)) -/// If max(P) == 0, degrades to joint(t) = C(t) * 0.5 (pure fragment apex). -pub fn find_joint_apex(cosine_profile: &[f32], precursor_trace: &[f32]) -> usize { - let max_p = precursor_trace.iter().copied().fold(0.0f32, f32::max); - - let mut best_val = f32::NEG_INFINITY; - let mut best_idx = 0usize; - - let n = cosine_profile.len().min(precursor_trace.len()); - for t in 0..n { - let p_factor = if max_p > 0.0 { - 0.5 + precursor_trace[t] / max_p - } else { - 0.5 - }; - let joint = cosine_profile[t] * p_factor; - if joint > best_val { - best_val = joint; - best_idx = t; - } - } - best_idx -} - -/// Compute all 11 apex-local features at the joint apex (METHODS.md Section 3.4). +/// Compute all 11 apex-local features at the apex (METHODS.md Section 3.4). /// /// `fragments` and `precursors` are the raw chromatogram data. /// `expected` contains both fragment and precursor predicted intensities. /// `cosine_profile` is C(t) = cos(t)^3 * I(t). /// `precursor_trace` is the summed precursor intensity trace. -/// `joint_apex` is the cycle index from `find_joint_apex`. +/// `joint_apex` is the apex cycle index (Pass 1's weighted-apex_profile pick). /// `n_cycles` is the total number of cycles in the extraction window. pub fn compute_apex_features( fragments: &MzMajorIntensityArray, @@ -1155,47 +1129,6 @@ mod tests { assert!(result.scribe_au > 0.0); } - // --- Test 6: joint apex --- - - #[test] - fn test_joint_apex_with_precursor() { - let n = 30; - // Fragment profile peaks at 15 - let cosine_profile: Vec = (0..n) - .map(|i| { - let x = (i as f32 - 15.0) / 3.0; - (-x * x / 2.0).exp() * 10.0 - }) - .collect(); - // Precursor trace also peaks at 15 - let precursor_trace: Vec = (0..n) - .map(|i| { - let x = (i as f32 - 15.0) / 4.0; - (-x * x / 2.0).exp() * 5.0 - }) - .collect(); - - let apex = find_joint_apex(&cosine_profile, &precursor_trace); - assert_eq!(apex, 15); - } - - #[test] - fn test_joint_apex_without_precursor() { - let n = 30; - let cosine_profile: Vec = (0..n) - .map(|i| { - let x = (i as f32 - 15.0) / 3.0; - (-x * x / 2.0).exp() * 10.0 - }) - .collect(); - // Zero precursor - let precursor_trace = vec![0.0f32; n]; - - let apex = find_joint_apex(&cosine_profile, &precursor_trace); - // Should still find the cosine peak at 15 - assert_eq!(apex, 15); - } - // --- Test 7: individual feature tests --- #[test]