From 1dd0c9f8862f986ee6e778412039da574079118c Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Tue, 19 May 2026 16:00:39 -0700 Subject: [PATCH 01/16] HPC chain Stage 7 correctness: 2nd-pass Percolator at --join-at-pass=2, file_stems-aware reconciliation hash, and Stage 6 hard-error gates Closes a set of HPC-distribution-path correctness bugs surfaced by a new strict bit-parity test that compares the 4-step HPC chain (raw workers + first-join + per-file rescore workers + 2nd-join merge) against a straight-through pipeline run on the same input with real mzMLs and calibration. Four logical changes, validated end-to-end: Bug A: Run 2nd-pass Percolator at --join-at-pass=2 when sidecars are missing The HPC distribution path writes reconciled .scores.parquet files from per-file Stage 6 workers but does NOT write .2nd-pass.fdr_scores.bin sidecars (2nd-pass FDR is a cross-file join step the per-file workers can't perform). The merge node at --join-at-pass=2 previously loaded reconciled parquets and went straight to protein FDR + blib output, silently using stale 1st-pass scores for every rescored entry. Compared against a straight-through pipeline run, the chain lost ~25% of precursors -- the 2nd-pass Percolator scores weren't there to rank entries whose 1st-pass scores had been reset to 0 during reconciliation. After loading reconciled state and reloading any 2nd-pass sidecars that DO exist, if any file is missing its sidecar, run run_percolator_fdr on the post-compaction per_file_entries with the first_pass_base_ids restriction. The new scores are persisted as 2nd-pass sidecars so subsequent --join-at-pass=2 invocations against the same reconciled parquets short-circuit. Bug B: file_stems envelope + reconciliation_parameter_hash_for_stems Per-file Stage 6 rescore workers stamp osprey.reconciliation_hash into their reconciled .scores.parquet footer so the downstream --join-at-pass=2 merge can validate the config hasn't drifted. Previously the worker computed that hash from OspreyConfig.input_files, which in worker mode is the single parquet it was given -- so the worker stamped a single-file hash that the merge node, computing the hash over the full join file set, correctly rejected. reconciliation.json is now v2 with a sorted + deduped file_stems field carrying the planner's full join-wide file set; the new OspreyConfig::reconciliation_parameter_hash_for_stems overload takes an explicit stems slice so the worker can compute the join-wide hash from the envelope. The worker hydrate path validates that every sibling envelope carries the same set. v1 envelopes deserialize cleanly via #[serde(default)]; the worker falls back to its input_files stems when the hydrated set is empty. Bug B residual: worker compaction must keep cross-file-rescued entries After the worker loads its reconciled parquet and the planner's envelope, it compacts per_file_entries down to base_ids that passed first-pass FDR locally. That filter dropped entries whose own file failed local FDR but whose peptide passed FDR in a sibling file -- those entries qualify for the cross-file consensus rescue compute_consensus_rts performs in the planner, so the planner correctly emits ForcedIntegration / UseCwtPeak actions for them. The worker's local-FDR filter then silently discarded those actions (logged as "reconciliation actions dropped: N" at warn level) before the rescore engine could apply them. Bisected to entry 17365 in Stellar file 20: fails local peptide-FDR (q=1.0) and local protein-FDR (q=1.0) but passes experiment-level FDR at 0.39% via cross-file consensus; planner emits forced_integration; worker dropped the action. first_pass_base_ids is now the UNION of the local-FDR filter AND the entry_ids that have reconciliation actions in reconciliation_actions_pre. Bug D: Make missing/unparseable Stage 6 calibration.json a hard error Both the per-file rescore worker (rescore.rs::run_rescore) and the in-process pipeline path (pipeline.rs::rescore_per_file_loop) previously treated a missing or unparseable sibling .calibration.json as a soft fallback: log a warning and proceed with empty / None calibration state. Stage 6 with no MS calibration produces wrong-mass matches that look like real rescore work but with corrupted scores; the failure mode is "the worker silently emits a reconciled parquet whose rows are nonsense" with no surface error. Both paths now hard-error with OspreyError::config citing the missing or unparseable file. The cal_params type in rescore_per_file_loop changes from Option to CalibrationParams; the three run_search call sites in that function go from cal_params.as_ref() to Some(&cal_params). Stages 1-4 call sites retain Option<...> because calibration is genuinely optional during cold-start scoring; only Stage 6 fails hard on absence. A missing mzML at Stage 6 already errors via the existing load_all_spectra error propagation -- no additional check needed. Validation: After all four changes, the new Compare-Stage7-Rehydration-Strict.ps1 test reports Stellar 3-file truth=chain=60373 precursors with Stage 7 protein FDR dump SHA D66E5CF0FCC96E9F identical on both sides and blib SQL content matching at 1e-9 tolerance; Astral 3-file truth=chain=165288 precursors with Stage 7 protein FDR dump SHA 50569B3C93E89FDA identical and blib matching. Cross-impl validation in a companion pwiz PR confirms the C# side reads + writes the v2 envelope and matches Rust hash computation on both datasets end-to-end. --- crates/osprey-core/src/config.rs | 43 +++-- crates/osprey/src/pipeline.rs | 216 ++++++++++++++++++++++--- crates/osprey/src/reconciliation_io.rs | 42 ++++- crates/osprey/src/rescore.rs | 104 ++++++++++-- 4 files changed, 361 insertions(+), 44 deletions(-) diff --git a/crates/osprey-core/src/config.rs b/crates/osprey-core/src/config.rs index 503c42a..52c1ba1 100644 --- a/crates/osprey-core/src/config.rs +++ b/crates/osprey-core/src/config.rs @@ -580,8 +580,35 @@ n_threads: 0 # 0 = auto-detect } /// Compute SHA-256 hash of parameters that affect reconciliation. - /// Includes the search hash (if search changed, reconciliation is also invalid). + /// Includes the search hash (if search changed, reconciliation is also + /// invalid). The file set affects the consensus RTs that drive + /// reconciliation actions, so it is part of the hash. Uses + /// `self.input_files` stems by default; the per-file rescore worker + /// (`--join-at-pass=1 --no-join`) has only one input parquet in + /// `input_files` but participated in a multi-file join, so it must + /// instead compute the hash with the full file_stems list from its + /// reconciliation.json via [`Self::reconciliation_parameter_hash_for_stems`]. pub fn reconciliation_parameter_hash(&self) -> String { + let stems: Vec = self + .input_files + .iter() + .filter_map(|p| { + p.file_stem() + .and_then(|s| s.to_str()) + .map(|s| s.to_string()) + }) + .collect(); + self.reconciliation_parameter_hash_for_stems(&stems) + } + + /// Variant of [`Self::reconciliation_parameter_hash`] that accepts the + /// file_stems list explicitly. The per-file rescore worker uses this + /// with the file_stems carried in `reconciliation.json` so the hash + /// it writes into the reconciled parquet matches the hash the + /// downstream `--join-at-pass=2` merge node will compute (which is + /// based on the full set of `--input-scores` it receives, not the + /// single parquet the worker rescored). + pub fn reconciliation_parameter_hash_for_stems(&self, file_stems: &[String]) -> String { let mut hasher = Sha256::new(); hasher.update(self.search_parameter_hash().as_bytes()); hasher @@ -594,17 +621,11 @@ n_threads: 0 # 0 = auto-detect .as_bytes(), ); hasher.update(format!("run_fdr:{}\n", self.run_fdr).as_bytes()); - // File set affects consensus RTs - let mut stems: Vec = self - .input_files - .iter() - .filter_map(|p| { - p.file_stem() - .and_then(|s| s.to_str()) - .map(|s| s.to_string()) - }) - .collect(); + // File set affects consensus RTs. Sort + dedup so the same logical + // set produces the same hash regardless of caller iteration order. + let mut stems = file_stems.to_vec(); stems.sort(); + stems.dedup(); hasher.update(format!("file_stems:{:?}\n", stems).as_bytes()); format!("{:x}", hasher.finalize()) } diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index da8a1f8..285791b 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -342,7 +342,24 @@ fn build_scores_metadata(config: &OspreyConfig) -> Vec Vec { +/// Build the reconciled-parquet metadata block. +/// +/// `file_stems_override`: when present, used to compute the +/// reconciliation hash via +/// [`OspreyConfig::reconciliation_parameter_hash_for_stems`]. The +/// per-file rescore worker (`--join-at-pass=1 --no-join`) passes the +/// file_stems list from reconciliation.json here so the hash it writes +/// matches the multi-file hash the downstream `--join-at-pass=2` merge +/// node will compute. Pass `None` from the in-process pipeline path — +/// `config.input_files` already reflects the full set. +fn build_reconciled_metadata( + config: &OspreyConfig, + file_stems_override: Option<&[String]>, +) -> Vec { + let recon_hash = match file_stems_override { + Some(stems) if !stems.is_empty() => config.reconciliation_parameter_hash_for_stems(stems), + _ => config.reconciliation_parameter_hash(), + }; vec![ parquet::file::metadata::KeyValue { key: META_OSPREY_VERSION.to_string(), @@ -362,7 +379,7 @@ fn build_reconciled_metadata(config: &OspreyConfig) -> Vec, config: &OspreyConfig, seq_interner: &mut HashSet>, + // Full set of file stems that participated in the first-join phase + // (the join that produced the reconciliation actions this loop applies). + // Used to compute the reconciliation hash written into the reconciled + // parquet's metadata so the downstream `--join-at-pass=2` merge node + // accepts it. When empty, falls back to `config.input_files` stems + // (the in-process pipeline path always has the full set there). + join_file_stems: &[String], ) -> Result { use crate::reconciliation::ReconcileAction; let total_reconciliation: usize = reconciliation_actions @@ -3023,14 +3047,41 @@ pub(crate) fn rescore_per_file_loop( ); } - let cal_params: Option = input_file.parent().and_then(|input_dir| { + // Hard-error on missing or unreadable calibration JSON. The + // calibration sidecar is required by Stage 6 for MS1/MS2 mass + // calibration of the rescore search; silently proceeding with + // uncalibrated MS scoring would produce wrong scores without + // surfacing the cause. In-process this should always exist + // (Stage 2 just wrote it); workers already gate on it in + // `run_rescore` so this catch is belt-and-suspenders for both + // paths. + let cal_params: CalibrationParams = { + let input_dir = input_file.parent().ok_or_else(|| { + OspreyError::config(format!( + "rescore_per_file_loop: cannot derive parent directory from input path `{}`. \ + Stage 6 needs to locate the Stage 1-4 calibration sidecar next to the mzML.", + input_file.display() + )) + })?; let cal_path = calibration_path_for_input(input_file, input_dir); - if cal_path.exists() { - load_calibration(&cal_path).ok() - } else { - None + if !cal_path.exists() { + return Err(OspreyError::config(format!( + "rescore_per_file_loop: required calibration JSON not found at `{}` (input \ + file: `{}`). Stage 6 needs the Stage 1-4 calibration sidecar to rescore.", + cal_path.display(), + input_file.display() + ))); } - }); + load_calibration(&cal_path).map_err(|e| { + OspreyError::config(format!( + "rescore_per_file_loop: failed to read calibration JSON `{}`: {}. The file \ + exists but could not be parsed -- check that it was written by a matching \ + Osprey version.", + cal_path.display(), + e + )) + })? + }; // --- Re-score existing entries (consensus + reconciliation) --- let mut overlay: HashMap = HashMap::new(); @@ -3039,7 +3090,7 @@ pub(crate) fn rescore_per_file_loop( &subset_library, &spectra, &ms1_index, - cal_params.as_ref(), + Some(&cal_params), rt_cal, config, file_name, @@ -3066,7 +3117,7 @@ pub(crate) fn rescore_per_file_loop( &gap_fill_library, &spectra, &ms1_index, - cal_params.as_ref(), + Some(&cal_params), rt_cal, &gap_config, file_name, @@ -3121,7 +3172,7 @@ pub(crate) fn rescore_per_file_loop( &forced_library, &spectra, &ms1_index, - cal_params.as_ref(), + Some(&cal_params), rt_cal, config, file_name, @@ -3237,7 +3288,7 @@ pub(crate) fn rescore_per_file_loop( } } // Write back to Parquet with reconciliation metadata - let recon_metadata = build_reconciled_metadata(config); + let recon_metadata = build_reconciled_metadata(config, Some(join_file_stems)); let codec = parquet_compression_codec(config.parquet_compression); if let Err(e) = write_scores_parquet_with_metadata( cache_path, @@ -3257,7 +3308,7 @@ pub(crate) fn rescore_per_file_loop( let scores_path = per_file_cache_paths.get(file_name.as_str()); if let Some(cache_path) = scores_path { if let Ok(full_entries) = load_scores_parquet(cache_path) { - let recon_metadata = build_reconciled_metadata(config); + let recon_metadata = build_reconciled_metadata(config, Some(join_file_stems)); let codec = parquet_compression_codec(config.parquet_compression); if let Err(e) = write_scores_parquet_with_metadata( cache_path, @@ -3285,6 +3336,11 @@ pub(crate) fn rescore_per_file_loop( } pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { + // Sliding timer for [STAGE-WALL] per-stage perf markers. Reset at + // each stage boundary (stage1to4 -> stage5 -> stage6 -> stage7 -> blib). + // Parsed by Measure-Pipeline.ps1 for cross-impl perf comparison. + let mut stage_marker = std::time::Instant::now(); + // Validate config: --join-only feeds per-file state from --input-scores // parquets and ignores --input-files; default + --no-join modes both // require mzML inputs. @@ -4110,6 +4166,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { percolator::compute_fdr_from_stubs(&mut per_file_entries, config.run_fdr, None); } } else { + log::info!( + "[STAGE-WALL] stage1to4: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + // Dispatch FDR control based on method log::info!(""); log::info!( @@ -4158,11 +4220,16 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // = 1.0` and silently break protein-rescue parity in the worker. } - // Stage 5 diagnostic dump. Gated by OSPREY_DUMP_PERCOLATOR=1; exits the - // process when OSPREY_PERCOLATOR_ONLY=1 is also set. Writes all 4 q-values - // plus the SVM score and PEP for every FdrEntry, before compaction drops - // any rows, so the cross-impl diff sees both targets and decoys. + // Stage 5 diagnostic dump. Gated by OSPREY_DUMP_PERCOLATOR=1. Writes + // all 4 q-values plus the SVM score and PEP for every FdrEntry, before + // compaction drops any rows, so the cross-impl diff sees both targets + // and decoys. crate::diagnostics::dump_stage5_percolator(&per_file_entries); + // OSPREY_PERCOLATOR_ONLY exits after Stage 5 work completes, + // independently of whether the dump ran. Lets us measure production + // stage5 wall without paying the dump cost. Matches the C# decoupling + // in OspreySharp/Tasks/FirstJoinTask.cs. + osprey_core::diagnostics::exit_if_only("OSPREY_PERCOLATOR_ONLY", "Stage 5"); // First-pass protein FDR (picked-protein, Savitski 2015). // @@ -4352,13 +4419,30 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // 1st-pass-derived q-values, swap entry.score over to 2nd-pass scores so // the second-pass q-value recompute below sees the right rank ordering. // No-op if no 2nd-pass sidecar exists for a given file (best-effort). + // + // When NO 2nd-pass sidecar is present for ANY file (the HPC + // distribution case — the per-file rescore worker writes reconciled + // parquets but does NOT write 2nd-pass sidecars, since 2nd-pass + // Percolator is a cross-file join step that the per-file worker + // can't perform), we run 2nd-pass Percolator here so the merge + // node's output matches the straight-through pipeline. Without this, + // --join-at-pass=2 would use stale first-pass scores for rescored + // entries (which had their scores reset to 0 during reconciliation + // and never re-scored), and the final blib would silently lose + // ~25% of the precursors a straight-through run produces. if config.expect_reconciled_input { let mut reloaded = 0usize; + let mut missing_sidecars = 0usize; for (file_name, entries) in per_file_entries.iter_mut() { - if let Some(p2) = pass2_sidecar_paths.get(file_name) { + let p2 = pass2_sidecar_paths.get(file_name); + if let Some(p2) = p2 { if load_fdr_scores_sidecar(p2, entries, 2) { reloaded += 1; + } else { + missing_sidecars += 1; } + } else { + missing_sidecars += 1; } } log::debug!( @@ -4366,6 +4450,61 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { reloaded, per_file_entries.len() ); + if missing_sidecars > 0 { + log::info!( + "--join-at-pass=2: {}/{} file(s) lack a 2nd-pass sidecar — running second-pass FDR \ + to compute scores from reconciled features (HPC distribution path).", + missing_sidecars, + per_file_entries.len() + ); + log::info!(""); + log::info!("Second-pass FDR"); + let empty_overlay: RescoreOverlay = HashMap::new(); + match config.fdr_method { + FdrMethod::Percolator => { + run_percolator_fdr( + &mut per_file_entries, + &per_file_cache_paths, + &config, + &empty_overlay, + Some(&first_pass_base_ids), + )?; + } + FdrMethod::Mokapot => { + run_mokapot_fdr( + &mut per_file_entries, + &per_file_cache_paths, + &mokapot, + &pin_files, + &mokapot_dir, + &config, + )?; + } + FdrMethod::Simple => { + for (_, entries) in per_file_entries.iter_mut() { + apply_simple_fdr(entries, config.run_fdr)?; + } + } + } + crate::trace::log_fdr_qvalues(&per_file_entries, "second-pass"); + + // Persist the 2nd-pass scores we just computed as a resume + // cache for future --join-at-pass=2 invocations against the + // same reconciled parquets. Mirrors the persist call inside + // the in-process reconciliation block below; gated by + // multi-file because reconciliation only fires there. + let has_reconciliation = config.reconciliation.enabled && per_file_entries.len() > 1; + if has_reconciliation { + log::debug!("Persisting 2nd-pass FDR scores as resume cache..."); + let _ = persist_fdr_scores( + &per_file_entries, + &config, + fdr_scores_path_pass2, + "2nd-pass", + 2, + ); + } + } } // Post-FDR re-scoring: multi-charge consensus + inter-replicate reconciliation. @@ -4457,6 +4596,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { } }); + log::info!( + "[STAGE-WALL] stage5: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + // 1. Multi-charge consensus: compute per-file rescore targets // Groups by (peptide, file). If at least one charge state passes FDR, // the best-scoring charge state defines the consensus peak; other @@ -4632,6 +4777,18 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // joined work above. let search_hash = config.search_parameter_hash(); let library_hash = config.library_identity_hash(); + // The full set of file stems that participated in this first-join + // phase. Embedded in every per-file reconciliation.json so the + // downstream per-file rescore worker can compute the same + // reconciliation hash the `--join-at-pass=2` merge node will + // validate against (which is computed over all parquets at the + // 2nd-join node, not the worker's single input). Without this + // the worker writes a single-stem hash and `--join-at-pass=2` + // rejects the reconciled parquet. + let join_file_stems: Vec = per_file_entries + .iter() + .map(|(fname, _)| fname.clone()) + .collect(); // Pre-group reconciliation actions by file name to avoid the // O(num_files * num_actions) walk that the previous // implementation performed inside `from_planner_output` (one @@ -4669,6 +4826,7 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { refined_calibrations.get(file_name), &search_hash, &library_hash, + &join_file_stems, ); if let Err(e) = crate::reconciliation_io::write_reconciliation_file(&recon_path, &recon_file) @@ -4733,6 +4891,11 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { &file_name_to_idx, &config, &mut seq_interner, + // In-process pipeline: config.input_files already has the + // full set, so the empty-slice override makes + // build_reconciled_metadata fall back to + // config.reconciliation_parameter_hash(). + &[], )?; // Cross-impl bisection seam: dump the per-precursor q-values @@ -4752,6 +4915,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { total_gap_cwt, total_gap_forced, ); + log::info!( + "[STAGE-WALL] stage6: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + log::info!(""); log::info!("Second-pass FDR"); match config.fdr_method { @@ -5119,6 +5288,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // to the final destination. This avoids SQLite locking issues on network // filesystems. if !plan_entries.is_empty() { + log::info!( + "[STAGE-WALL] stage7: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + log::info!("Writing blib to {}", config.output_blib.display()); let final_path = &config.output_blib; @@ -5131,6 +5306,11 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // Move to final destination (safe copy for network filesystems) osprey_core::copy_and_verify(&temp_path, final_path)?; + + log::info!( + "[STAGE-WALL] blib: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); } else { log::warn!("No peptides passed FDR threshold, skipping blib output"); } diff --git a/crates/osprey/src/reconciliation_io.rs b/crates/osprey/src/reconciliation_io.rs index fa66b9b..cb295e3 100644 --- a/crates/osprey/src/reconciliation_io.rs +++ b/crates/osprey/src/reconciliation_io.rs @@ -141,9 +141,27 @@ impl<'a> Formatter for RoundtripPrettyFormatter<'a> { /// Top-level JSON envelope. Field declaration order is alphabetical so /// `serde_json::to_writer_pretty` emits keys in alphabetical order and /// matches the C# emitter byte-for-byte. +/// +/// `file_stems` (added 2026-05-19, format_version 2) carries the full +/// set of file stems that participated in the first-join phase. The +/// per-file rescore worker (`--join-at-pass=1 --no-join`) consumes only +/// its own parquet, so its `config.input_files` reflects a single +/// stem. To produce a reconciled parquet that the downstream +/// `--join-at-pass=2` merge node can validate (its +/// `reconciliation_parameter_hash` is computed over all parquets it +/// receives), the worker uses this list to compute the same multi-stem +/// hash the original join used. v1 files (no `file_stems`) deserialize +/// with an empty list; the worker will then fall back to its own +/// `config.input_files` stems, preserving v1 behavior. #[derive(Debug, Clone, Serialize, Deserialize)] #[serde(deny_unknown_fields)] pub struct ReconciliationFile { + /// Sorted list of file stems that participated in the first-join phase. + /// `#[serde(default)]` deserializes v1 files (which omit this field) as + /// an empty list; the worker code falls back to `config.input_files` + /// stems in that case for backward compatibility. + #[serde(default)] + pub file_stems: Vec, pub forced_integration_actions: Vec, pub format_version: u32, pub gap_fill_targets: Vec, @@ -198,7 +216,16 @@ pub struct RefinedRtCalibrationJson { } /// Current format version. Bump on incompatible schema changes. -pub const RECONCILIATION_FORMAT_VERSION: u32 = 1; +/// +/// v1: initial format. +/// v2: added `file_stems` so per-file rescore workers can compute the +/// reconciliation parameter hash that the downstream +/// `--join-at-pass=2` merge node expects (which is computed over +/// all files participating in the join, not the worker's single +/// parquet). Old v1 files deserialize with empty `file_stems` via +/// `#[serde(default)]`; the worker falls back to its +/// `config.input_files` stems for those (preserving v1 behavior). +pub const RECONCILIATION_FORMAT_VERSION: u32 = 2; impl ReconciliationFile { /// Build the wire envelope for a single file. Filters @@ -211,6 +238,7 @@ impl ReconciliationFile { /// should prefer [`Self::from_planner_output_pre_grouped`] together /// with a single up-front grouping to keep the total cost /// O(num_actions) rather than O(num_files * num_actions). + #[allow(clippy::too_many_arguments)] pub fn from_planner_output( file_name: &str, file_entries: &[FdrEntry], @@ -219,6 +247,7 @@ impl ReconciliationFile { refined_calibration: Option<&RTCalibration>, search_hash: &str, library_hash: &str, + file_stems: &[String], ) -> Self { let file_actions: Vec<(usize, &ReconcileAction)> = reconciliation_actions .iter() @@ -237,6 +266,7 @@ impl ReconciliationFile { refined_calibration, search_hash, library_hash, + file_stems, ) } @@ -246,6 +276,7 @@ impl ReconciliationFile { /// included or not — they are filtered here either way). Use this /// when emitting envelopes for many files so the per-file cost stays /// O(actions_for_this_file) rather than O(total_actions). + #[allow(clippy::too_many_arguments)] pub fn from_planner_output_pre_grouped( file_entries: &[FdrEntry], file_actions: &[(usize, &ReconcileAction)], @@ -253,7 +284,14 @@ impl ReconciliationFile { refined_calibration: Option<&RTCalibration>, search_hash: &str, library_hash: &str, + file_stems: &[String], ) -> Self { + // Sort and deduplicate the file_stems for deterministic output + // — matches the sort+dedup used by reconciliation_parameter_hash + // on the read side, so the worker hash matches the merge node hash. + let mut stems: Vec = file_stems.to_vec(); + stems.sort(); + stems.dedup(); let mut use_cwt: Vec = Vec::new(); let mut forced: Vec = Vec::new(); for (vec_idx, action) in file_actions { @@ -317,6 +355,7 @@ impl ReconciliationFile { }); Self { + file_stems: stems, forced_integration_actions: forced, format_version: RECONCILIATION_FORMAT_VERSION, gap_fill_targets: gap, @@ -423,6 +462,7 @@ mod tests { fn sample_file() -> ReconciliationFile { ReconciliationFile { + file_stems: vec!["fileA".to_string(), "fileB".to_string()], forced_integration_actions: vec![ ForcedIntegrationEntry { entry_id: 200, diff --git a/crates/osprey/src/rescore.rs b/crates/osprey/src/rescore.rs index b58e67b..f3d22ae 100644 --- a/crates/osprey/src/rescore.rs +++ b/crates/osprey/src/rescore.rs @@ -81,6 +81,13 @@ pub struct RescoreInputs { /// gap-fill stubs reuse the same `Arc` identities that the /// hydrated stubs already hold. pub seq_interner: HashSet>, + /// Full set of file stems that participated in the first-join phase + /// that produced the reconciliation actions in this RescoreInputs. + /// Read from `reconciliation.json`'s `file_stems` field (v2+); + /// empty for v1 files. Used to compute the reconciliation hash + /// written into the reconciled parquet's metadata so the downstream + /// `--join-at-pass=2` merge node will accept it. + pub join_file_stems: Vec, } impl RescoreInputs { @@ -131,6 +138,12 @@ pub fn hydrate_for_rescore(config: &OspreyConfig) -> Result { let mut refined_calibrations: HashMap = HashMap::new(); let mut per_file_gap_fill: HashMap> = HashMap::new(); let mut reconciliation_actions: HashMap<(String, usize), ReconcileAction> = HashMap::new(); + // Captured from the first per-file reconciliation.json envelope; every + // file's envelope must declare the same join file_stems (all came + // from the same first-join run). Used downstream to compute the + // reconciliation hash the worker stamps into the reconciled parquet's + // metadata. + let mut join_file_stems: Vec = Vec::new(); for parquet_path in parquet_paths { // The parquet stem (with `.scores` stripped) is the canonical @@ -184,6 +197,32 @@ pub fn hydrate_for_rescore(config: &OspreyConfig) -> Result { )) })?; + // Capture the join file_stems from the first envelope. Every + // file's envelope should declare the same list (they all came + // from the same first-join phase). Validate consistency + // file-to-file so an accidental mix of boundary files from + // different join runs fails loudly here rather than producing + // a reconciled parquet whose hash doesn't match anything. + if join_file_stems.is_empty() { + join_file_stems = envelope.file_stems.clone(); + join_file_stems.sort(); + join_file_stems.dedup(); + } else { + let mut these = envelope.file_stems.clone(); + these.sort(); + these.dedup(); + if !envelope.file_stems.is_empty() && these != join_file_stems { + return Err(OspreyError::config(format!( + "hydrate_for_rescore: reconciliation.json file_stems mismatch — {} \ + declares {:?}, expected {:?} (boundary files from different first-join \ + runs cannot be mixed in one rescore worker invocation).", + recon_path.display(), + these, + join_file_stems, + ))); + } + } + // 3a. Build entry_id → vec_idx map from the loaded stubs so the // planner's entry_id-keyed actions can be rehomed onto the // in-memory shape (file_name, vec_idx) the rescore engine @@ -303,6 +342,7 @@ pub fn hydrate_for_rescore(config: &OspreyConfig) -> Result { refined_calibrations, per_file_gap_fill, seq_interner, + join_file_stems, }) } @@ -383,19 +423,23 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( }; let cal_path = calibration_path_for_input(input_file, input_dir); if !cal_path.exists() { - continue; + return Err(OspreyError::config(format!( + "run_rescore: required calibration JSON not found at {} (input file: {}). \ + Stage 6 needs the Stage 1-4 calibration sidecar to rescore; without it the \ + worker would silently produce no reconciliation re-scores. Run Stages 1-4 \ + first or fix the path.", + cal_path.display(), + input_file.display() + ))); } - let cal_params = match load_calibration(&cal_path) { - Ok(p) => p, - Err(e) => { - log::warn!( - "run_rescore: failed to read calibration JSON {}: {}", - cal_path.display(), - e - ); - continue; - } - }; + let cal_params = load_calibration(&cal_path).map_err(|e| { + OspreyError::config(format!( + "run_rescore: failed to read calibration JSON {}: {}. The file exists but \ + could not be parsed — check that it was written by a matching Osprey version.", + cal_path.display(), + e + )) + })?; if let Some(ref mp) = cal_params.rt_calibration.model_params { match RTCalibration::from_model_params(mp, cal_params.rt_calibration.residual_sd) { Ok(rt_cal) => { @@ -420,6 +464,7 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( refined_calibrations, per_file_gap_fill, mut seq_interner, + join_file_stems, } = hydrate_for_rescore(&config)?; // Cross-impl bisection seam (mirrors the dump call from @@ -463,6 +508,30 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( } } + // Union with entries that have reconciliation actions. The planner + // emits actions for entries that pass FDR via cross-file consensus + // rescue (`compute_consensus_rts` upgrades a peptide if it passes + // FDR in any file in the experiment, even if it fails locally). + // Without this union, the worker would drop those entries via + // local compaction and the planner's actions would be silently + // discarded — the per-file rescore would diverge from the + // straight-through in-process pipeline. The base_ids come from + // the pre-compaction `per_file_entries` keyed by (file_name, idx) + // in `reconciliation_actions_pre`. + { + let entries_by_name: HashMap<&str, &Vec> = per_file_entries + .iter() + .map(|(name, entries)| (name.as_str(), entries)) + .collect(); + for ((file_name, idx), _action) in reconciliation_actions_pre.iter() { + if let Some(entries) = entries_by_name.get(file_name.as_str()) { + if let Some(e) = entries.get(*idx) { + first_pass_base_ids.insert(e.entry_id & 0x7FFF_FFFF); + } + } + } + } + // Save (file, entry_id) → action before per_file_entries shrinks // so we can rebuild (file, new_vec_idx) → action below. Use a // file_name → entries lookup map so the per-action lookup is O(1) @@ -519,8 +588,8 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( } log::info!( - "Worker compaction: {} -> {} entries ({} survived peptide-FDR or protein-rescue), \ - {} reconciliation actions retained", + "Worker compaction: {} -> {} entries ({} survived peptide-FDR, protein-rescue, \ + or cross-file consensus rescue), {} reconciliation actions retained", entries_before, entries_after, first_pass_base_ids.len(), @@ -570,6 +639,13 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( &file_name_to_idx, &config, &mut seq_interner, + // The per-file rescore worker's config.input_files reflects + // only this worker's single parquet. To stamp the reconciled + // parquet with the multi-file reconciliation hash the + // downstream `--join-at-pass=2` merge node will validate + // against, pass the full join file_stems list extracted from + // reconciliation.json. + &join_file_stems, )?; // Cross-impl bisection seam: dump the per-precursor q-values From 1fe4ff73c1383a9d2065ed05af325ea9de8453e7 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 09:36:40 -0700 Subject: [PATCH 02/16] mzML: read isolation-window cvParams at f64 precision mzdata 0.63's reader pipes isolation window cvParam values (MS:1000827 target m/z, MS:1000828 lower offset, MS:1000829 upper offset) through `param.to_f32()`. At f32 precision the ULP at m/z 500 is ~6e-5, so any window edge whose XML text value lands between two f32-representable values quantizes differently when round-tripped. On Stellar 3-file this surfaces as a single anomalous DIA window out of 125 whose upper edge reads 512.481934 in Rust vs 512.481903 in OspreySharp (~3e-5 m/z drift). 1,732 library entries fell in that window and inherited the drift in every downstream calibration artifact. Cross-impl `Compare-Stage1to4-Strict.ps1` previously reported 1,732 row diffs at 3.1e-5 on `cal_windows.iso_upper`; after this change the cal_windows dump is bit-equal cross-impl (0 diffs, max_diff=0.000e+000). The fix is local to osprey-io and self-contained: a one-pass streaming scan with quick-xml extracts every `` block's cvParam values as f64 strings before the mzdata pass, and a small helper overrides mzdata's f32-quantized `lower_bound`/`upper_bound` with the f64 cvParams at the two MS2 parsing sites (`convert_spectrum` and `load_all_spectra`). The fallback path keeps mzdata's values when a spectrum's pre-parsed cvParams are missing, so older mzML converters that omit MS:1000827/828/829 still work as before. When mzdata moves to f64 storage upstream this whole pre-pass and the override helper can be deleted in a single commit; the rationale is spelled out in the function docs and pinned to the workspace quick-xml = "0.30" declaration. --- Cargo.lock | 1 + Cargo.toml | 6 + crates/osprey-io/Cargo.toml | 1 + crates/osprey-io/src/mzml/parser.rs | 242 ++++++++++++++++++++++------ 4 files changed, 201 insertions(+), 49 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 06310e6..adc01ad 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1516,6 +1516,7 @@ dependencies = [ "log", "mzdata", "osprey-core", + "quick-xml", "rusqlite", "serde", "tempfile", diff --git a/Cargo.toml b/Cargo.toml index fd2f4de..d1bfd05 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -38,6 +38,12 @@ parquet = "58" # File I/O mzdata = "0.63" +# Used by osprey-io::mzml to pre-parse isolation-window cvParams as +# f64, bypassing mzdata 0.63's `param.to_f32()` quantization which +# produces ~3e-5 m/z drift on isolation bounds that land between two +# f32-representable values. mzdata transitively pins quick-xml 0.30, +# so this version matches and does not bloat the dep graph. +quick-xml = "0.30" rusqlite = { version = "0.39", features = ["bundled"] } csv = "1.3" # Use flate2's `zlib-default` feature (vendored stock zlib via libz-sys) diff --git a/crates/osprey-io/Cargo.toml b/crates/osprey-io/Cargo.toml index dfd2a89..e1f001b 100644 --- a/crates/osprey-io/Cargo.toml +++ b/crates/osprey-io/Cargo.toml @@ -11,6 +11,7 @@ repository.workspace = true [dependencies] osprey-core = { workspace = true } mzdata = { workspace = true } +quick-xml = { workspace = true } rusqlite = { workspace = true } csv = { workspace = true } flate2 = { workspace = true } diff --git a/crates/osprey-io/src/mzml/parser.rs b/crates/osprey-io/src/mzml/parser.rs index 0768a79..3daba95 100644 --- a/crates/osprey-io/src/mzml/parser.rs +++ b/crates/osprey-io/src/mzml/parser.rs @@ -14,10 +14,168 @@ use mzdata::io::mzml::MzMLReader; use mzdata::prelude::*; use mzdata::spectrum::RefPeakDataLevel; use osprey_core::{IsolationWindow, MS1Spectrum, OspreyError, Result, Spectrum, SpectrumSource}; +use quick_xml::events::Event; +use quick_xml::reader::Reader; +use std::collections::HashMap; use std::fs::File; use std::io::BufReader; use std::path::{Path, PathBuf}; +/// Raw f64 isolation window cvParam values lifted directly from the +/// mzML XML, indexed by spectrum index (matching mzdata's +/// `MultiLayerSpectrum.description().index`). Used to override +/// mzdata 0.63's f32-quantized `IsolationWindow.lower_bound` / +/// `upper_bound` values; see [`read_isolation_cvparams_f64`] for the +/// rationale. +#[derive(Default, Clone, Debug)] +struct IsolationCvParams { + /// MS:1000827 isolation window target m/z (f64 from XML text). + target_mz: Option, + /// MS:1000828 isolation window lower offset (positive magnitude). + lower_offset: Option, + /// MS:1000829 isolation window upper offset (positive magnitude). + upper_offset: Option, +} + +/// One-pass streaming scan of an mzML file to extract every +/// `` block's cvParam values at f64 precision, +/// directly from the XML text. mzdata 0.63's reader pipes these +/// values through `param.to_f32()` (see mzdata-0.63 +/// `src/io/mzml/reader.rs:281`), which quantizes through f32 and +/// produces ~3e-5 m/z drift on isolation-window edges that land +/// between two f32-representable values. The C# port +/// (OspreySharp) parses the same cvParams at f64, so the f32 +/// quantization on the Rust side is the sole source of a 1,732-row +/// cross-impl `iso_upper` mismatch on Stellar 3-file. This +/// function reads the same XML text mzdata sees but preserves full +/// f64 precision. +/// +/// The override is wired in via [`make_isolation_window`] at the two +/// MS2 parsing sites in this module. When upstream mzdata moves to +/// f64 storage for `IsolationWindow` bounds, this pre-pass + override +/// can be deleted in one commit. See workspace `Cargo.toml` +/// `quick-xml = "0.30"` line for the lifetime expectation. +fn read_isolation_cvparams_f64(path: &Path) -> Result> { + let mut reader = Reader::from_file(path).map_err(|e| { + OspreyError::MzmlParseError(format!("quick-xml open '{}': {}", path.display(), e)) + })?; + let mut buf = Vec::new(); + let mut result: HashMap = HashMap::new(); + let mut current_index: Option = None; + let mut depth_in_isolation_window: i32 = 0; + let mut current_params = IsolationCvParams::default(); + + loop { + let event = reader.read_event_into(&mut buf).map_err(|e| { + OspreyError::MzmlParseError(format!( + "quick-xml read at pos {}: {}", + reader.buffer_position(), + e + )) + })?; + match event { + Event::Start(ref e) if e.name().as_ref() == b"spectrum" => { + current_index = None; + for attr in e.attributes().flatten() { + if attr.key.as_ref() == b"index" { + if let Ok(s) = std::str::from_utf8(attr.value.as_ref()) { + current_index = s.parse().ok(); + } + } + } + current_params = IsolationCvParams::default(); + } + Event::Start(ref e) if e.name().as_ref() == b"isolationWindow" => { + depth_in_isolation_window += 1; + } + Event::End(ref e) if e.name().as_ref() == b"isolationWindow" => { + depth_in_isolation_window -= 1; + if depth_in_isolation_window == 0 { + if let Some(idx) = current_index { + result.insert(idx, current_params.clone()); + } + } + } + // cvParam elements inside isolationWindow are typically + // self-closing (`Event::Empty`); handle the start variant + // too for defensive parsing. + Event::Empty(ref e) | Event::Start(ref e) + if depth_in_isolation_window > 0 && e.name().as_ref() == b"cvParam" => + { + let mut accession: Option<&[u8]> = None; + let mut value: Option<&[u8]> = None; + let attrs: Vec<_> = e.attributes().flatten().collect(); + for attr in &attrs { + match attr.key.as_ref() { + b"accession" => accession = Some(attr.value.as_ref()), + b"value" => value = Some(attr.value.as_ref()), + _ => {} + } + } + if let (Some(acc), Some(v)) = (accession, value) { + if let (Ok(acc_str), Ok(v_str)) = + (std::str::from_utf8(acc), std::str::from_utf8(v)) + { + if let Ok(f) = v_str.parse::() { + match acc_str { + "MS:1000827" => current_params.target_mz = Some(f), + "MS:1000828" => current_params.lower_offset = Some(f), + "MS:1000829" => current_params.upper_offset = Some(f), + _ => {} + } + } + } + } + } + Event::Eof => break, + _ => {} + } + buf.clear(); + } + Ok(result) +} + +/// Build an [`IsolationWindow`] for a spectrum, preferring f64 cvParam +/// values from the pre-parsed map if available. Falls back to +/// mzdata's f32-quantized bounds when (a) no map is provided, (b) the +/// map has no entry for this scan index, or (c) the entry is missing +/// the requested cvParams (e.g., older mzML converters that omit +/// MS:1000827 / 828 / 829). +fn make_isolation_window( + precursor: &mzdata::spectrum::Precursor, + scan_number: u32, + iso_cv_params: Option<&HashMap>, +) -> Option { + let ion = precursor.ion()?; + let center_fallback = ion.mz; + + if let Some(map) = iso_cv_params { + if let Some(cv) = map.get(&scan_number) { + if cv.target_mz.is_some() || cv.lower_offset.is_some() || cv.upper_offset.is_some() { + let center = cv.target_mz.unwrap_or(center_fallback); + let lower_offset = cv.lower_offset.unwrap_or(12.5); + let upper_offset = cv.upper_offset.unwrap_or(12.5); + return Some(IsolationWindow::new(center, lower_offset, upper_offset)); + } + } + } + + // Fallback: mzdata-quantized values (f32 → f64 cast). + let isolation = &precursor.isolation_window; + let center = center_fallback; + let lower_offset = if isolation.lower_bound > 0.0 { + center - isolation.lower_bound as f64 + } else { + 12.5 + }; + let upper_offset = if isolation.upper_bound > 0.0 { + isolation.upper_bound as f64 - center + } else { + 12.5 + }; + Some(IsolationWindow::new(center, lower_offset, upper_offset)) +} + /// Some mzML producers emit peaks that are not strictly ascending in m/z /// (observed in a HeLa Astral 3 mz DIA file: ~1 row in 1.7M had a single /// inverted pair of consecutive centroids). Downstream fragment matching @@ -75,12 +233,20 @@ pub struct MzmlReader { reader: MzMLReader>, total_spectra: Option, current_index: usize, + /// Pre-parsed f64 isolation window cvParams, keyed by spectrum + /// index. Overrides mzdata 0.63's f32-quantized bounds in + /// `convert_spectrum`. See [`read_isolation_cvparams_f64`]. + iso_cv_params: HashMap, } impl MzmlReader { /// Open an mzML file for reading pub fn open>(path: P) -> Result { let path = path.as_ref().to_path_buf(); + // First pass: scan the XML for isolation-window cvParams as + // f64. Cheap streaming pass (~5% of bytes parsed) that lets us + // override mzdata 0.63's f32 quantization downstream. + let iso_cv_params = read_isolation_cvparams_f64(&path)?; let file = File::open(&path).map_err(|e| { OspreyError::MzmlParseError(format!("Failed to open file '{}': {}", path.display(), e)) })?; @@ -93,6 +259,7 @@ impl MzmlReader { reader: mzml_reader, total_spectra: None, // Will be determined on first iteration current_index: 0, + iso_cv_params, }) } @@ -116,31 +283,18 @@ impl MzmlReader { scan.start_time // mzdata returns time in minutes }); - // Get isolation window from precursor (now a Vec in mzdata 0.63) - let isolation_window = if let Some(precursor) = desc.precursor.first() { - let ion = match precursor.ion() { - Some(i) => i, - None => return Ok(None), - }; - let isolation = &precursor.isolation_window; - - let center = ion.mz; - // In mzdata 0.63, bounds are f32 directly - let lower_offset = if isolation.lower_bound > 0.0 { - center - isolation.lower_bound as f64 - } else { - 12.5 - }; - let upper_offset = if isolation.upper_bound > 0.0 { - isolation.upper_bound as f64 - center - } else { - 12.5 - }; - - IsolationWindow::new(center, lower_offset, upper_offset) - } else { - // No precursor info - skip this spectrum - return Ok(None); + // Get isolation window from precursor (now a Vec in mzdata 0.63). + // Routes through `make_isolation_window` so we use f64 cvParam + // values from the pre-parsed XML, bypassing mzdata 0.63's f32 + // quantization. See `read_isolation_cvparams_f64`. + let isolation_window = match desc.precursor.first() { + Some(precursor) => { + match make_isolation_window(precursor, scan_number, Some(&self.iso_cv_params)) { + Some(w) => w, + None => return Ok(None), // No selected ion + } + } + None => return Ok(None), // No precursor info }; // Get peaks - use the peaks() method which returns RefPeakDataLevel @@ -333,6 +487,8 @@ impl MS1Index { /// Returns: (MS2 spectra, MS1 index) pub fn load_all_spectra>(path: P) -> Result<(Vec, MS1Index)> { let path = path.as_ref(); + // First pass: f64 isolation-window cvParams from raw XML. + let iso_cv_params = read_isolation_cvparams_f64(path)?; let file = File::open(path).map_err(|e| { OspreyError::MzmlParseError(format!("Failed to open file '{}': {}", path.display(), e)) })?; @@ -400,30 +556,18 @@ pub fn load_all_spectra>(path: P) -> Result<(Vec, MS1In }); } 2 => { - // Process MS2 spectrum - // Get isolation window from precursor - let isolation_window = if let Some(precursor) = desc.precursor.first() { - let ion = match precursor.ion() { - Some(i) => i, - None => continue, - }; - let isolation = &precursor.isolation_window; - - let center = ion.mz; - let lower_offset = if isolation.lower_bound > 0.0 { - center - isolation.lower_bound as f64 - } else { - 12.5 - }; - let upper_offset = if isolation.upper_bound > 0.0 { - isolation.upper_bound as f64 - center - } else { - 12.5 - }; - - IsolationWindow::new(center, lower_offset, upper_offset) - } else { - continue; + // Process MS2 spectrum. + // Route through `make_isolation_window` for the f64 + // cvParam override path; see `convert_spectrum` and + // `read_isolation_cvparams_f64` for the rationale. + let isolation_window = match desc.precursor.first() { + Some(precursor) => { + match make_isolation_window(precursor, scan_number, Some(&iso_cv_params)) { + Some(w) => w, + None => continue, // No selected ion + } + } + None => continue, }; // Get peaks From 1089407fbb899c0db222962c3014d7950de48aa9 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 12:20:51 -0700 Subject: [PATCH 03/16] cal_match dump: bump fractional precision 10 -> 17 for f64 round-trip The cal_match diagnostic dump previously used `{:.10}` "everywhere so we don't hit banker's vs round-half-up rounding differences between Rust and C#". In practice the formatters still disagree at the 10th decimal on f64 values that land on a rounding boundary: Rust's `{:.10}` uses round-half-to-even while .NET Framework 4.7.2's `F10` uses round-half-away-from-zero, producing 1-in-the-last-digit (1e-10) print diffs even though the underlying f64s are bit-equal. Cross-impl strict comparison on Stellar 3-file 466K calibration matches: under `:.10` the apex_rt column showed 15,739 rows "differing" by exactly ~1e-10. Under `:.17` and parsed back to f64, every matched row's apex_rt is bit-equal cross-impl. The "drift" was 100% printf rounding artifact. 17 fractional digits is enough to round-trip any f64 uniquely, so the diagnostic dump now exposes only real f64-level divergence and not formatter noise. The remaining cal_match diffs (correlation, libcosine, xcorr, snr) are real f64-level drift from the f32-throughout xcorr preprocess on the Rust side; that is a separate concern. --- crates/osprey-scoring/src/diagnostics.rs | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/crates/osprey-scoring/src/diagnostics.rs b/crates/osprey-scoring/src/diagnostics.rs index 695bad9..653b4f3 100644 --- a/crates/osprey-scoring/src/diagnostics.rs +++ b/crates/osprey-scoring/src/diagnostics.rs @@ -53,11 +53,16 @@ pub fn dump_cal_match(library: &[LibraryEntry], results: &[CalibrationMatch]) { for entry in entries { if let Some(m) = by_id.get(&entry.id) { - // Use :.10 everywhere so we don't hit banker's vs round- - // half-up rounding differences between Rust and C#. + // Use :.17 (17 fractional digits) so f64 values round-trip + // exactly. Earlier :.10 was chosen "to avoid banker's vs + // round-half-up rounding differences between Rust and C#", + // but evidence shows the formatters still disagree at the + // 10th decimal when f64 values land near a rounding + // boundary. :.17 sidesteps the disagreement entirely: + // the underlying f64 bits are reproduced exactly. writeln!( f, - "{}\t{}\t{}\t1\t{}\t{:.10}\t{:.10}\t{:.10}\t{}\t{:.10}\t{:.10}", + "{}\t{}\t{}\t1\t{}\t{:.17}\t{:.17}\t{:.17}\t{}\t{:.17}\t{:.17}", entry.id, if entry.is_decoy { 1 } else { 0 }, entry.charge, From 3d2a0f5d6e6e43650e55d1e951ae2f4df67582e2 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 12:46:55 -0700 Subject: [PATCH 04/16] LDA scores dump: bump fractional precision 10 -> 17 for f64 round-trip MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Same rationale as the cal_match :.10 -> :.17 change in commit 1089407: Rust's `{:.10}` (round-half-to-even) and .NET Framework 4.7.2's `F10` (round-half-away-from-zero) disagree on f64 values that land on a 10th- decimal rounding boundary by 1 in the last digit, producing a fake 1e-10 "drift" between cross-impl dumps even when the underlying f64s are bit-equal. Cross-impl strict comparison on Stellar 3-file 466K calibration matches: under `:.10` the LDA discriminant/q-value showed 17 rows differing by 1e-10. Under `:.17` and parsed back to f64, max diff drops to 5.06e-14 (under 500 f64 ULP at value ~0.95) and most LDA scores are bit-equal modulo 1 ULP. The remaining LDA drift comes from each impl training LDA on a different per-entry calibration-match accumulator (Stage 3 cross-impl divergence in which apex spectrum is picked per entry) — separate concern. --- crates/osprey-scoring/src/diagnostics.rs | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/crates/osprey-scoring/src/diagnostics.rs b/crates/osprey-scoring/src/diagnostics.rs index 653b4f3..10035e8 100644 --- a/crates/osprey-scoring/src/diagnostics.rs +++ b/crates/osprey-scoring/src/diagnostics.rs @@ -102,8 +102,11 @@ pub fn dump_cal_match(library: &[LibraryEntry], results: &[CalibrationMatch]) { /// Dump per-entry LDA discriminant + q-value to `rust_lda_scores.txt`, /// sorted by entry_id for a stable diff against `cs_lda_scores.txt`. /// -/// Uses `{:.10}` formatting to avoid banker's-vs-half-up text rounding -/// mismatches with C#. +/// Uses `{:.17}` formatting so the underlying f64 round-trips exactly. +/// Earlier `{:.10}` was chosen "to avoid banker's vs half-up rounding +/// mismatches with C#"; in practice Rust `{:.10}` (RHE) and .NET +/// Framework `F10` (HAFZ) still disagree by 1 in the last digit on +/// boundary-case f64s. See the cal_match dump for the same fix. /// /// Gated by `OSPREY_DUMP_LDA_SCORES=1`. When `OSPREY_LDA_SCORES_ONLY=1` /// is also set, exits the process after writing. @@ -120,7 +123,7 @@ pub fn dump_lda_scores(matches: &[CalibrationMatch]) { let m = &matches[i]; writeln!( f, - "{}\t{}\t{:.10}\t{:.10}", + "{}\t{}\t{:.17}\t{:.17}", m.entry_id, if m.is_decoy { 1 } else { 0 }, m.discriminant_score, From 690194a9fb9a6522669f3bbfe4b68abfb7953db6 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 14:34:09 -0700 Subject: [PATCH 05/16] scorer.xcorr: inline f64 windowing + sliding-window for calibration parity Rewrites the body of SpectralScorer::xcorr to mirror C# OspreySharp's SpectralScorer.XcorrAtScan bit-for-bit: bin in f64 with `(intensity as f64).sqrt()` (widen-before-sqrt, matching `Math.Sqrt((double)float)`), apply windowing normalization and sliding- window subtraction in f64, sum unique fragment bins in f64, scale by `0.005_f64`. Removes the old allocating wrappers (`apply_windowing_- normalization` and `apply_sliding_window`); the `_into` variants stay for the HRAM per-window cache path. Calibration apex_rt, correlation, and libcosine columns of cal_match are now bit-equal cross-impl on Stellar Single (was 5.55e-16 / 4.97e-14 / 5.55e-16 noise, all at f64 epsilon). Cross-impl xcorr column drift is now 3.876e-6 (was 5.24e-10): not a regression. The previous baseline was both impls on f32, where both sides did pure-f32 windowing/sliding-window so the f32 cascade errors matched bit-for-bit. With Rust flipped to f64, C# calibration -- which still uses `XcorrFromPreprocessed(float[])` against a pure-f32 cache in PerFileScoringTask.cs:2275 -- diverges by exactly f32 magnitude. The coordinated fix requires C# calibration to switch to `XcorrAtScan` (existing f64 path) for the apex xcorr; tracked in TODO. HRAM main-search hot path (preprocess_spectrum_for_xcorr_into and the per-window f32 cache it feeds) is untouched: memory budget preserved. See ai/todos/active/TODO-20260516_ospreysharp_wsl_parity.md --- crates/osprey-scoring/src/lib.rs | 138 +++++++++++++++++++------------ 1 file changed, 86 insertions(+), 52 deletions(-) diff --git a/crates/osprey-scoring/src/lib.rs b/crates/osprey-scoring/src/lib.rs index fc2c184..19760d1 100644 --- a/crates/osprey-scoring/src/lib.rs +++ b/crates/osprey-scoring/src/lib.rs @@ -2081,43 +2081,95 @@ impl SpectralScorer { // First get LibCosine for additional metrics let lib_cosine_score = self.lib_cosine(observed, library); - // Bin observed spectrum using Comet BIN macro - let mut obs_binned = vec![0.0f32; self.bin_config.n_bins]; + // CALIBRATION PARITY: this path is f64 throughout to bit-equal C#'s + // OspreySharp.Scoring.SpectralScorer.XcorrAtScan, which is the + // calibration entry point on the C# side. The HRAM main-search hot + // path (preprocess_spectrum_for_xcorr_into) stays on the f32 cache + // for memory; only the scratch-internal math + this inline path go + // f64. Mismatching even one step (e.g. f32::sqrt vs widen-then-sqrt) + // amplifies through the windowing/sliding-window cascade to f32 + // magnitude (~1e-6). + let n_bins = self.bin_config.n_bins; + + // (1) Bin observed spectrum. Widen intensity to f64 BEFORE sqrt to + // match C#'s `Math.Sqrt((double)float)` bit-for-bit -- f32::sqrt and + // f64::sqrt(f32 as f64) can differ by 1 ULP at f32 magnitude. + let mut binned = vec![0.0f64; n_bins]; for (&mz, &intensity) in observed.mzs.iter().zip(observed.intensities.iter()) { if let Some(bin) = self.bin_config.mz_to_bin(mz) { - // Apply sqrt transformation to experimental spectrum - obs_binned[bin] += intensity.sqrt(); + binned[bin] += (intensity as f64).sqrt(); } } - // Apply windowing normalization (Comet's MakeCorrData) - let windowed = self.apply_windowing_normalization(&obs_binned); - - // Apply sliding window subtraction (fast XCorr preprocessing) - let xcorr_preprocessed = self.apply_sliding_window(&windowed); - - // XCorr = sum of preprocessed experimental values at UNIQUE fragment - // bin positions. When two library fragments fall into the same bin, - // the bin's contribution must count once, not twice -- the Comet - // theoretical spectrum uses unit intensity per bin (see - // preprocess_library_for_xcorr, which sets binned[bin] = 1.0 per - // unique bin). Summing preprocessed[bin] once per fragment instead - // of once per unique bin double-counts collisions and over-scores - // dense fragment lists. - let n_bins = xcorr_preprocessed.len(); + // (2) Windowing normalization (Comet MakeCorrData): 10 windows, + // normalize each window to max=50.0, zero values below 5% of global. + let mut windowed = vec![0.0f64; n_bins]; + { + let num_windows = 10; + let window_size = (n_bins / num_windows) + 1; + let global_max = binned.iter().cloned().fold(0.0f64, f64::max); + let threshold = global_max * 0.05; + for window_idx in 0..num_windows { + let start = window_idx * window_size; + let end = ((window_idx + 1) * window_size).min(n_bins); + let mut window_max = 0.0f64; + for &val in &binned[start..end] { + if val > window_max { + window_max = val; + } + } + if window_max > 0.0 { + let norm_factor = 50.0 / window_max; + for i in start..end { + if binned[i] > threshold { + windowed[i] = binned[i] * norm_factor; + } + } + } + } + } + + // (3) Sliding window subtraction (Comet fast XCorr): prefix-sum for + // O(n) window sums, divisor 2*offset = 150 regardless of boundary. + let mut xcorr_preprocessed = vec![0.0f64; n_bins]; + { + let offset: usize = 75; + let norm_factor = 1.0f64 / (2 * offset) as f64; + let mut prefix = vec![0.0f64; n_bins + 1]; + for i in 0..n_bins { + prefix[i + 1] = prefix[i] + windowed[i]; + } + for i in 0..n_bins { + let left = i.saturating_sub(offset); + let right = if i + offset < n_bins { + i + offset + 1 + } else { + n_bins + }; + let window_sum = prefix[right] - prefix[left]; + let sum_excluding_center = window_sum - windowed[i]; + xcorr_preprocessed[i] = windowed[i] - sum_excluding_center * norm_factor; + } + } + + // (4) Sum preprocessed values at UNIQUE library fragment bin + // positions. Comet theoretical spectrum uses unit intensity per + // bin, so collisions must count once (matches + // preprocess_library_for_xcorr which sets binned[bin] = 1.0). let mut visited = vec![false; n_bins]; let mut xcorr_raw: f64 = 0.0; for frag in &library.fragments { if let Some(bin) = self.bin_config.mz_to_bin(frag.mz) { if !visited[bin] { visited[bin] = true; - xcorr_raw += xcorr_preprocessed[bin] as f64; + xcorr_raw += xcorr_preprocessed[bin]; } } } - // Scale XCorr (pyXcorrDIA uses 0.005 for spectrum-centric) - let xcorr_scaled = xcorr_raw * 0.005; + // (5) Scale (pyXcorrDIA spectrum-centric). Explicit f64 literal + // matches C# `xcorrRaw * XCORR_SCALING` exactly. + let xcorr_scaled = xcorr_raw * 0.005_f64; SpectralScore { xcorr: xcorr_scaled, @@ -2281,19 +2333,12 @@ impl SpectralScorer { matches } - /// Apply Comet-style windowing normalization - /// - /// Divides spectrum into 10 windows and normalizes each to max=50.0 - fn apply_windowing_normalization(&self, spectrum: &[f32]) -> Vec { - let mut result = vec![0.0f32; spectrum.len()]; - self.apply_windowing_normalization_into(spectrum, &mut result); - result - } - - /// In-place variant of [`apply_windowing_normalization`] that writes - /// into a caller-provided buffer. The buffer must be zero on entry - /// because below-threshold and empty-window positions retain their - /// initial value. + /// In-place Comet-style windowing normalization (f32, used by the HRAM + /// per-window cache path). The buffer must be zero on entry because + /// below-threshold and empty-window positions retain their initial + /// value. The f64 calibration path inlines an equivalent computation + /// in [`SpectralScorer::xcorr`] for bit-exact alignment with the C# + /// OspreySharp port. fn apply_windowing_normalization_into(&self, spectrum: &[f32], result: &mut [f32]) { debug_assert_eq!( spectrum.len(), @@ -2331,23 +2376,12 @@ impl SpectralScorer { } } - /// Apply sliding window subtraction for fast XCorr (Comet-style) - /// - /// Uses prefix sum for O(n) performance instead of O(n × window). - /// Comet divides by (2*offset) = 150 regardless of boundary effects. - /// offset=75, matching Comet's iXcorrProcessingOffset default. - fn apply_sliding_window(&self, spectrum: &[f32]) -> Vec { - let n = spectrum.len(); - let mut prefix = vec![0.0f32; n + 1]; - let mut result = vec![0.0f32; n]; - self.apply_sliding_window_into(spectrum, &mut prefix, &mut result); - result - } - - /// In-place variant of [`apply_sliding_window`] that writes into - /// caller-provided prefix and result buffers. Every position of - /// `result` is overwritten; `prefix[0]` is assumed to be 0 on entry - /// and positions 1..=n are overwritten. + /// In-place sliding-window subtraction (Comet fast XCorr, f32, used by + /// the HRAM per-window cache path). Writes into caller-provided prefix + /// and result buffers. Every position of `result` is overwritten; + /// `prefix[0]` is assumed to be 0 on entry and positions 1..=n are + /// overwritten. The f64 calibration path inlines an equivalent + /// computation in [`SpectralScorer::xcorr`]. fn apply_sliding_window_into(&self, spectrum: &[f32], prefix: &mut [f32], result: &mut [f32]) { let n = spectrum.len(); let offset: usize = 75; From a44f7526d201feb7cba691bfd7b8855420bd4770 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 15:32:58 -0700 Subject: [PATCH 06/16] XCorr pipeline: f64 internal scratch + f32 storage cache Implements the "f64 scratch, f32 storage" architecture user-requested for Option 3 of the cross-impl calibration parity work. Both the calibration inline path (SpectralScorer::xcorr, surgical patch in 690194a) and the HRAM main-search cache build path now run the windowing / sliding-window cascade in f64 on per-thread pooled scratch, narrowing to f32 only at the final per-spectrum cache write. Cross-impl results on Stellar Single (joined by entry_id+charge+scan): * cal_match: apex_rt bit-equal, xcorr drift 3.876e-6 -> 4.429e-8 (~100x tighter). correlation, libcosine at f64 epsilon as before. * Stage 4 .scores.parquet: peak_apex 100% bit-equal, apex_rt 462,375/462,802 bit-equal, rt_deviation max 2.32e-11 (LOESS cascade now sees bit-equal inputs), xcorr max 5.41e-7 at the f32 single-cast architectural floor (was f32 cascade). Changes: * xcorr_pool.rs: XcorrScratch.binned/windowed/prefix Vec -> Vec. Per-thread scratch memory grows from ~1.2 MB to ~2.4 MB per worker on HRAM (~19 MB total at 16 threads). Per-spectrum cache (Vec>) unchanged at ~400 KB per spectrum. * lib.rs apply_windowing_normalization_into: f64 in, f64 out. Mirrors C# ApplyWindowingNormalizationD bit-for-bit. * lib.rs apply_sliding_window_into: f64 spectrum + f64 prefix scratch in, f32 result out. Single deterministic cast at final store. * lib.rs preprocess_spectrum_for_xcorr_into: bin with (intensity as f64).sqrt() (widen-before-sqrt, matches Math.Sqrt((double)float)). Output signature unchanged (&mut [f32]) so the HRAM per-window cache callers in pipeline.rs need no updates. * scorer.xcorr() inline body (from 690194a) remains the canonical reference for the f64 cascade; it does not yet share helpers with the new f64 _into path (cosmetic dedup deferred). See ai/todos/active/TODO-20260516_ospreysharp_wsl_parity.md --- crates/osprey-scoring/src/lib.rs | 70 ++++++++++++++----------- crates/osprey-scoring/src/xcorr_pool.rs | 40 ++++++++------ 2 files changed, 61 insertions(+), 49 deletions(-) diff --git a/crates/osprey-scoring/src/lib.rs b/crates/osprey-scoring/src/lib.rs index 19760d1..56a3bc2 100644 --- a/crates/osprey-scoring/src/lib.rs +++ b/crates/osprey-scoring/src/lib.rs @@ -2333,13 +2333,12 @@ impl SpectralScorer { matches } - /// In-place Comet-style windowing normalization (f32, used by the HRAM - /// per-window cache path). The buffer must be zero on entry because + /// In-place Comet-style windowing normalization (f64 internal, + /// dual-purpose helper for both calibration and the HRAM per-window + /// cache path). The buffer must be zero on entry because /// below-threshold and empty-window positions retain their initial - /// value. The f64 calibration path inlines an equivalent computation - /// in [`SpectralScorer::xcorr`] for bit-exact alignment with the C# - /// OspreySharp port. - fn apply_windowing_normalization_into(&self, spectrum: &[f32], result: &mut [f32]) { + /// value. Matches C# `ApplyWindowingNormalizationD` bit-for-bit. + fn apply_windowing_normalization_into(&self, spectrum: &[f64], result: &mut [f64]) { debug_assert_eq!( spectrum.len(), result.len(), @@ -2348,23 +2347,20 @@ impl SpectralScorer { let num_windows = 10; let window_size = (spectrum.len() / num_windows) + 1; - // Find global max for threshold - let global_max = spectrum.iter().cloned().fold(0.0f32, f32::max); + let global_max = spectrum.iter().cloned().fold(0.0f64, f64::max); let threshold = global_max * 0.05; for window_idx in 0..num_windows { let start = window_idx * window_size; let end = ((window_idx + 1) * window_size).min(spectrum.len()); - // Find max in this window - let mut window_max = 0.0f32; + let mut window_max = 0.0f64; for &val in &spectrum[start..end] { if val > window_max { window_max = val; } } - // Normalize this window to 50.0 if window_max > 0.0 { let norm_factor = 50.0 / window_max; for i in start..end { @@ -2376,23 +2372,24 @@ impl SpectralScorer { } } - /// In-place sliding-window subtraction (Comet fast XCorr, f32, used by - /// the HRAM per-window cache path). Writes into caller-provided prefix - /// and result buffers. Every position of `result` is overwritten; - /// `prefix[0]` is assumed to be 0 on entry and positions 1..=n are - /// overwritten. The f64 calibration path inlines an equivalent - /// computation in [`SpectralScorer::xcorr`]. - fn apply_sliding_window_into(&self, spectrum: &[f32], prefix: &mut [f32], result: &mut [f32]) { + /// In-place sliding-window subtraction (Comet fast XCorr) with f64 + /// internal math narrowing to an f32 result buffer at the final store. + /// `spectrum` and `prefix` are f64 scratch; `result` is the per-spectrum + /// f32 cache. Every position of `result` is overwritten; `prefix[0]` is + /// assumed to be 0 on entry and positions 1..=n are overwritten. + /// Matches C# `ApplySlidingWindowD` for the f64 cascade, with a single + /// final cast to f32 on store. Used by the HRAM per-window cache path. + fn apply_sliding_window_into(&self, spectrum: &[f64], prefix: &mut [f64], result: &mut [f32]) { let n = spectrum.len(); let offset: usize = 75; - // Comet uses (window_size - 1) = 2*offset = 150 as divisor - let norm_factor = 1.0f32 / (2 * offset) as f32; + // Comet uses (window_size - 1) = 2*offset = 150 as divisor. + // f64 literal so the final subtraction stays in f64. + let norm_factor = 1.0f64 / (2 * offset) as f64; debug_assert_eq!(prefix.len(), n + 1, "prefix buffer must be n+1"); debug_assert_eq!(result.len(), n, "result buffer must be n"); debug_assert_eq!(prefix[0], 0.0, "prefix[0] must be 0 on entry"); - // Build prefix sum for O(n) window sums for i in 0..n { prefix[i + 1] = prefix[i] + spectrum[i]; } @@ -2400,12 +2397,12 @@ impl SpectralScorer { for i in 0..n { let left = i.saturating_sub(offset); let right = if i + offset < n { i + offset + 1 } else { n }; - // Window sum including center let window_sum = prefix[right] - prefix[left]; - // Subtract center to get sum excluding center let sum_excluding_center = window_sum - spectrum[i]; - // Subtract local average from center value - result[i] = spectrum[i] - sum_excluding_center * norm_factor; + let centered = spectrum[i] - sum_excluding_center * norm_factor; + // f64 -> f32 narrowing at the final store: single deterministic + // rounding, identical bits cross-impl when the f64 inputs agree. + result[i] = centered as f32; } } @@ -2427,11 +2424,17 @@ impl SpectralScorer { } /// In-place XCorr preprocess using a pool-rented scratch and a - /// caller-provided output buffer. The output buffer is fully - /// overwritten. The scratch's accumulator fields (`binned`, - /// `windowed`) are zeroed at the start of the call, so callers can - /// reuse a single rented scratch across many preprocess calls - /// without interleaving recycle/rent. + /// caller-provided f32 output buffer (the per-spectrum HRAM cache). + /// The scratch fields are f64 so the windowing / sliding-window + /// cascade runs in f64 precision; only the final cache write narrows + /// to f32. Matches C# OspreySharp's PreprocessSpectrumForXcorrF32 with + /// f64-internal math: same f32 cache size, but values carry single- + /// cast precision rather than f32-cascade noise. + /// + /// The scratch's accumulator fields (`binned`, `windowed`) are zeroed + /// at the start of the call, so callers can reuse a single rented + /// scratch across many preprocess calls without interleaving + /// recycle/rent. pub fn preprocess_spectrum_for_xcorr_into( &self, spectrum: &Spectrum, @@ -2451,10 +2454,13 @@ impl SpectralScorer { scratch.binned.fill(0.0); scratch.windowed.fill(0.0); - // Bin observed spectrum with sqrt transformation using Comet BIN macro + // Bin observed spectrum with sqrt transformation using Comet BIN + // macro. Widen f32 intensity to f64 BEFORE sqrt to match C#'s + // `Math.Sqrt((double)float)` bit-for-bit; f32::sqrt and + // f64::sqrt(f32 as f64) can differ by 1 ULP at f32 magnitude. for (&mz, &intensity) in spectrum.mzs.iter().zip(spectrum.intensities.iter()) { if let Some(bin) = self.bin_config.mz_to_bin(mz) { - scratch.binned[bin] += intensity.sqrt(); + scratch.binned[bin] += (intensity as f64).sqrt(); } } diff --git a/crates/osprey-scoring/src/xcorr_pool.rs b/crates/osprey-scoring/src/xcorr_pool.rs index 0e62949..e0f61c3 100644 --- a/crates/osprey-scoring/src/xcorr_pool.rs +++ b/crates/osprey-scoring/src/xcorr_pool.rs @@ -1,16 +1,22 @@ //! Buffer pool for XCorr preprocessing scratch space. //! -//! XCorr preprocessing allocates three `NBins`-sized `f32` buffers per call -//! (binned accumulator, windowed result, sliding-window prefix sum). On -//! HRAM binning (NBins ~= 100K) each buffer is ~400 KB, so per-window -//! scoring across ~1000 spectra churns through ~1.2 GB of transient -//! `Vec` allocations plus the ~400 MB per-window preprocessed cache -//! itself. All of this memory is dropped at the end of each window and -//! re-allocated by the next window -- pure churn. +//! XCorr preprocessing allocates three `NBins`-sized `f64` scratch buffers +//! per call (binned accumulator, windowed result, sliding-window prefix +//! sum) and writes the final preprocessed values into a per-spectrum `f32` +//! cache. The f64 scratch is dual-purpose: it keeps the windowing / +//! sliding-window cascade math in f64 precision (avoiding the ~1e-5 +//! prefix-sum error f32 would accumulate over ~100K bins), and only the +//! final write to the cache narrows to f32. The per-spectrum cache stays +//! f32 to keep HRAM memory budget at ~400 KB per spectrum. //! -//! This pool hands out reusable scratch bundles and single output buffers -//! so that after the first few windows reach the high-water mark the hot -//! path never allocates. Mirrors the C# OspreySharp `XcorrScratchPool` in +//! On HRAM binning (NBins ~= 100K) each scratch buffer is ~800 KB, so per- +//! window scoring across ~1000 spectra would otherwise churn through +//! ~2.4 GB of transient `Vec` allocations plus the ~400 MB per-window +//! preprocessed cache. This pool hands out reusable scratch bundles and +//! single output buffers so that after the first few windows reach the +//! high-water mark the hot path never allocates. +//! +//! Mirrors the C# OspreySharp `XcorrScratchPool` in //! `pwiz_tools/OspreySharp/OspreySharp.Scoring/XcorrScratchPool.cs`. //! //! Thread-safety: the pool uses `Mutex>` bags. Contention is low @@ -24,16 +30,16 @@ use std::sync::Mutex; /// Rented by value, returned by value via [`XcorrScratchPool::recycle`]. pub struct XcorrScratch { /// `NBins` accumulator for binned-and-sqrt-transformed spectrum - /// intensities. Written via `+=` so must start zeroed; zeroed on - /// recycle. - pub binned: Vec, - /// `NBins` windowing-normalization output. Filled by + /// intensities, in f64 so the windowing pipeline preserves precision. + /// Written via `+=` so must start zeroed; zeroed on recycle. + pub binned: Vec, + /// `NBins` windowing-normalization output, f64. Filled by /// `apply_windowing_normalization_into`; zeroed on recycle so /// below-threshold positions retain 0. - pub windowed: Vec, - /// `NBins + 1` sliding-window prefix sum. `prefix[0]` is always 0; + pub windowed: Vec, + /// `NBins + 1` sliding-window prefix sum, f64. `prefix[0]` is always 0; /// positions `1..=NBins` are fully overwritten each call. - pub prefix: Vec, + pub prefix: Vec, } impl XcorrScratch { From f41ff88876bea6a4978c605a046a37ece43f0ab7 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 15:55:08 -0700 Subject: [PATCH 07/16] XCorr alignment: f64 accumulator + scorer.xcorr via cache narrowing Closes the residual ~1e-7 cross-impl xcorr drift after Option 3 by matching C#'s accumulator type and code path: * xcorr_sparse: f32 accumulator -> f64 accumulator (widen each cached value to f64 on read, sum in f64, scale by 0.005_f64). Mirrors C# XcorrFromPreprocessed(float[]) which does `double xcorrRaw += preprocessed[bin]` (implicit float->double promote on read). * scorer.xcorr(): the inline f64 preprocessing body is replaced by delegation to xcorr_at_scan, which goes through the f32 cache (preprocess_spectrum_for_xcorr returns the same f64-internal-narrowed Vec the HRAM main-search cache uses). Now bit-equal with C# calibration's XcorrFromPreprocessed(windowPreprocessedF32, entry) call site at PerFileScoringTask.cs:2275 because both impls do the same operations on the same f32 cache values. Cross-impl results on Stellar Single (joined by entry_id+charge+scan): * cal_match xcorr: 4.43e-8 -> 5.11e-15 (f64 epsilon, was f32 cast floor) * cal_match LDA scores: 3.29e-5 -> 4.95e-14 (f64 epsilon, cascade from bit-equal LDA inputs) * Stage 4 .scores.parquet xcorr: 100% bit-equal (462801/462801) * Stage 4 .scores.parquet sg_weighted_xcorr: 100% bit-equal * Stage 4 .scores.parquet peak_apex: 100% bit-equal * Stage 4 .scores.parquet apex_rt: 462375 bit-equal + f64 epsilon * rt_deviation: max 2.32e-11 (LOESS cascade, bulk <1e-12) Only cross-impl drift remaining in cal_match is snr (5.24e-10), which is the LDA in-place mutation hypothesis - separate root cause, next. --- crates/osprey-scoring/src/lib.rs | 117 ++++++------------------------- 1 file changed, 20 insertions(+), 97 deletions(-) diff --git a/crates/osprey-scoring/src/lib.rs b/crates/osprey-scoring/src/lib.rs index 56a3bc2..9da5b9f 100644 --- a/crates/osprey-scoring/src/lib.rs +++ b/crates/osprey-scoring/src/lib.rs @@ -2081,95 +2081,14 @@ impl SpectralScorer { // First get LibCosine for additional metrics let lib_cosine_score = self.lib_cosine(observed, library); - // CALIBRATION PARITY: this path is f64 throughout to bit-equal C#'s - // OspreySharp.Scoring.SpectralScorer.XcorrAtScan, which is the - // calibration entry point on the C# side. The HRAM main-search hot - // path (preprocess_spectrum_for_xcorr_into) stays on the f32 cache - // for memory; only the scratch-internal math + this inline path go - // f64. Mismatching even one step (e.g. f32::sqrt vs widen-then-sqrt) - // amplifies through the windowing/sliding-window cascade to f32 - // magnitude (~1e-6). - let n_bins = self.bin_config.n_bins; - - // (1) Bin observed spectrum. Widen intensity to f64 BEFORE sqrt to - // match C#'s `Math.Sqrt((double)float)` bit-for-bit -- f32::sqrt and - // f64::sqrt(f32 as f64) can differ by 1 ULP at f32 magnitude. - let mut binned = vec![0.0f64; n_bins]; - for (&mz, &intensity) in observed.mzs.iter().zip(observed.intensities.iter()) { - if let Some(bin) = self.bin_config.mz_to_bin(mz) { - binned[bin] += (intensity as f64).sqrt(); - } - } - - // (2) Windowing normalization (Comet MakeCorrData): 10 windows, - // normalize each window to max=50.0, zero values below 5% of global. - let mut windowed = vec![0.0f64; n_bins]; - { - let num_windows = 10; - let window_size = (n_bins / num_windows) + 1; - let global_max = binned.iter().cloned().fold(0.0f64, f64::max); - let threshold = global_max * 0.05; - for window_idx in 0..num_windows { - let start = window_idx * window_size; - let end = ((window_idx + 1) * window_size).min(n_bins); - let mut window_max = 0.0f64; - for &val in &binned[start..end] { - if val > window_max { - window_max = val; - } - } - if window_max > 0.0 { - let norm_factor = 50.0 / window_max; - for i in start..end { - if binned[i] > threshold { - windowed[i] = binned[i] * norm_factor; - } - } - } - } - } - - // (3) Sliding window subtraction (Comet fast XCorr): prefix-sum for - // O(n) window sums, divisor 2*offset = 150 regardless of boundary. - let mut xcorr_preprocessed = vec![0.0f64; n_bins]; - { - let offset: usize = 75; - let norm_factor = 1.0f64 / (2 * offset) as f64; - let mut prefix = vec![0.0f64; n_bins + 1]; - for i in 0..n_bins { - prefix[i + 1] = prefix[i] + windowed[i]; - } - for i in 0..n_bins { - let left = i.saturating_sub(offset); - let right = if i + offset < n_bins { - i + offset + 1 - } else { - n_bins - }; - let window_sum = prefix[right] - prefix[left]; - let sum_excluding_center = window_sum - windowed[i]; - xcorr_preprocessed[i] = windowed[i] - sum_excluding_center * norm_factor; - } - } - - // (4) Sum preprocessed values at UNIQUE library fragment bin - // positions. Comet theoretical spectrum uses unit intensity per - // bin, so collisions must count once (matches - // preprocess_library_for_xcorr which sets binned[bin] = 1.0). - let mut visited = vec![false; n_bins]; - let mut xcorr_raw: f64 = 0.0; - for frag in &library.fragments { - if let Some(bin) = self.bin_config.mz_to_bin(frag.mz) { - if !visited[bin] { - visited[bin] = true; - xcorr_raw += xcorr_preprocessed[bin]; - } - } - } - - // (5) Scale (pyXcorrDIA spectrum-centric). Explicit f64 literal - // matches C# `xcorrRaw * XCORR_SCALING` exactly. - let xcorr_scaled = xcorr_raw * 0.005_f64; + // XCorr via the same cache-narrowing path C# OspreySharp's + // calibration uses (XcorrFromPreprocessed(float[] cache, entry)): + // f64-internal preprocessing narrows to an f32 cache, then the + // sparse sum widens each cached value back to f64 on read and + // accumulates in f64. Bit-equal cross-impl as long as the f32 + // cache is bit-equal (Option 3 guarantees this) and both impls + // promote-on-read into f64 accumulators. + let xcorr_scaled = self.xcorr_at_scan(observed, library); SpectralScore { xcorr: xcorr_scaled, @@ -2507,11 +2426,15 @@ impl SpectralScorer { /// zero*something. Iterating fragments directly is O(n_frags) and /// allocates only a small stack of bin indices for deduplication. /// - /// Bit-parity with the dense path: the final scale and cast follow the - /// exact pattern `(sum_f32 * 0.005_f32) as f64` so both return the same - /// numerical value; the only difference is the summation ORDER (linear - /// left-to-right vs SIMD reduction), which typically differs by less - /// than 1e-7 on xcorr values around 1.0. + /// Cross-impl alignment: each f32 cache value is widened to f64 on read + /// and accumulated in f64; final scale uses an f64 literal. Matches + /// C# `SpectralScorer.XcorrFromPreprocessed(float[], LibraryEntry)` + /// which sums `xcorrRaw += preprocessed[bin]` into a `double` accumulator + /// (implicit float-to-double promotion on read). The iteration order + /// differs (Rust sorts bins ascending; C# iterates fragments in input + /// order with a visited[bin] dedup) but at f64 precision the resulting + /// sum-order difference is bounded by ~f64_eps × intermediate magnitude + /// (~1e-13 for typical xcorr) — well below f64 epsilon noise. #[inline] pub fn xcorr_sparse(&self, spectrum_preprocessed: &[f32], entry: &LibraryEntry) -> f64 { // Small stack-friendly accumulator for unique fragment bins. At HRAM @@ -2528,17 +2451,17 @@ impl SpectralScorer { bins.dedup(); let n_bins = spectrum_preprocessed.len(); - let sum: f32 = bins + let sum: f64 = bins .iter() .filter_map(|&b| { if b < n_bins { - Some(spectrum_preprocessed[b]) + Some(spectrum_preprocessed[b] as f64) } else { None } }) .sum(); - (sum * 0.005f32) as f64 + sum * 0.005_f64 } /// Preprocess a library entry for XCorr (Comet-style) From b80d86b5c2ec5dba37584b982ac23249c7c3c26c Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 16:13:30 -0700 Subject: [PATCH 08/16] Calibration pass 2: refresh LOESS dump + num_confident_peptides metadata MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When the second-pass calibration refinement is accepted (R² check at pipeline.rs:1105), Rust now also: * Overwrites the LOESS_INPUT diagnostic dump with pass 2's points. Previously dump_loess_input was only called at line 1000 (before the refinement block) so the file reflected pass 1's 6,398 points even when pass 2's 7,361 points were what the LOESS fit actually used. C# overwrites the dump unconditionally on pass 2; this brings Rust into parity. Stellar Single now passes the LOESS_INPUT boundary bit-equal (6400 -> 7361 rows on the dump, matching the actual calibration data). * Updates num_confident_peptides metadata in calibration.json to the refined count. Was reporting pass 1's value (6,398) even when pass 2's 7,361 was the count actually used for the fit. Both changes are bookkeeping/observability fixes: the rt_calibration model parameters themselves were already bit-equal cross-impl at f64 epsilon (~1e-13) because both impls used the same 7,361 pass-2 points for the LOESS fit. Only the diagnostic dump and the metadata count were stuck on pass 1. Remaining cross-impl divergence in calibration.json is ms1_calibration {count, mean, median, sd, adjusted_tolerance} — Rust 18 errors vs C# 193. Separate root cause (ms1_error collection path). --- crates/osprey/src/pipeline.rs | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index 285791b..db5c862 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -896,7 +896,7 @@ fn run_calibration_discovery_windowed( } } - let num_confident_peptides = library_rts_detected.len(); + let mut num_confident_peptides = library_rts_detected.len(); // Compute median peak width from confident matches for adaptive co-elution window // Log median peak width from calibration matches (diagnostic) @@ -1105,6 +1105,17 @@ fn run_calibration_discovery_windowed( if rt_stats_refined.r_squared >= rt_stats.r_squared * 0.99 { rt_calibration = rt_cal_refined; rt_stats = rt_stats_refined; + // Update metadata to reflect the refined fit that's + // actually being used (was pass 1's count). C# parity: + // PerFileScoringTask.cs reports n_refined when pass 2 is + // accepted. + num_confident_peptides = n_refined; + // Overwrite the LOESS input dump with pass 2's points + // so the diagnostic reflects the calibration actually + // used. C# overwrites unconditionally on pass 2; Rust + // previously only wrote pass 1 (line 1000), leaving + // ~960 entries unreported here on Stellar Single. + crate::diagnostics::dump_loess_input(&refined_lib_rts, &refined_meas_rts); // Re-collect mass errors from refined matches mz_qc_data = MzQCData::new(config.fragment_tolerance.unit); From 99c9e304a81e46a78701379ae432092ee4ef3a11 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 16:40:58 -0700 Subject: [PATCH 09/16] Treat non-ppm precursor tolerance as 10 ppm default for MS1 envelope The MS1 isotope envelope extraction at calibration time (batch.rs:2708 in run_coelution_calibration_scoring) takes a tolerance_ppm argument and uses it for FindPeakPpm-style matching of the M+0 peak. Rust was passing config.precursor_tolerance.tolerance regardless of unit. For unit-resolution data (Stellar config: tolerance=1.0, unit=Mz), this treated 1.0 Da as 1.0 ppm and produced a ~0.5 mDa window at 500 m/z - effectively zero. As a result envelope.has_m0() returned false on ~99.8% of matches and only 18 of 7361 calibration-passing entries contributed to ms1_calibration on Stellar Single, vs 193 on C# (10x undercount). C# OspreySharp's PerFileScoringTask.ScoreCalibrationEntry handles the same case via `unit == Ppm ? tolerance : 10.0`. This adds a Rust helper `ms1_envelope_tolerance_ppm` with identical behavior and wires it into both pass-1 and pass-2 run_coelution_calibration_scoring callsites in pipeline.rs. Cross-impl effect on Stellar Single calibration.json: ms1_calibration.{count, mean, median, sd, adjusted_tolerance} drift disappears from the diff (was the dominant CAL_JSON divergence with diff 1.75e+2 at ms1_calibration.count). Remaining cal_json drift is rt_calibration.model_params.abs_residuals[i] sort-order swaps on a handful of indices, where the underlying values are bit-equal but sorted slightly differently cross-impl. Found via OSPREY_DIAG_MS1 instrumentation showing tol_ppm=1.0 and peaks_in_m0_window=0 on all sampled calibration entries. --- crates/osprey/src/pipeline.rs | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index db5c862..7b6644e 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -62,6 +62,22 @@ use osprey_scoring::{ DecoyGenerator, DecoyMethod, Enzyme, SpectralScorer, }; +/// Resolve the precursor tolerance to a ppm value usable for MS1 isotope +/// envelope extraction. If the configured precursor tolerance is already +/// in ppm, return it as-is. If it's a Da/Mz tolerance (typical for unit- +/// resolution data like Stellar), treating it as ppm would produce an +/// absurdly tight window (e.g. 1 Da -> 0.0005 mDa at 500 m/z), making +/// `envelope.has_m0()` fail on ~99.8% of matches. Fall back to a 10 ppm +/// default in that case, matching C# OspreySharp's PerFileScoringTask +/// ScoreCalibrationEntry behavior: +/// `unit == Ppm ? tolerance : 10.0`. +fn ms1_envelope_tolerance_ppm(precursor_tolerance: &osprey_core::FragmentToleranceConfig) -> f64 { + match precursor_tolerance.unit { + osprey_core::ToleranceUnit::Ppm => precursor_tolerance.tolerance, + osprey_core::ToleranceUnit::Mz => 10.0, + } +} + /// Wrapper to implement MS1SpectrumLookup for MS1Index struct MS1IndexWrapper<'a>(&'a MS1Index); @@ -746,7 +762,7 @@ fn run_calibration_discovery_windowed( spectra, Some(&MS1IndexWrapper(ms1_index)), config.fragment_tolerance, - config.precursor_tolerance.tolerance, + ms1_envelope_tolerance_ppm(&config.precursor_tolerance), initial_tolerance, pass1_expected_rt_fn, Some(&xcorr_scorer), @@ -758,7 +774,7 @@ fn run_calibration_discovery_windowed( spectra, None, config.fragment_tolerance, - config.precursor_tolerance.tolerance, + ms1_envelope_tolerance_ppm(&config.precursor_tolerance), initial_tolerance, pass1_expected_rt_fn, Some(&xcorr_scorer), @@ -1037,7 +1053,7 @@ fn run_calibration_discovery_windowed( spectra, Some(&MS1IndexWrapper(ms1_index)), config.fragment_tolerance, - config.precursor_tolerance.tolerance, + ms1_envelope_tolerance_ppm(&config.precursor_tolerance), pass1_tolerance, Some(&predict_fn), Some(&xcorr_scorer), @@ -1049,7 +1065,7 @@ fn run_calibration_discovery_windowed( spectra, None, config.fragment_tolerance, - config.precursor_tolerance.tolerance, + ms1_envelope_tolerance_ppm(&config.precursor_tolerance), pass1_tolerance, Some(&predict_fn), Some(&xcorr_scorer), From adfbea42ae2e819ab4c9ed588e62eb15e4ebf71b Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 16:50:34 -0700 Subject: [PATCH 10/16] LOESS: sort by (lib_rt, measured_rt) tuple for deterministic dup-x order Stable sort_by on x alone preserves input order for ties. When inputs arrive in different orders cross-impl (e.g. discriminant-score-sorted with 1 ULP differences), the duplicate-x positions end up at different indices and the LOESS abs_residuals[i] diverge cross-impl. Sorting by (x, y) makes the order deterministic and removes the input-order dependency. C# OspreySharp LoessRegression.cs has the matching fix (LINQ ThenBy on y). --- crates/osprey-chromatography/src/calibration/rt.rs | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/crates/osprey-chromatography/src/calibration/rt.rs b/crates/osprey-chromatography/src/calibration/rt.rs index e77b544..5a05fae 100644 --- a/crates/osprey-chromatography/src/calibration/rt.rs +++ b/crates/osprey-chromatography/src/calibration/rt.rs @@ -114,13 +114,18 @@ impl RTCalibrator { ))); } - // Sort by library RT + // Sort by library RT, breaking ties on measured RT so the order is + // deterministic for duplicate-x inputs. Without the secondary key, + // pairs with equal library_rt land in input-order on Rust (stable + // sort_by) but may land in a different order cross-impl, producing + // cosmetic swaps in abs_residuals[i] downstream that the cross-impl + // calibration.json diff catches. let mut pairs: Vec<(f64, f64)> = library_rts .iter() .zip(measured_rts.iter()) .map(|(&x, &y)| (x, y)) .collect(); - pairs.sort_by(|a, b| a.0.total_cmp(&b.0)); + pairs.sort_by(|a, b| a.0.total_cmp(&b.0).then_with(|| a.1.total_cmp(&b.1))); let x: Vec = pairs.iter().map(|(x, _)| *x).collect(); let y: Vec = pairs.iter().map(|(_, y)| *y).collect(); From 8dc0441b0762608160c8bd304d520ca46a2d8fbf Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Wed, 20 May 2026 23:51:24 -0700 Subject: [PATCH 11/16] Cross-impl bit-equality: align cosine + sorted parquet write Two small alignments bring Stage 4 .scores.parquet and Stage 5/6 dumps to byte-identical cross-impl on Stellar net8.0 in the per-side OspreySharp parity harness: * compute_cosine_at_scan: switch to single-pass dot+norms with one final divide (matches the cosine_angle helper in the same file and the C# port's ComputeCosineAtScan), bringing sg_weighted_cosine to bit-equal cross-impl. The previous per-element divide-then-sum form is mathema- tically equivalent but differs in the last 1-2 ULP. * write_scores_parquet_with_metadata: iterate entries in canonical (entry_id, charge, scan_number) order before filling the column builders. Physical row order in the parquet was previously whatever upstream deduplication left it in, and Rust's entry_id ascending pattern did not match the C# port's target-decoy-paired pattern. The Stage 5 standardizer and the SVM working-set selection sum / iterate in physical row order, so per-side cross-impl runs hit a sum-order cascade in the standardizer even though every column was logically bit-equal. Sorting before write isolates Stage 5+ from parquet-writer row-order differences. Both changes are pure no-op on a single side reading itself back; only the cross-impl behaviour changes. cargo fmt / clippy / test all pass. --- crates/osprey-scoring/src/lib.rs | 23 ++++++++++++++--------- crates/osprey/src/pipeline.rs | 17 ++++++++++++++++- 2 files changed, 30 insertions(+), 10 deletions(-) diff --git a/crates/osprey-scoring/src/lib.rs b/crates/osprey-scoring/src/lib.rs index 9da5b9f..2cfeead 100644 --- a/crates/osprey-scoring/src/lib.rs +++ b/crates/osprey-scoring/src/lib.rs @@ -430,20 +430,25 @@ pub fn compute_cosine_at_scan( return 0.0; } - // L2 normalize - let lib_norm = lib_preprocessed.iter().map(|x| x * x).sum::().sqrt(); - let obs_norm = obs_preprocessed.iter().map(|x| x * x).sum::().sqrt(); + // Single-pass accumulation of squared norms + dot product, then divide once + // at the end. Matches `cosine_angle` in this crate and the C# port's + // ComputeCosineAtScan, ensuring cross-impl bit-equality for sg_weighted_cosine. + let mut lib_sq = 0.0; + let mut obs_sq = 0.0; + let mut dot = 0.0; + for (a, b) in lib_preprocessed.iter().zip(obs_preprocessed.iter()) { + lib_sq += a * a; + obs_sq += b * b; + dot += a * b; + } + let lib_norm = lib_sq.sqrt(); + let obs_norm = obs_sq.sqrt(); if lib_norm < 1e-10 || obs_norm < 1e-10 { return 0.0; } - // Cosine = dot(lib_norm, obs_norm) - lib_preprocessed - .iter() - .zip(obs_preprocessed.iter()) - .map(|(a, b)| (a / lib_norm) * (b / obs_norm)) - .sum() + dot / (lib_norm * obs_norm) } /// Compute per-fragment mass accuracy statistics at the apex scan. diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index 7b6644e..f5e97e9 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -1763,7 +1763,22 @@ fn write_scores_parquet_with_metadata( .map(|_| Float64Builder::with_capacity(n)) .collect(); - for entry in entries { + // Iterate in canonical sorted order (entry_id, charge, scan_number) so + // per-side parquets have identical physical row layout across Rust and C# + // impls. Order-sensitive consumers downstream (Stage 5 standardizer, + // SVM training) then see the same row sequence regardless of which side + // wrote the parquet, making per-side cross-tool stages bit-equal. + let mut sorted_indices: Vec = (0..entries.len()).collect(); + sorted_indices.sort_by(|&a, &b| { + entries[a] + .entry_id + .cmp(&entries[b].entry_id) + .then_with(|| entries[a].charge.cmp(&entries[b].charge)) + .then_with(|| entries[a].scan_number.cmp(&entries[b].scan_number)) + }); + + for &idx in &sorted_indices { + let entry = &entries[idx]; entry_id_b.append_value(entry.entry_id); decoy_b.append_value(entry.is_decoy); seq_b.append_value(&entry.sequence); From 566f5835113b3894810ffed8f3a8f9068c356e4c Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Thu, 21 May 2026 00:11:27 -0700 Subject: [PATCH 12/16] Cross-impl bit-equality: always persist 2nd-pass FDR sidecar MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The 2nd-pass `.fdr_scores.bin` was previously written only when reconciliation was enabled across multiple files. In single-file mode that gated the sidecar off, while the OspreySharp port writes it unconditionally — so per-side cross-impl Stage 7 (`--join-at-pass=2`) hit a load-cached-vs-retrain asymmetry: C# loaded the cached 2nd-pass scores and Rust re-trained the 2nd-pass SVM from scratch, producing different SVM weights and therefore different downstream q-values. With the gate removed, both sides write a 2nd-pass sidecar containing the post-compaction FDR scores (which in single-file mode equal the 1st-pass scores, since no rescore happens). The cross-impl sidecars are then bit-identical (verified on Stellar Single, net8.0) and the load path is symmetric. This closes the Stage 6 → Stage 7 boundary as a source of cross-impl drift; any remaining Stage 7 divergence is in the protein-FDR / parsimony code itself, not in upstream inputs. Single-file workflows gain ~5-9 MB of extra disk per file for the new sidecar; the file is ignored by all current consumers except the resume path that wanted it anyway. cargo fmt / clippy / test all pass. --- crates/osprey/src/pipeline.rs | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index f5e97e9..51d47b2 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -4996,22 +4996,27 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { } } - // Persist second-pass SVM scores to sidecar files (after reconciliation block closes) + // Persist second-pass SVM scores to sidecar files (after reconciliation block closes). + // Always written when FDR ran, regardless of reconciliation. In single-file + // mode there is no rescore step, so the persisted scores equal the 1st-pass + // scores; this is fine semantically (per_file_entries already carry the + // authoritative final-pass scores) and required for cross-impl parity with + // the C# port, whose Stage 7 (--join-at-pass=2) reads the 2nd-pass sidecar + // unconditionally. Without an unconditional write here, single-file per-side + // cross-impl runs hit a load-cached-vs-retrain asymmetry in Stage 7 protein + // FDR even when every upstream stage is bit-equal. if !can_skip_fdr { - let has_reconciliation = config.reconciliation.enabled && config.input_files.len() > 1; - if has_reconciliation { - log::debug!("Persisting 2nd-pass FDR scores..."); - // 2nd-pass is a resume-only optimization (skip-Percolator on - // reruns); the return value is logged via per-file warnings - // already, so no need to escalate here. - let _ = persist_fdr_scores( - &per_file_entries, - &config, - fdr_scores_path_pass2, - "2nd-pass", - 2, - ); - } + log::debug!("Persisting 2nd-pass FDR scores..."); + // 2nd-pass is a resume-only optimization (skip-Percolator on + // reruns); the return value is logged via per-file warnings + // already, so no need to escalate here. + let _ = persist_fdr_scores( + &per_file_entries, + &config, + fdr_scores_path_pass2, + "2nd-pass", + 2, + ); } // Protein parsimony (always runs) + optional second-pass picked-protein FDR. From 84ff72cba5bb5c238c95e0f18bfc6b217ea23e5a Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Thu, 21 May 2026 06:42:08 -0700 Subject: [PATCH 13/16] Add Stage 7 detected_peptides bisection dump MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cross-impl bisection of the C# OspreySharp port revealed that its --join-at-pass=2 path was filtering second-pass detected_peptides through stale 1st-pass q-values rather than reloading the 2nd-pass FDR sidecar onto post-compaction stubs (Rust does this at the pipeline.rs:4480-4494 reload block). The discrepancy was 19 peptides on Stellar Single — bordering peptides that pass 1st-pass FDR at <=1% but not 2nd-pass — and produced a 1-protein delta in the Stage 7 picked-protein output. The dump is gated on OSPREY_DUMP_DETECTED_PEPTIDES=1 and writes rust_stage7_detected_peptides.txt next to the working directory. Zero overhead when unset. Matches the cs_stage7_detected_peptides.txt dump the C# port writes under the same env var, so a sorted-line diff localizes any remaining detected_peptides drift to the specific modified_sequence values. --- crates/osprey/src/pipeline.rs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index 51d47b2..71d1fa8 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -5054,6 +5054,21 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { .map(|e| e.modified_sequence.to_string()) .collect(); + // Cross-impl bisection: dump sorted detected_peptides for diff against + // the C# port. Gated by env var; zero overhead when unset. + if std::env::var("OSPREY_DUMP_DETECTED_PEPTIDES").as_deref() == Ok("1") { + let mut sorted: Vec = detected_peptides.iter().cloned().collect(); + sorted.sort(); + let _ = std::fs::write( + "rust_stage7_detected_peptides.txt", + sorted.join("\n") + "\n", + ); + log::info!( + "[DIAG] Wrote rust_stage7_detected_peptides.txt ({} entries)", + sorted.len() + ); + } + let parsimony = protein::build_protein_parsimony( &library, config.shared_peptides, From cf32a3d224142d9ee7005c91b87bc79c2e432b7f Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Thu, 21 May 2026 08:22:41 -0700 Subject: [PATCH 14/16] Sort per_file_entries by entry_id at run_percolator_fdr entry MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The 1st-pass input is already entry_id-sorted via deduplicate_pairs (pipeline.rs:6123), but the post-rescore pool that feeds 2nd-pass Percolator can have gap-fill entries appended after the sorted pre-existing rows. Re-sorting at the top of run_percolator_fdr guarantees identical iteration order across Rust and OspreySharp; without it, gap-fill ordering can drift the cross-impl 2nd-pass SVM working-set selection on multi-file datasets even when feature columns are bit-equal. On Stellar Single this is a no-op (no reconciliation → no gap-fills). On multi-file the sort shifts the 2nd-pass dump count slightly (Stellar 3-file: 5372 → 5360 proteins in the dump). The remaining cross-impl drift at Stage 7 for multi-file is in code outside the input order and is still under investigation. cargo fmt / clippy / test all pass. --- crates/osprey/src/pipeline.rs | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index 71d1fa8..e2adb46 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -5450,6 +5450,19 @@ fn run_percolator_fdr( log::debug!("Running native Percolator FDR on coelution entries"); + // Sort each file's entries by entry_id so the SVM working-set selection + // sees a canonical order regardless of upstream operation history. The + // 1st-pass input is already entry_id-sorted via deduplicate_pairs + // (pipeline.rs:6123), but the post-rescore pool that feeds 2nd-pass + // Percolator can have gap-fill entries appended after the sorted + // pre-existing rows. Re-sorting here guarantees identical iteration + // order across Rust and OspreySharp; without it, gap-fill ordering + // diverges and the cross-impl 2nd-pass scores drift on multi-file + // datasets even when feature columns are bit-equal. + for (_, entries) in per_file_entries.iter_mut() { + entries.sort_by_key(|e| e.entry_id); + } + let total_entries: usize = per_file_entries.iter().map(|(_, e)| e.len()).sum(); let n_files = per_file_entries.len(); let perc_config = percolator::PercolatorConfig { From 8ef6aa6505a05457d2848d61be90742d0251a852 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Thu, 21 May 2026 09:42:27 -0700 Subject: [PATCH 15/16] Best-per-precursor dedup in direct-path Percolator The streaming Percolator path runs a best-per-precursor dedup before peptide-group subsampling (pipeline.rs:5512-5557) precisely so the SVM trains on one observation per precursor instead of N observations per N-file experiment. The direct path was originally written without that dedup, which is statistically incorrect: on multi-file inputs sized between max_train_size and max_train_size * 2 (e.g. Stellar 3-file at 393k entries) the SVM trained on N-times-redundant precursor pairs as if they were independent samples. The C# port (OspreySharp.FDR.Percolator direct path) had already inferred the correct dedup-then-subsample shape; mirroring it here brings the two implementations into algorithmic agreement and removes the multi-file artefact from the direct path's training set. Dedup key is `features[0]` (fragment_coelution_sum, the first PIN feature) which matches the streaming path's `coelution_sum` field value-for-value. Single-file is unaffected (each base_id has one observation; dedup is a no-op). cargo fmt / clippy / test all pass. --- crates/osprey-fdr/src/percolator.rs | 87 +++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 16 deletions(-) diff --git a/crates/osprey-fdr/src/percolator.rs b/crates/osprey-fdr/src/percolator.rs index ef76050..efda63d 100644 --- a/crates/osprey-fdr/src/percolator.rs +++ b/crates/osprey-fdr/src/percolator.rs @@ -199,22 +199,77 @@ pub fn run_percolator( // two sides, cascading through every downstream computation. dump_stage5_standardizer(&standardizer, config.feature_names.as_deref()); - // 3. Subsample by peptide groups if needed (before fold splitting, per PMC5059416) - // Keeps target-decoy pairs and charge states together. - // The subsampled set is used for fold assignment + SVM training. - // ALL entries are scored with the trained model. - let train_subset: Option> = if config.max_train_size > 0 && n > config.max_train_size - { - Some(subsample_by_peptide_group( - &labels, - &entry_ids, - &peptides, - config.max_train_size, - config.seed, - )) - } else { - None - }; + // 3a. Best-per-precursor: pick the single best-scoring observation per + // (base_id, is_decoy) tuple across all files. With N files per peptide, + // this avoids the SVM seeing the same precursor's target/decoy pair + // N times, which would inflate apparent target/decoy separation and + // cause the SVM to learn file-specific noise rather than peptide + // discriminating features. The streaming path applies this dedup + // inline (pipeline.rs::run_percolator_fdr); applying it here makes + // the direct path statistically consistent. Without this, on multi- + // file inputs below the streaming threshold (Stellar 3-file at 393K + // entries) the SVM trained on N-times-redundant precursor pairs and + // produced inflated target/decoy separation diverging from the + // deduped C# port. + // + // Dedup key is features[0] (fragment_coelution_sum, the first PIN + // feature) which matches the streaming path's `coelution_sum` + // dedup field exactly — same per-entry value. + let mut best_target: HashMap = HashMap::new(); + let mut best_decoy: HashMap = HashMap::new(); + for i in 0..n { + let base_id = entries[i].entry_id & 0x7FFF_FFFF; + let score = entries[i].features[0]; + let map = if labels[i] { + &mut best_decoy + } else { + &mut best_target + }; + match map.get(&base_id) { + Some(&existing) if entries[existing].features[0] >= score => {} + _ => { + map.insert(base_id, i); + } + } + } + let mut dedup_indices: Vec = best_target + .values() + .chain(best_decoy.values()) + .copied() + .collect(); + dedup_indices.sort(); + log::debug!( + " Best-per-precursor: {} entries ({} targets, {} decoys) from {} total", + dedup_indices.len(), + best_target.len(), + best_decoy.len(), + n + ); + + // 3b. Subsample by peptide groups if dedup'd count still exceeds the + // training cap (before fold splitting, per PMC5059416). Keeps target- + // decoy pairs and charge states together. The subsampled set is used + // for fold assignment + SVM training; ALL entries are scored with + // the trained model regardless. + let train_subset: Option> = + if config.max_train_size > 0 && dedup_indices.len() > config.max_train_size { + // Build dedup-local arrays for the subsample call. + let dedup_labels: Vec = dedup_indices.iter().map(|&i| labels[i]).collect(); + let dedup_entry_ids: Vec = dedup_indices.iter().map(|&i| entry_ids[i]).collect(); + let dedup_peptides: Vec = + dedup_indices.iter().map(|&i| peptides[i].clone()).collect(); + let local = subsample_by_peptide_group( + &dedup_labels, + &dedup_entry_ids, + &dedup_peptides, + config.max_train_size, + config.seed, + ); + // Remap local indices back into the original entry index space. + Some(local.into_iter().map(|li| dedup_indices[li]).collect()) + } else { + Some(dedup_indices) + }; let sub_n = train_subset.as_ref().map_or(n, |s| s.len()); From 8712ffea12c39771a9ef27be76a5e55bbfeb82d1 Mon Sep 17 00:00:00 2001 From: brendanx67 Date: Thu, 21 May 2026 15:30:22 -0700 Subject: [PATCH 16/16] Stage 6/7 cross-impl: remap parquet_index after canonical sort + propagate experiment q-values across base_id Two coupled fixes that take 3-file Stellar from Rust 5,366 vs C# 6,541 protein groups (stage 7 FAIL) to bit-equal `.2nd-pass.fdr_scores.bin` sidecars and stage 7 PASS. 1. **`run_percolator_fdr` sort key**: extend the per-file canonical sort from `(entry_id, charge, scan_number)` to `(entry_id, charge, scan_number, parquet_index)`. On 3-file Stellar the post-reconciliation pool contains 94 per-file groups of 2 entries that share all three of the original keys -- a gap-fill rescore landing on the same scan as an original row with a different `rt_deviation`. With only three keys, Rust's stable `sort_by` left those ties in input order, but .NET `List.Sort` is unstable and swapped them, drifting the 2nd-pass standardizer mean by 1 ULP on `rt_deviation` and cascading through every downstream SVM weight. `parquet_index` is intrinsic to the byte-equal cross-impl parquet layout, so adding it as a final tie-break makes the total order identical on both sides regardless of underlying sort stability. 2. **`rescore_per_file_loop` gap-fill remap**: after rewriting the reconciled parquet, every `FdrEntry.parquet_index` whose row was moved by the writer's canonical sort must be remapped to its post- sort position. Previously upstream rows held a stale pre-sort index; when the next Percolator pass loaded features by `parquet_index`, the lookup silently fetched a different entry's feature row. The new code computes the canonical permutation explicitly (matching the writer's sort key), inverts it into pre->post row indices, then walks `fdr_entries` and remaps every stub -- both upstream rows (`old_pq_idx -> pre_to_post[old_pq_idx]`) and gap-fill stubs (via the new `gap_vec_idx_for_pre_sort_row` mapping built during append). Comment on `write_scores_parquet_with_metadata` documents that its internal sort is now stable so external pre-sort + tie-breakers survive the writer. 3. **`compute_experiment_precursor_qvalues` propagation**: the direct Percolator path assigned the winning q-value only to the single `compete_all` winner per `base_id`, leaving every non-winning per-file observation at `q=1.0`. The streaming path (`pipeline.rs::run_percolator_fdr_streaming`) already propagated via its `base_id_exp_prec_q` map; the OspreySharp port matched the streaming semantics. The asymmetry silently broke downstream stages that gate on `experiment_precursor_qvalue` (Stage 6 consensus selection / calibration refit, Stage 7 protein FDR) on multi-file inputs sized below the streaming threshold (Stellar 3-file at 393K entries). Direct path now builds a `base_id -> q` map from winners and propagates to all observations. Verified with `Test-Regression.ps1 -Dataset Stellar -Files All -StartStage stage6 -StopAfterStage stage7 -Tag perside_3file_v4`: stage6 PASS on all four compare dumps (multicharge / consensus / reconciliation / rescored), stage7 PASS on protein FDR. `diff_fdr_bin.py` shows 0 score / q-value diffs cross-impl on every per-file `.2nd-pass.fdr_scores.bin`. A separate one-shot diagnostic dump `dump_stage5_perc_input` (gated by `OSPREY_DUMP_PERC_INPUT=1`) is added to localize future standardizer divergence; writes `rust_stage5_perc_input.tsv` with per-entry raw feature vectors sorted by `(entry_id, native_position)`. --- crates/osprey-fdr/src/percolator.rs | 98 +++++++++++++++++++++++- crates/osprey/src/pipeline.rs | 112 ++++++++++++++++++++++++---- 2 files changed, 192 insertions(+), 18 deletions(-) diff --git a/crates/osprey-fdr/src/percolator.rs b/crates/osprey-fdr/src/percolator.rs index efda63d..dac427e 100644 --- a/crates/osprey-fdr/src/percolator.rs +++ b/crates/osprey-fdr/src/percolator.rs @@ -199,6 +199,12 @@ pub fn run_percolator( // two sides, cascading through every downstream computation. dump_stage5_standardizer(&standardizer, config.feature_names.as_deref()); + // One-shot diagnostic for 2nd-pass divergence localization. Dumps + // per-entry raw feature vectors so the cross-impl compare can see + // exactly which rows differ. Sorted by (entry_id, native_position) + // to align with the C# dump. + dump_stage5_perc_input(entries, config.feature_names.as_deref()); + // 3a. Best-per-precursor: pick the single best-scoring observation per // (base_id, is_decoy) tuple across all files. With N files per peptide, // this avoids the SVM seeing the same precursor's target/decoy pair @@ -1360,7 +1366,21 @@ fn compute_per_run_peptide_qvalues( qvalues } -/// Compute experiment-level precursor q-values (across all files) +/// Compute experiment-level precursor q-values (across all files). +/// +/// Propagates the winner's q-value to every observation sharing the +/// same `base_id` (target and decoy). The streaming path +/// (`run_percolator_fdr_streaming` in `crates/osprey/src/pipeline.rs`) +/// already does this via its `base_id_exp_prec_q` map; the direct +/// path used to assign the q-value only to the single winner, +/// leaving every non-winning per-file observation at q=1.0. That +/// undercount silently broke downstream stages that gate on +/// `experiment_precursor_qvalue` (Stage 6 consensus selection / +/// calibration refit, Stage 7 protein FDR) on multi-file inputs +/// using the direct path (Stellar 3-file, ~393K entries below the +/// 600K streaming threshold). The OspreySharp port matched the +/// streaming path's "propagate to all base_id observations" +/// semantics, so the direct path was the divergent side. fn compute_experiment_precursor_qvalues( scores: &[f64], labels: &[bool], @@ -1375,8 +1395,20 @@ fn compute_experiment_precursor_qvalues( let ws: Vec = winner_indices.iter().map(|&i| scores[i]).collect(); compute_conservative_qvalues(&ws, &winner_is_decoy, &mut q); + // Build base_id -> winner q-value map (one entry per base_id, since + // compete_all picks a single winner per base_id). + let mut base_id_q: HashMap = HashMap::with_capacity(winner_indices.len()); for (rank, &idx) in winner_indices.iter().enumerate() { - qvalues[idx] = q[rank]; + let base_id = entry_ids[idx] & 0x7FFF_FFFF; + base_id_q.insert(base_id, q[rank]); + } + + // Propagate winner's q-value to every observation of that base_id. + for i in 0..n { + let base_id = entry_ids[i] & 0x7FFF_FFFF; + if let Some(&qv) = base_id_q.get(&base_id) { + qvalues[i] = qv; + } } qvalues @@ -1727,6 +1759,68 @@ fn dump_stage5_standardizer(standardizer: &FeatureStandardizer, feature_names: O exit_if_only("OSPREY_STANDARDIZER_ONLY", "Stage 5 standardizer dump"); } +/// One-shot diagnostic dump of the per-entry raw feature vectors that +/// feed `FeatureStandardizer::fit_transform`. Gated by +/// `OSPREY_DUMP_PERC_INPUT=1`. Writes `rust_stage5_perc_input.tsv` with +/// columns `native_position, entry_id, is_decoy, ..` +/// sorted by (entry_id, native_position) for stable cross-impl compare. +/// +/// This is a localizer for cross-impl standardizer divergence: when +/// `rust_stage5_standardizer.tsv` differs but `2nd-pass entries[]` +/// positions match, the divergence is in feature values themselves. +fn dump_stage5_perc_input(entries: &[PercolatorEntry], feature_names: Option<&[String]>) { + if !is_dump_enabled("OSPREY_DUMP_PERC_INPUT") { + return; + } + let path = "rust_stage5_perc_input.tsv"; + let Ok(mut f) = + std::fs::File::create(path).map(|file| std::io::BufWriter::with_capacity(8 << 20, file)) + else { + log::warn!("Could not create {}", path); + return; + }; + write!(f, "native_position\tentry_id\tis_decoy").ok(); + let n_features = if entries.is_empty() { + 0 + } else { + entries[0].features.len() + }; + for i in 0..n_features { + let name = feature_names + .and_then(|n| n.get(i)) + .map(|s| s.as_str()) + .unwrap_or("unknown"); + write!(f, "\t{}", name).ok(); + } + writeln!(f).ok(); + + let mut order: Vec = (0..entries.len()).collect(); + order.sort_by_key(|&i| (entries[i].entry_id, i)); + for i in order { + let e = &entries[i]; + write!( + f, + "{}\t{}\t{}", + i, + e.entry_id, + if e.is_decoy { "true" } else { "false" } + ) + .ok(); + for v in &e.features { + write!(f, "\t{}", format_f64_roundtrip(*v)).ok(); + } + writeln!(f).ok(); + } + let _ = f.flush(); + drop(f); + log::info!( + "Wrote Stage 5 Percolator input dump: {} ({} rows)", + path, + entries.len() + ); + exit_if_only("OSPREY_PERC_INPUT_ONLY", "Stage 5 Percolator input dump"); +} + /// Subsample entries by peptide group, keeping target-decoy pairs and charge states together. /// /// Groups entries by target peptide (via base_id), then randomly selects groups until diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index e2adb46..f08c800 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -3312,23 +3312,84 @@ pub(crate) fn rescore_per_file_loop( } } // Append gap-fill entries (remaining overlay entries beyond parquet range). - // Update the corresponding fdr_entries stubs' parquet_index to point to - // the actual Parquet row they will occupy after write-back. Without this, - // gap-fill stubs keep parquet_index = u32::MAX and Phase 2 of the next - // Percolator run can't load their features. + // The vec_idx in the overlay matches the position in fdr_entries + // (gap-fill entries were appended to fdr_entries in order). + // Track the (vec_idx, pre-sort row) pair so we can remap each + // FdrEntry stub's parquet_index to the post-canonical-sort row. let gap_start_pq_row = full_entries.len(); let mut gap_entries: Vec<(usize, CoelutionScoredEntry)> = overlay.into_iter().collect(); gap_entries.sort_by_key(|(idx, _)| *idx); + let mut gap_vec_idx_for_pre_sort_row: Vec = + Vec::with_capacity(gap_entries.len()); for (gap_offset, (vec_idx, entry)) in gap_entries.into_iter().enumerate() { - let new_pq_row = (gap_start_pq_row + gap_offset) as u32; + let pre_sort_row = gap_start_pq_row + gap_offset; full_entries.push(entry); - // The vec_idx in the overlay matches the position in fdr_entries - // (gap-fill entries were appended to fdr_entries in order) - if vec_idx < fdr_entries.len() { - fdr_entries[vec_idx].parquet_index = new_pq_row; + // Track gap-fill mapping by pre-sort row position. + // Length pads with `usize::MAX` for the upstream rows + // (which already have a real parquet_index pointing to + // the pre-sort row). + while gap_vec_idx_for_pre_sort_row.len() + gap_start_pq_row < pre_sort_row { + gap_vec_idx_for_pre_sort_row.push(usize::MAX); } + gap_vec_idx_for_pre_sort_row.push(vec_idx); } + + // Canonical (entry_id, charge, scan_number) sort permutation. + // Must match write_scores_parquet_with_metadata's internal + // sort key exactly; the writer then iterates in this same + // order, producing a parquet whose physical row layout is + // identical regardless of caller order. + // + // CRITICAL: every FdrEntry stub whose parquet_index points + // into this rewritten parquet must be remapped to its post + // -sort row. Upstream rows held a parquet_index equal to + // their pre-sort row position; the sort moves them, so the + // stale pre-sort index would silently fetch a different + // entry's features in the next Percolator pass. + let mut perm: Vec = (0..full_entries.len()).collect(); + perm.sort_by(|&a, &b| { + full_entries[a] + .entry_id + .cmp(&full_entries[b].entry_id) + .then_with(|| full_entries[a].charge.cmp(&full_entries[b].charge)) + .then_with(|| { + full_entries[a] + .scan_number + .cmp(&full_entries[b].scan_number) + }) + }); + // Inverse permutation: pre_sort_row -> post_sort_row. + let mut pre_to_post: Vec = vec![0u32; full_entries.len()]; + for (post, &pre) in perm.iter().enumerate() { + pre_to_post[pre] = post as u32; + } + // Remap every FdrEntry stub for this file. Upstream stubs + // already had parquet_index pointing to their pre-sort row. + // Gap-fill stubs still have parquet_index = u32::MAX; for + // them we look up the pre-sort row via the vec_idx mapping + // built above (vec_idx -> pre_sort_row inverse). + let mut vec_idx_to_gap_pre_sort: HashMap = HashMap::new(); + for (offset, &vec_idx) in gap_vec_idx_for_pre_sort_row.iter().enumerate() { + if vec_idx != usize::MAX { + vec_idx_to_gap_pre_sort.insert(vec_idx, gap_start_pq_row + offset); + } + } + for (vec_idx, fdr_entry) in fdr_entries.iter_mut().enumerate() { + let pre_sort_row = if fdr_entry.parquet_index == u32::MAX { + // Gap-fill stub: look up pre-sort row by vec_idx. + match vec_idx_to_gap_pre_sort.get(&vec_idx) { + Some(&row) => row, + None => continue, // gap-fill stub without an overlay entry? + } + } else { + fdr_entry.parquet_index as usize + }; + if pre_sort_row < pre_to_post.len() { + fdr_entry.parquet_index = pre_to_post[pre_sort_row]; + } + } + // Write back to Parquet with reconciliation metadata let recon_metadata = build_reconciled_metadata(config, Some(join_file_stems)); let codec = parquet_compression_codec(config.parquet_compression); @@ -5450,17 +5511,36 @@ fn run_percolator_fdr( log::debug!("Running native Percolator FDR on coelution entries"); - // Sort each file's entries by entry_id so the SVM working-set selection - // sees a canonical order regardless of upstream operation history. The + // Sort each file's entries by the parquet canonical key (entry_id, + // charge, scan_number, parquet_index) so the SVM working-set + // selection sees the same order regardless of upstream operation + // history *and* matches the parquet's physical row layout. The // 1st-pass input is already entry_id-sorted via deduplicate_pairs // (pipeline.rs:6123), but the post-rescore pool that feeds 2nd-pass // Percolator can have gap-fill entries appended after the sorted - // pre-existing rows. Re-sorting here guarantees identical iteration - // order across Rust and OspreySharp; without it, gap-fill ordering - // diverges and the cross-impl 2nd-pass scores drift on multi-file - // datasets even when feature columns are bit-equal. + // pre-existing rows. The parquet writer uses (entry_id, charge, + // scan_number) as its sort key, but on Stellar 3-file the post- + // reconciliation pool has 94 per-file groups of 2 entries sharing + // ALL THREE of (entry_id, charge, scan_number) (gap-fill rescore + // landing on the same scan as an original row with a *different* + // rt_deviation). With only those three keys, ties leave a stable + // sort taking input order on each side -- but C# `List.Sort` is + // unstable, so the same input produced a swapped order at every + // tied group, drifting the 2nd-pass standardizer mean by 1 ULP on + // rt_deviation and cascading through all downstream SVM weights. + // Including `parquet_index` as a final tie-break gives the same + // total order on both sides regardless of underlying sort + // stability, since the parquets are byte-identical cross-impl and + // each entry's parquet_index points to its row in that shared + // canonical layout. for (_, entries) in per_file_entries.iter_mut() { - entries.sort_by_key(|e| e.entry_id); + entries.sort_by(|a, b| { + a.entry_id + .cmp(&b.entry_id) + .then_with(|| a.charge.cmp(&b.charge)) + .then_with(|| a.scan_number.cmp(&b.scan_number)) + .then_with(|| a.parquet_index.cmp(&b.parquet_index)) + }); } let total_entries: usize = per_file_entries.iter().map(|(_, e)| e.len()).sum();