diff --git a/pwiz_tools/OspreySharp/Osprey-workflow.html b/pwiz_tools/OspreySharp/Osprey-workflow.html index 8529c2eac0a..9e344ad0bb0 100644 --- a/pwiz_tools/OspreySharp/Osprey-workflow.html +++ b/pwiz_tools/OspreySharp/Osprey-workflow.html @@ -502,10 +502,15 @@

Osprey / OspreySharp DIA pipeline workflow

byte-identical RefSpectraPeaks blobs).

- End-to-end performance (3-file regression, post-sort-fix - binaries): - OspreySharp ties or beats Osprey on most stages of both datasets. - The table below is the headline; per-stage commentary follows it. + End-to-end performance (3-file regression, Windows + native, C: SSD source data, AV-scan excluded; all 8 cells + median-of-3 refreshed 2026-05-24): + Measured on the test/integration-bucket3 branch + (PRs #38-#45 stacked, including the decoy gap-fill exclusion fix + AND the Welford → sum/n revert that recovered ~21 s on Stellar + Rust stage1to4) immediately after 3-file Astral straight-through + cross-impl bit-equality landed at 1e-9. Variance is <5 s on + Stellar and <90 s on Astral across the three repeats.

@@ -520,12 +525,12 @@

Osprey / OspreySharp DIA pipeline workflow

- - - - - - + + + + + +
stage1to41:301:080.76× (C# faster)8:2611:171.34×
stage5 (1st-pass FDR + plan)1:441:531.09× (~tied)3:302:450.79× (C# faster)
stage6 (rescore + gap-fill)4:103:340.86× (C# faster)12:4517:301.37×
stage7 (2nd-pass FDR + protein)0:030:031.0× (tied)0:160:130.81× (C# faster)
blib write0:110:080.71× (C# faster)0:550:450.81× (C# faster)
Total7:386:460.89× (C# faster)~25:51~32:351.26×
stage1to41:421:130.72x (C# faster)7:546:240.81x (C# faster)
stage5 (1st-pass FDR + plan)1:142:041.68x1:275:123.59x
stage6 (rescore + gap-fill)0:340:270.81x (C# faster)7:103:490.53x (C# faster)
stage7 (2nd-pass FDR + protein)0:430:501.17x2:543:071.07x
blib write0:190:020.14x (C# faster)2:020:070.06x (C# faster)
Total4:334:381.02x (~tied)21:2918:410.87x (C# faster)

@@ -533,44 +538,236 @@

Osprey / OspreySharp DIA pipeline workflow

+

+ WSL Linux end-to-end performance — ext4 inside + /home (median-of-3 2026-05-24): + Same workstation as the Windows benchmark above, WSL2 Ubuntu + 22.04, .NET 8.0.421, Rust 1.95.0. The WSL ext4.vhdx + is on the C: SSD (relocated this session via + wsl --export/--import from a D: HDD + location, where it had drifted at some point between the prior + 2026-05-18 measurement and now). Test data is staged inside + the VHDX at /home/brendanx/test/osprey-runs, so + every read and write goes through Linux ext4 with no 9P drvfs + crossing. This is the “Linux-native” configuration: + Rust is 23-25% faster than C# in this regime (Stellar 3:56 vs + 4:52; Astral 15:38 vs 19:31), matching the historical pattern + of native Linux ext4 favouring the Rust toolchain. Cross-OS + output parity is gated separately + by Test-Snapshot.ps1 at 1e-6 + tolerance on four median-polish columns + (median_polish_cosine, _residual_ratio, + _min_fragment_r2, _residual_correlation) + because both runtimes delegate Math.Log and + Math.Exp to the host libm (Linux glibc vs Windows ucrt), + which differ at f64 ULP level. All other Stage 1-4 columns are + bit-equal cross-OS, and the cross-OS divergence in stage5 percolator + TSV output is bounded to a handful of formatted-float characters + per ~200 MB file. +

+ + + + + + + + + + + + + + + + + + + + +
StageStellar RustStellar C#C#/RustAstral RustAstral C#C#/Rust
stage1to41:441:190.76x (C# faster)6:406:471.02x (~tied)
stage5 (1st-pass FDR + plan)1:052:091.98x1:145:174.25x
stage6 (rescore + gap-fill)0:210:271.27x3:543:530.99x (~tied)
stage7 (2nd-pass FDR + protein)0:390:521.33x2:553:241.17x
blib write0:060:030.54x (C# faster)0:530:090.18x (C# faster)
Total3:564:521.23x (Rust faster)15:3819:311.25x (Rust faster)
+

+ WSL Linux end-to-end performance — 9P drvfs to + Windows /mnt/c (median-of-3 2026-05-24): + Same WSL distro and binaries, but with test data living on the + Windows side at C:\test\osprey-runs (accessed as + /mnt/c/test/osprey-runs from inside WSL). Every + read and write crosses the 9P protocol bridge between Linux and + Windows NTFS — this is the configuration developers get by + default when they keep their data on the Windows filesystem and + just run the pipeline from a WSL shell. 9P adds substantial + per-syscall overhead (240 vs 14 files/s for fsync'd metadata + writes; raw I/O bench below) but the absolute totals are still + workable, and the C#/Rust ratio collapses to near 1.0× on + both datasets — C# and Rust are essentially tied here, and + on Stellar C# is marginally faster. The 9P penalty applies to + both implementations equally, so the cross-impl story is + unchanged: C# is not a performance sacrifice on the storage path + most users have. +

+ + + + + + + + + + + + + + + + + + + + +
StageStellar RustStellar C#C#/RustAstral RustAstral C#C#/Rust
stage1to43:112:030.65x (C# faster)12:519:490.76x (C# faster)
stage5 (1st-pass FDR + plan)1:062:122.00x1:175:234.20x
stage6 (rescore + gap-fill)0:530:390.75x (C# faster)5:415:020.88x (C# faster)
stage7 (2nd-pass FDR + protein)0:480:511.07x (~tied)3:083:301.11x
blib write0:060:091.40x0:560:240.43x (C# faster)
Total6:055:560.97x (~tied)23:5424:091.01x (~tied)
+

+ A raw 9P / ext4 I/O bench from this session + (ai/.tmp/wsl-io-bench/summary.md) explains the + spread — ext4 inside the VHDX is dramatically faster than + drvfs for both sequential writes and the small-fsync metadata + pattern that stage5/6 dump writers used to hit: +

+ + + + + + + + + + + + + + + + +
WSL targetBackingSeq writeSeq read (O_DIRECT)500×4KiB fsync
/mnt/cC: SSD via 9P drvfs356 MB/s431 MB/s190 files/s
/mnt/dD: HDD via 9P drvfs215 MB/s140 MB/s30 files/s
/homeext4 / VHDX / C: SSD2064 MB/s ‡12328 MB/s †cached
/dev/shmtmpfs (RAM)3333 MB/s10085 MB/scached
+

+ † Cache-influenced read: the 4 GiB test file is fully cached + after the immediately preceding write on a 64 GB workstation; + O_DIRECT on 9p drvfs does not always bypass the + Windows-side NTFS cache. Sustained read on the underlying HDD + platters is ~200 MB/s. + ‡ fdatasync forces Linux dirty pages out to + /dev/sdd (the VHDX virtual disk), not all the way + through the Windows VHDX layer to the underlying SSD; large + sustained writes that exceed the page cache would converge to + SSD write speed. +

+

+ Reference for past closure (2026-05-17 to 2026-05-21 work): + Two parallelism / allocation fixes closed an earlier WSL gap: + (a) parallelizing the c-value loop in + PercolatorFdr.GridSearchC (was serial; Rust used + par_iter) — ~1:30 savings on Stellar and ~0:15 on + Astral stage5; + (b) wiring SpectralScorer.PreprocessSpectrumForXcorrInto + to use its passed-in XcorrScratch (it was nominally + scratch-aware but actually allocated four float[NBins] + arrays per call then copied) — ~5:14 savings on Astral stage6 and + ~1:20 on Astral stage1to4. Both changes are behavior-preserving + (cross-impl Test-Snapshot remained bit-exact). Stellar runs through + UnitStrategy which is not affected by the f32 scratch + path, so the Astral-targeted fix produces no regression there. +

+

+ Methodology notes. + (1) Storage layout. All numbers in the three tables + above are on C: SSD physical media; both Windows + C:\test\osprey-runs and the WSL ext4.vhdx + (now backing both /home and /mnt/c + via different paths) live on the same SK hynix 2 TB NVMe. + Defender real-time scan exclusions cover both + C:\test and the Windows mount points the WSL distro + touches; without those exclusions, every disk-bound stage paid + a 5-10% scan overhead. + (2) Rust dump-writer patch (still load-bearing). Rust's + stage-5/6 diagnostic-dump writers in osprey v26.6.0 use unbuffered + std::fs::File::create: each writeln! + crosses to the kernel. On WSL drvfs/9P that took ~14 min for + the 74 MB stage5 percolator dump alone (vs ~1 min for the + percolator algorithm itself); BufWriter eliminates this on either + OS. The numbers above reflect that patch (committed via PR #36) + plus an explicit f.flush(); drop(f); before + exit_if_only. + (3) Welford → sum/n MS2 calibration revert. An earlier + attempt to close the Astral 3-file cross-impl bit-equality gate + used Welford's online running mean for MS2 calibration; this + avoided LLVM/JIT vectorisation divergence at the cost of a + 4-op-per-element loop that defeated SIMD on both sides + (~21 s on Stellar Rust stage1to4, more on Astral). The same + bit-equality goal is now met by a plain sum / n with + a deterministic (base_id, entry_id) sort upstream and + IEEE 754 strict mode (Rust default, .NET default) preventing + re-association. Stellar + Astral 3-file end-to-end pass at 1e-9 + cross-impl with the simpler implementation, and the perf table + above is on the post-revert binaries. + (4) C# Astral parallelism. C# Astral was run with + OSPREY_MAX_PARALLEL_FILES=1 — the Windows reference + above runs three files concurrently in stage 1-4, which on this + 64 GB workstation requires ~60 GB peak under WSL (which fits + ~30 GB of overhead in the Hyper-V + ext4.vhdx machinery). Without + the cap, Stage 1-4 is OOM-killed in WSL. Rust naturally processes + files sequentially and is unaffected. A fully fair + WSL-vs-Windows comparison would set + OSPREY_MAX_PARALLEL_FILES=1 on both sides too. +

BLIB writer: OspreySharp's BlibWriter uses DotNetZip / Ionic.Zlib level 6 (the same library Skyline's diff --git a/pwiz_tools/OspreySharp/OspreySharp.Chromatography/LoessRegression.cs b/pwiz_tools/OspreySharp/OspreySharp.Chromatography/LoessRegression.cs index e55610a9823..9f4b5c02c14 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Chromatography/LoessRegression.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Chromatography/LoessRegression.cs @@ -246,13 +246,19 @@ public static LoessModel Fit(double[] x, double[] y, double bandwidth = 0.3, if (x.Length < 2) throw new ArgumentException("Need at least 2 data points"); - // Sort by x (stable, matching Rust's slice::sort_by). Array.Sort - // with Comparison is unstable (introsort) and reorders ties - // differently than Rust for duplicate x values, causing LOESS - // divergence on data with repeated x (e.g. multi-charge peptides - // sharing a library RT). + // Sort by x with a secondary key on y, so the order is + // deterministic for duplicate x values (e.g. multi-charge + // peptides sharing a library RT). LINQ OrderBy / ThenBy are + // stable, but the upstream input order can vary across impls + // when discriminant-score-sorted ties get filtered at 1 ULP, + // so we cannot rely on input-order tiebreaking for cross-impl + // bit-equality. Matches Rust osprey-chromatography + // calibration/rt.rs which sorts by (x, y) tuple. int n = x.Length; - int[] order = Enumerable.Range(0, n).OrderBy(i => x[i]).ToArray(); + int[] order = Enumerable.Range(0, n) + .OrderBy(i => x[i]) + .ThenBy(i => y[i]) + .ToArray(); double[] sortedX = new double[n]; double[] sortedY = new double[n]; diff --git a/pwiz_tools/OspreySharp/OspreySharp.Chromatography/MzCalibration.cs b/pwiz_tools/OspreySharp/OspreySharp.Chromatography/MzCalibration.cs index 4dbb93743a9..1676d1e9b5b 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Chromatography/MzCalibration.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Chromatography/MzCalibration.cs @@ -257,27 +257,51 @@ private static MzCalibrationResult CalculateSingleCalibration(double[] errors, int n = errors.Length; - // Calculate mean - double sum = 0; + // Naive sum/n with consistent left-to-right accumulation. + // Mirror of osprey-chromatography/src/calibration/mass.rs. + // + // For ppm mass-error inputs the running sum stays well below + // the f64 precision wall (worst case ~18 k errors × ~0.22 ppm + // ≈ 4 000, f64 ULP ~5e-13, ~9 digits headroom for the mean — + // well clear of the 1e-9 cross-impl gate). + // + // The caller sorts the error list by (base_id, entry_id) + // before this method runs, matching the Rust caller's + // deterministic ordering. With deterministic order and .NET's + // default IEEE 754 strict mode (no `-ffast-math`, no + // associative re-ordering), the JIT cannot vectorise the + // dependent reduction `sum += errors[i]`, so the loop stays + // serial-scalar and bit-equal to the Rust mirror which is + // similarly constrained by Rust's default IEEE 754 mode. + // + // The earlier Welford-Knuth recurrence was bit-equal to Rust + // but cost ~21 s extra on Stellar Rust stage1to4 (4 ops per + // element vs 1) for stability we did not need at ppm scale. + double sum = 0.0; for (int i = 0; i < n; i++) + { sum += errors[i]; + } double mean = sum / n; - // Calculate median + // Sample variance via a second left-to-right walk over + // (x - mean)^2 contributions. Two passes total still beat + // Welford's single-pass 4-op recurrence on per-element cost. + double sumSqDev = 0.0; + for (int i = 0; i < n; i++) + { + double d = errors[i] - mean; + sumSqDev += d * d; + } + double m2 = sumSqDev; + + // Calculate median (sort-based, unaffected by the Welford change) double median = LoessRegression.Median(errors); - // Calculate standard deviation (sample SD) + // Sample standard deviation: sqrt(M2 / (n-1)). double variance = 0; if (n > 1) - { - double sumSq = 0; - for (int i = 0; i < n; i++) - { - double d = errors[i] - mean; - sumSq += d * d; - } - variance = sumSq / (n - 1); - } + variance = m2 / (n - 1); double sd = Math.Sqrt(variance); return new MzCalibrationResult diff --git a/pwiz_tools/OspreySharp/OspreySharp.Chromatography/RTCalibration.cs b/pwiz_tools/OspreySharp/OspreySharp.Chromatography/RTCalibration.cs index c52e272ba2d..a15fcc29c1d 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Chromatography/RTCalibration.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Chromatography/RTCalibration.cs @@ -117,7 +117,16 @@ public RTCalibration Fit(double[] libraryRts, double[] measuredRts) // for each duplicate and the subsequent LOESS fit diverges from // Rust. LINQ OrderBy is stable and matches Rust. int n = libraryRts.Length; - int[] order = Enumerable.Range(0, n).OrderBy(i => libraryRts[i]).ToArray(); + // Sort by (libraryRt, measuredRt) so the outer x/y arrays match + // the inner sort inside LoessRegression.Fit (which also uses + // (x, y) for duplicate-x determinism). Without the secondary + // key, outer y[i] and inner fitted[i] can end up corresponding + // to different data points at duplicate-x positions, producing + // swapped abs_residuals[i] in the calibration model. + int[] order = Enumerable.Range(0, n) + .OrderBy(i => libraryRts[i]) + .ThenBy(i => measuredRts[i]) + .ToArray(); double[] x = new double[n]; double[] y = new double[n]; diff --git a/pwiz_tools/OspreySharp/OspreySharp.Core/OspreyConfig.cs b/pwiz_tools/OspreySharp/OspreySharp.Core/OspreyConfig.cs index 86d16cd70ca..cd6a1fa3d90 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Core/OspreyConfig.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Core/OspreyConfig.cs @@ -387,6 +387,33 @@ internal static string EscapeForRustDebug(string s) /// hash invariant to invocation order). /// public string ReconciliationParameterHash() + { + var stems = new List(InputFiles?.Count ?? 0); + if (InputFiles != null) + { + foreach (var path in InputFiles) + { + string stem = Path.GetFileNameWithoutExtension(path); + if (!string.IsNullOrEmpty(stem)) + stems.Add(stem); + } + } + return ReconciliationParameterHashForStems(stems); + } + + ///

+ /// Compute the reconciliation parameter hash for an explicit set of + /// file stems. Used by per-file Stage 6 rescore workers, whose + /// only carries this worker's single + /// parquet — the hash that the downstream --join-at-pass=2 + /// merge node expects is computed over ALL files in the join, so + /// the worker must read the full set from the planner's + /// reconciliation.json envelope and pass it in here. The + /// stems are sorted + deduped internally so the hash is invariant + /// to caller ordering. Mirrors Rust + /// OspreyConfig::reconciliation_parameter_hash_for_stems. + /// + public string ReconciliationParameterHashForStems(IReadOnlyList fileStems) { using (var sha256 = SHA256.Create()) { @@ -401,17 +428,31 @@ public string ReconciliationParameterHash() // Mirror Rust's `format!("file_stems:{:?}\n", stems)` output // exactly. {:?} on Vec yields ["a", "b"] with the // brackets and double-quoted, comma-space-separated values. - var stems = new List(InputFiles?.Count ?? 0); - if (InputFiles != null) + // Stems are sorted + deduped here so the hash matches the + // Rust side, which also sorts + dedups before hashing. + var stems = new List(fileStems?.Count ?? 0); + if (fileStems != null) { - foreach (var path in InputFiles) + foreach (var stem in fileStems) { - string stem = Path.GetFileNameWithoutExtension(path); if (!string.IsNullOrEmpty(stem)) stems.Add(stem); } } stems.Sort(StringComparer.Ordinal); + // Dedup in place (stems is sorted, so duplicates are + // adjacent). Rust does `dedup()` on a sorted Vec; same here. + int write = 0; + for (int read = 0; read < stems.Count; read++) + { + if (read == 0 || !string.Equals(stems[read], stems[read - 1], StringComparison.Ordinal)) + { + stems[write++] = stems[read]; + } + } + if (write < stems.Count) + stems.RemoveRange(write, stems.Count - write); + var stemsList = new StringBuilder("["); for (int i = 0; i < stems.Count; i++) { diff --git a/pwiz_tools/OspreySharp/OspreySharp.FDR/FdrDiagnostics.cs b/pwiz_tools/OspreySharp/OspreySharp.FDR/FdrDiagnostics.cs new file mode 100644 index 00000000000..38dd00ed356 --- /dev/null +++ b/pwiz_tools/OspreySharp/OspreySharp.FDR/FdrDiagnostics.cs @@ -0,0 +1,125 @@ +/* + * Original author: Brendan MacLean , + * MacCoss Lab, Department of Genome Sciences, UW + * AI assistance: Claude Code (Claude Opus 4.7) + * + * Based on osprey (https://github.com/MacCossLab/osprey) + * by Michael J. MacCoss, MacCoss Lab, Department of Genome Sciences, UW + * + * Copyright 2026 University of Washington - Seattle, WA + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +// Cross-impl bisection dumps for FDR-level diagnostics. Functions in this +// class are env-var-gated by static bools (evaluated once at class load) +// and no-op in production runs. They write small TSVs to the current +// working directory so they can be diffed against the Rust osprey-fdr +// crate's matching dumps in osprey-fdr/src/diagnostics.rs. +// +// This is a per-project diagnostics class for OspreySharp.FDR; the +// top-level project has its own OspreyDiagnostics (which cannot be +// referenced from here due to layering). Naming kept distinct +// (FdrDiagnostics vs Diagnostics in OspreySharp.Core) to avoid +// collision when both namespaces are imported. + +using System; +using System.Collections.Generic; +using System.Globalization; +using System.IO; +using pwiz.OspreySharp.Core; + +namespace pwiz.OspreySharp.FDR +{ + public static class FdrDiagnostics + { + private static bool IsOne(string name) + { + return Environment.GetEnvironmentVariable(name) == @"1"; + } + + /// + /// OSPREY_DUMP_STAGE7_WINNERS: dump the full cumulative-FDR winners + /// list (target + decoy together) to cs_stage7_winners.tsv after the + /// sort in ComputeProteinFdr. Columns: rank, score, is_decoy, + /// raw_qvalue, monotonic_qvalue. The existing + /// WriteStage7ProteinFdrDump (in OspreyDiagnostics) emits only target + /// winners' scores; decoy-winner scores are not exposed there, + /// hiding cross-impl divergences driven by decoy-winner scores or + /// sort-position interleaving in the cumulative sweep. + /// + public static readonly bool DumpStage7Winners = IsOne(@"OSPREY_DUMP_STAGE7_WINNERS"); + + /// + /// OSPREY_DUMP_BEST_PEPTIDE_SCORES: dump the per-modseq aggregated + /// best-score map from CollectBestPeptideScores to + /// cs_best_peptide_scores.tsv. Surfaces the protein-FDR input set so + /// upstream aggregation divergences (e.g. different per-peptide max + /// scores from compaction asymmetry) can be diffed directly. + /// + public static readonly bool DumpBestPeptideScores = IsOne(@"OSPREY_DUMP_BEST_PEPTIDE_SCORES"); + + /// + /// Write cs_stage7_winners.tsv. Caller passes a (score, is_decoy) + /// tuple list in sort order plus the parallel q-value arrays. Check + /// first to skip the LINQ projection + /// on the disabled-dump path. + /// + public static void WriteStage7WinnersDump( + IList<(double Score, bool IsDecoy)> winners, + double[] rawQvalues, + double[] monotonicQvalues) + { + const string path = @"cs_stage7_winners.tsv"; + var inv = CultureInfo.InvariantCulture; + using (var sw = new StreamWriter(path)) + { + sw.WriteLine("rank\tscore\tis_decoy\traw_qvalue\tmonotonic_qvalue"); + for (int i = 0; i < winners.Count; i++) + { + sw.WriteLine(string.Format(inv, "{0}\t{1}\t{2}\t{3}\t{4}", + i, + Diagnostics.FormatF64Roundtrip(winners[i].Score), + winners[i].IsDecoy ? "true" : "false", + Diagnostics.FormatF64Roundtrip(rawQvalues[i]), + Diagnostics.FormatF64Roundtrip(monotonicQvalues[i]))); + } + } + } + + /// + /// Write cs_best_peptide_scores.tsv. Rows sorted by + /// modified_sequence for stable cross-impl diff. + /// + public static void WriteBestPeptideScoresDump(Dictionary best) + { + const string path = @"cs_best_peptide_scores.tsv"; + var inv = CultureInfo.InvariantCulture; + var keys = new List(best.Keys); + keys.Sort(StringComparer.Ordinal); + using (var sw = new StreamWriter(path)) + { + sw.WriteLine("modified_sequence\tscore\tis_decoy\tbest_qvalue"); + foreach (var seq in keys) + { + var ps = best[seq]; + sw.WriteLine(string.Format(inv, "{0}\t{1}\t{2}\t{3}", + seq, + Diagnostics.FormatF64Roundtrip(ps.Score), + ps.IsDecoy ? "true" : "false", + Diagnostics.FormatF64Roundtrip(ps.BestQvalue))); + } + } + } + } +} diff --git a/pwiz_tools/OspreySharp/OspreySharp.FDR/PercolatorFdr.cs b/pwiz_tools/OspreySharp/OspreySharp.FDR/PercolatorFdr.cs index f9fdf809c17..23b046321be 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.FDR/PercolatorFdr.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.FDR/PercolatorFdr.cs @@ -248,12 +248,34 @@ public static PercolatorResults RunPercolator( } } + // One-shot diagnostic for 2nd-pass divergence localization. + // Gated by OSPREY_DUMP_PERC_INPUT=1; exits via + // OSPREY_PERC_INPUT_ONLY=1. Dumps the raw per-entry feature + // vectors fed into the standardizer so cross-impl compare + // can pinpoint which rows differ. + if (string.Equals(Environment.GetEnvironmentVariable(@"OSPREY_DUMP_PERC_INPUT"), @"1")) + { + WriteStage5PercInputDump(entries, config.FeatureNames); + if (string.Equals(Environment.GetEnvironmentVariable(@"OSPREY_PERC_INPUT_ONLY"), @"1")) + { + Console.Error.WriteLine(@"[BISECT] OSPREY_PERC_INPUT_ONLY set - aborting after dump"); + Environment.Exit(0); + } + } + // 3a. Best-per-precursor: pick the single best-scoring observation per // (base_id, isDecoy) 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. Matches the Rust streaming Percolator path. + // discriminating features. Mirrors the streaming Percolator path's + // dedup step (RunPercolatorStreaming); the Rust direct path + // historically omitted this step, but on multi-file inputs sized + // below the streaming threshold (Stellar 3-file at 393k entries) + // the omission produced a statistically incorrect training set + // that treated multi-file repeats of the same precursor as + // independent samples. Rust was patched to match this dedup + // (osprey-fdr/src/percolator.rs::run_percolator direct path). int[] bestPerPrecursor = SelectBestPerPrecursor(labels, entryIds, entries); int dedupTargets = 0, dedupDecoys = 0; @@ -400,17 +422,32 @@ public static PercolatorResults RunPercolator( foldTrainIndices[fold] = list.ToArray(); } + // One scratch pool for the whole outer-fold Parallel.For and + // every nested GridSearchC.Parallel.For below. Initial size + // = the subsampled training set; grid_search inner SVMs may + // need a different (typically smaller) capacity, handled by + // SvmTrainScratch.EnsureCapacity on rent. Pool grows + // organically to the parallel-worker high-water mark and + // arrays stay in gen-2 LOH for the rest of the run. + var svmScratchPool = new SvmTrainScratchPool(subN, nFeatures); + // Train all folds in parallel. Each fold reads from the shared // subFeatures matrix (read-only) and produces an independent model. - // Matches the Rust implementation's into_par_iter() over folds. + // Mirrors Rust's into_par_iter(). Use OspreyParallel.For (explicit + // dedicated threads) rather than TPL Parallel.For: the TPL + // TaskReplicator was throttling effective parallelism to ~2.5x + // on HRAM Astral (vs Rust rayon's ~9x) even with the same + // per-call cost. Explicit threads remove the ThreadPool + // scheduling variable. var swTrain = Stopwatch.StartNew(); - System.Threading.Tasks.Parallel.For(0, config.NFolds, fold => + OspreyParallel.For(0, config.NFolds, config.NFolds, fold => { var swFold = Stopwatch.StartNew(); int iters; foldModels[fold] = TrainFold( subFeatures, subLabels, subEntryIds, subPeptides, - foldTrainIndices[fold], initialScores, config, trainFdr, out iters); + foldTrainIndices[fold], initialScores, config, trainFdr, + svmScratchPool, out iters); foldIterations[fold] = iters; swFold.Stop(); foldElapsed[fold] = swFold.Elapsed.TotalSeconds; @@ -729,7 +766,7 @@ public static PercolatorResults ScorePopulationAndComputeFdr( var pepOrder = new int[nWinners]; for (int k = 0; k < nWinners; k++) pepOrder[k] = k; - Array.Sort(pepOrder, (a, b) => // Array.Sort OK: TDC's CompeteAll already produced one winner per base_id, so each base_id appears at most once in pepOrder — no ties. + Array.Sort(pepOrder, (a, b) => // Array.Sort OK: TDC's CompeteAll already produced one winner per base_id, so each base_id appears at most once in pepOrder -- no ties. { uint ba = entryIds[winnerIndices[a]] & BASE_ID_MASK; uint bb = entryIds[winnerIndices[b]] & BASE_ID_MASK; @@ -813,11 +850,19 @@ private static LinearSvmClassifier TrainFold( double[] initialScores, PercolatorConfig config, double trainFdr, + SvmTrainScratchPool svmScratchPool, out int bestIteration) { int nFeatures = stdFeatures.Cols; var currentScores = (double[])initialScores.Clone(); + // Rent one scratch for this outer fold's sequential Train calls + // (the final per-iteration Train at the bottom of the loop). The + // inner parallel grid search rents its own scratches from the + // same pool. Return at the end of the fold. + var foldScratch = svmScratchPool != null ? svmScratchPool.Rent() : null; + try { + var bestModel = LinearSvmClassifier.Train( Matrix.Zeros(0, nFeatures), new bool[0], 1.0, config.Seed); bestIteration = 0; @@ -859,7 +904,19 @@ private static LinearSvmClassifier TrainFold( for (int i = 0; i < svmIndices.Count; i++) svmGlobalIndices[i] = trainIndices[svmIndices[i]]; - var svmFeatures = ExtractRows(stdFeatures, svmGlobalIndices); + // svmFeatures is live from here through the Train call below + // (used by Train + DecisionFunction). Use foldScratch.TrainData + // to avoid an 8+ MB LOH allocation per TrainFold iteration. + Matrix svmFeatures; + if (foldScratch != null) + { + foldScratch.EnsureExtractCapacity(svmGlobalIndices.Length, nFeatures); + svmFeatures = ExtractRowsInto(stdFeatures, svmGlobalIndices, foldScratch.TrainData); + } + else + { + svmFeatures = ExtractRows(stdFeatures, svmGlobalIndices); + } var svmLabels = new bool[svmIndices.Count]; var svmEntryIds = new uint[svmIndices.Count]; for (int i = 0; i < svmIndices.Count; i++) @@ -878,21 +935,33 @@ private static LinearSvmClassifier TrainFold( double bestC = GridSearchC( svmFeatures, svmLabels, svmEntryIds, config.CValues, svmFoldAssignments, config.NFolds, - config.Seed, trainFdr); + config.Seed, trainFdr, svmScratchPool); // iii. Train SVM with best C var model = LinearSvmClassifier.Train( - svmFeatures, svmLabels, bestC, config.Seed); + svmFeatures, svmLabels, bestC, config.Seed, foldScratch); // iv. Score ALL training set entries with new model - var trainFeatures = ExtractRows(stdFeatures, trainIndices); + // trainFeatures is live just for the DecisionFunction call; + // foldScratch.TestData is not used elsewhere in this iteration + // (svmFeatures is in TrainData), so reuse it here. + Matrix trainFeatures; + if (foldScratch != null) + { + foldScratch.EnsureExtractCapacity(trainIndices.Length, nFeatures); + trainFeatures = ExtractRowsInto(stdFeatures, trainIndices, foldScratch.TestData); + } + else + { + trainFeatures = ExtractRows(stdFeatures, trainIndices); + } var newTrainScores = model.DecisionFunction(trainFeatures); for (int i = 0; i < trainIndices.Length; i++) currentScores[trainIndices[i]] = newTrainScores[i]; // v. Count passing targets - int nPassing = CountPassing(newTrainScores, trainLabels, trainEntryIds, trainFdr); + int nPassing = CountPassing(newTrainScores, trainLabels, trainEntryIds, trainFdr, foldScratch); if (nPassing > bestPassing) { @@ -912,6 +981,11 @@ private static LinearSvmClassifier TrainFold( bestIteration = Math.Max(bestIteration, 1); return bestModel; + + } finally { + if (foldScratch != null && svmScratchPool != null) + svmScratchPool.Return(foldScratch); + } } // ============================================================ @@ -1090,25 +1164,226 @@ private static void ComputeQvalues( public static int CountPassing( double[] scores, bool[] labels, uint[] entryIds, double fdrThreshold) { - var allIndices = new int[scores.Length]; + return CountPassing(scores, labels, entryIds, fdrThreshold, null); + } + + /// + /// Overload that reuses pre-allocated buffers from a + /// . Pass null + /// to allocate per-call (the legacy path). For the hot Percolator + /// path (CountPassing is called ~570x per grid-search session), + /// passing scratch eliminates ~400 KB of per-call LOH allocation + /// (int[scores.Length] + double[winners]) plus the + /// CompeteFromIndices internal allocations via the scratch-aware + /// helper below. + /// + public static int CountPassing( + double[] scores, bool[] labels, uint[] entryIds, double fdrThreshold, + SvmTrainScratch scratch) + { + if (scratch == null) + { + // Allocating path -- preserved verbatim for callers + // that don't have a scratch (tests, non-hot sites). + var allIndices = new int[scores.Length]; + for (int i = 0; i < scores.Length; i++) + allIndices[i] = i; + + int[] wi; + double[] ws; + bool[] wd; + CompeteFromIndices(scores, labels, entryIds, allIndices, out wi, out ws, out wd); + + var qValues = new double[wi.Length]; + ComputeQvalues(ws, wd, qValues); + + int count = 0; + for (int rank = 0; rank < wi.Length; rank++) + { + if (!labels[wi[rank]] && qValues[rank] <= fdrThreshold) + count++; + } + return count; + } + + scratch.EnsureCountPassingCapacity(scores.Length); + int[] allIdx = scratch.CountPassingIndices; for (int i = 0; i < scores.Length; i++) - allIndices[i] = i; + allIdx[i] = i; - int[] wi; - double[] ws; - bool[] wd; - CompeteFromIndices(scores, labels, entryIds, allIndices, out wi, out ws, out wd); + int winnerCount = CompeteFromIndicesInto( + scores, labels, entryIds, allIdx, scores.Length, scratch); - var qValues = new double[wi.Length]; - ComputeQvalues(ws, wd, qValues); + double[] qVals = scratch.CountPassingQvalues; + // ComputeQvalues operates on a winner-sized slice; pass the + // prefix of the pooled arrays (Compute reads scores[i] for + // i in [0, n), assuming n = winnerCount). + ComputeQvaluesInto( + scratch.CompetitionWinnerScores, scratch.CompetitionWinnerIsDecoy, + qVals, winnerCount); - int count = 0; - for (int rank = 0; rank < wi.Length; rank++) + int[] winIdx = scratch.CompetitionWinnerIndices; + int passCount = 0; + for (int rank = 0; rank < winnerCount; rank++) { - if (!labels[wi[rank]] && qValues[rank] <= fdrThreshold) - count++; + if (!labels[winIdx[rank]] && qVals[rank] <= fdrThreshold) + passCount++; + } + return passCount; + } + + /// + /// Scratch-pooled internal variant of . + /// Writes winners into 's three + /// CompetitionWinner* arrays (prefix [0..returned count) is + /// active). Same algorithm as the allocating version; only the + /// output destination differs. Returns the active winner count. + /// + private static int CompeteFromIndicesInto( + double[] scores, bool[] labels, uint[] entryIds, + int[] indices, int indicesCount, + SvmTrainScratch scratch) + { + // Allocate the small per-call dictionaries / list at full + // expected capacity to avoid rehash growth. Could be pooled + // on scratch in a follow-up; the n*p allocations above are + // the bigger LOH issue. + var targets = new Dictionary>(indicesCount / 2); + var decoys = new Dictionary>(indicesCount / 2); + + for (int ii = 0; ii < indicesCount; ii++) + { + int idx = indices[ii]; + uint baseId = entryIds[idx] & BASE_ID_MASK; + double s = scores[idx]; + if (labels[idx]) + { + KeyValuePair existing; + if (decoys.TryGetValue(baseId, out existing)) + { + if (s > existing.Value) + decoys[baseId] = new KeyValuePair(idx, s); + } + else + { + decoys[baseId] = new KeyValuePair(idx, s); + } + } + else + { + KeyValuePair existing; + if (targets.TryGetValue(baseId, out existing)) + { + if (s > existing.Value) + targets[baseId] = new KeyValuePair(idx, s); + } + else + { + targets[baseId] = new KeyValuePair(idx, s); + } + } + } + + // Walk pairs into local struct array (parallel-array layout + // avoids the per-element Tuple class allocation that the + // public CompeteFromIndices pays). + int maxWinners = targets.Count + decoys.Count; + scratch.EnsureCountPassingCapacity(maxWinners); + int[] winIdx = scratch.CompetitionWinnerIndices; + double[] winScores = scratch.CompetitionWinnerScores; + bool[] winDecoy = scratch.CompetitionWinnerIsDecoy; + // baseIds for tie-break ordering; reuse CountPassingIndices + // as a uint[] surrogate (interpret bits). Cleaner: small + // separate buffer; for now allocate per-call (small). + var winBaseIds = new uint[maxWinners]; + + int n = 0; + foreach (var kvp in targets) + { + uint baseId = kvp.Key; + int tIdx = kvp.Value.Key; + double tScore = kvp.Value.Value; + KeyValuePair de; + if (decoys.TryGetValue(baseId, out de)) + { + if (tScore > de.Value) + { winIdx[n] = tIdx; winScores[n] = tScore; winDecoy[n] = false; winBaseIds[n] = baseId; n++; } + else + { winIdx[n] = de.Key; winScores[n] = de.Value; winDecoy[n] = true; winBaseIds[n] = baseId; n++; } + } + else + { + winIdx[n] = tIdx; winScores[n] = tScore; winDecoy[n] = false; winBaseIds[n] = baseId; n++; + } + } + foreach (var kvp in decoys) + { + if (!targets.ContainsKey(kvp.Key)) + { + winIdx[n] = kvp.Value.Key; winScores[n] = kvp.Value.Value; + winDecoy[n] = true; winBaseIds[n] = kvp.Key; n++; + } + } + + // Sort: score desc, then baseId asc. Build index permutation + // then permute the parallel arrays. Sorting an int[] of + // length n with a comparison delegate beats the previous + // List>.Sort because no per-element boxing was + // required to populate the list. + var perm = new int[n]; + for (int i = 0; i < n; i++) perm[i] = i; + // The tie-break key (winBaseIds) is unique per row -- post-deduplication + // best-per-precursor selection above guarantees one row per (base_id, isDecoy) + // tuple -- so the comparator never returns 0 for distinct rows and introsort's + // instability is moot. Exemption comment must be on the Array.Sort line itself + // for the regex in CodeInspectionTest.TestNoUnstableArraySort to recognize it. + Array.Sort(perm, (a, b) => // Array.Sort OK: unique baseId tie-break makes comparator total + { + int cmp = winScores[b].CompareTo(winScores[a]); + if (cmp != 0) return cmp; + return winBaseIds[a].CompareTo(winBaseIds[b]); + }); + + // Apply permutation in-place via scratch swap arrays. Reuse + // the still-spare prefix of CountPassingQvalues as a double + // swap buffer; for int and bool we need small temp arrays. + var tmpIdx = new int[n]; + var tmpScores = new double[n]; + var tmpDecoy = new bool[n]; + for (int i = 0; i < n; i++) + { + tmpIdx[i] = winIdx[perm[i]]; + tmpScores[i] = winScores[perm[i]]; + tmpDecoy[i] = winDecoy[perm[i]]; + } + Array.Copy(tmpIdx, winIdx, n); + Array.Copy(tmpScores, winScores, n); + Array.Copy(tmpDecoy, winDecoy, n); + return n; + } + + /// + /// Variant of that operates on the + /// active prefix [0..n) of pre-allocated arrays. + /// + private static void ComputeQvaluesInto( + double[] scores, bool[] isDecoy, double[] qValuesOut, int n) + { + // Bit-identical to ComputeQvalues but operates on prefix [0..n). + int nTarget = 0; + int nDecoy = 0; + for (int i = 0; i < n; i++) + { + if (isDecoy[i]) nDecoy++; + else nTarget++; + qValuesOut[i] = nTarget > 0 ? (double)nDecoy / nTarget : 1.0; + } + double qMin = 1.0; + for (int i = n - 1; i >= 0; i--) + { + qMin = Math.Min(qMin, qValuesOut[i]); + qValuesOut[i] = qMin; } - return count; } /// @@ -1219,13 +1494,30 @@ private static void FindBestInitialFeature( private static double GridSearchC( Matrix features, bool[] labels, uint[] entryIds, double[] cValues, int[] foldAssignments, int nFolds, - ulong seed, double fdrThreshold) + ulong seed, double fdrThreshold, + SvmTrainScratchPool svmScratchPool) { - double bestC = cValues[0]; - int bestTotal = 0; - - foreach (double c in cValues) - { + // Evaluate each candidate C in parallel. Mirrors Rust's + // c_values.par_iter() in osprey-ml/src/svm.rs::grid_search_c. + // Each C is independent (no shared mutable state during + // training); the per-C totalPassing is stored by index so + // the tie-break below is deterministic. OspreyParallel.For + // (explicit threads) replaces TPL Parallel.For for the same + // reason as the outer loop above. + var totalPassingByC = new int[cValues.Length]; + OspreyParallel.For(0, cValues.Length, cValues.Length, ci => + { + // Rent one scratch per parallel c-value; reused across + // the inner sequential nFolds Train calls. Returned at + // end of this parallel body. + var localScratch = svmScratchPool != null ? svmScratchPool.Rent() : null; + // Ensure ExtractRowsInto buffers can hold the larger of + // train/test sizes (= labels.Length, the parent set, + // which is the upper bound on either subset). + if (localScratch != null) + localScratch.EnsureExtractCapacity(labels.Length, features.Cols); + try { + double c = cValues[ci]; int totalPassing = 0; for (int fold = 0; fold < nFolds; fold++) { @@ -1242,13 +1534,22 @@ private static double GridSearchC( if (trainIdx.Count == 0 || testIdx.Count == 0) continue; - var trainFeatures = ExtractRows(features, trainIdx.ToArray()); + Matrix trainFeatures, testFeatures; + if (localScratch != null) + { + trainFeatures = ExtractRowsInto(features, trainIdx.ToArray(), localScratch.TrainData); + testFeatures = ExtractRowsInto(features, testIdx.ToArray(), localScratch.TestData); + } + else + { + trainFeatures = ExtractRows(features, trainIdx.ToArray()); + testFeatures = ExtractRows(features, testIdx.ToArray()); + } var trainLabels = new bool[trainIdx.Count]; for (int i = 0; i < trainIdx.Count; i++) trainLabels[i] = labels[trainIdx[i]]; - var model = LinearSvmClassifier.Train(trainFeatures, trainLabels, c, seed); - var testFeatures = ExtractRows(features, testIdx.ToArray()); + var model = LinearSvmClassifier.Train(trainFeatures, trainLabels, c, seed, localScratch); var testScores = model.DecisionFunction(testFeatures); var testLabels = new bool[testIdx.Count]; var testEntryIds = new uint[testIdx.Count]; @@ -1258,16 +1559,28 @@ private static double GridSearchC( testEntryIds[i] = entryIds[testIdx[i]]; } - totalPassing += CountPassing(testScores, testLabels, testEntryIds, fdrThreshold); + totalPassing += CountPassing(testScores, testLabels, testEntryIds, fdrThreshold, localScratch); + } + totalPassingByC[ci] = totalPassing; + } finally { + if (localScratch != null && svmScratchPool != null) + svmScratchPool.Return(localScratch); } + }); - if (totalPassing > bestTotal) + // Tie-break: first index with the maximum totalPassing wins, + // matching the strict `>` semantics of the prior serial loop + // and the corresponding Rust path. + double bestC = cValues[0]; + int bestTotal = totalPassingByC[0]; + for (int ci = 1; ci < cValues.Length; ci++) + { + if (totalPassingByC[ci] > bestTotal) { - bestTotal = totalPassing; - bestC = c; + bestTotal = totalPassingByC[ci]; + bestC = cValues[ci]; } } - return bestC; } @@ -1508,7 +1821,7 @@ private static double[] ComputeExperimentPrecursorQvalues( // Propagate the winner's q-value to all observations sharing the // same base_id (both target and decoy sides). Matches Rust's // base_id_exp_prec_q HashMap at osprey-fdr/src/percolator.rs:2168 - // — without this, non-winning per-file observations of a + // -- without this, non-winning per-file observations of a // multi-file precursor stay at q=1.0 and downstream stages that // gate on experiment_precursor_qvalue (Stage 6 calibration refit // and reconciliation) miss the bulk of the consensus pool. @@ -1833,7 +2146,20 @@ private static void WriteStage5SubsampleDump( var order = new int[n]; for (int i = 0; i < n; i++) order[i] = i; - Array.Sort(order, (a, b) => entries[a].EntryId.CompareTo(entries[b].EntryId)); // Array.Sort OK: EntryId is unique per entry, so no ties + // EntryId is NOT unique in the 2nd-pass entries[] vector -- + // a single (base_id, charge) precursor observed across N + // files contributes N entries with the same EntryId, and + // post-reconciliation gap-fill can add yet more duplicates + // at the same (EntryId, Charge, ScanNumber). Tie-break on + // the input index a (native_position) so the dump order + // is deterministic AND matches Rust's stable + // sort_by_key(|&i| entries[i].entry_id), which preserves + // native_position order at duplicate EntryIds. + Array.Sort(order, (a, b) => // Array.Sort OK: tie-break on native_position (the input index a/b) makes the comparator total + { + int c = entries[a].EntryId.CompareTo(entries[b].EntryId); + return c != 0 ? c : a.CompareTo(b); + }); using (var sw = new StreamWriter(path)) { @@ -1939,6 +2265,59 @@ private static void WriteStage5StandardizerDump( Console.Error.WriteLine(@"Wrote Stage 5 standardizer dump: {0} ({1} features)", path, means.Length); } + /// + /// One-shot diagnostic dump of the raw per-entry feature vectors + /// fed into FeatureStandardizer.FitTransform. Mirrors Rust + /// dump_stage5_perc_input. Writes cs_stage5_perc_input.tsv with + /// columns native_position, entry_id, is_decoy, <features...> + /// sorted by (entry_id, native_position). + /// + private static void WriteStage5PercInputDump( + IList entries, + string[] featureNames) + { + const string path = @"cs_stage5_perc_input.tsv"; + var inv = CultureInfo.InvariantCulture; + int nFeatures = entries.Count > 0 ? entries[0].Features.Length : 0; + using (var sw = new StreamWriter(path)) + { + sw.NewLine = "\n"; + sw.Write(@"native_position entry_id is_decoy"); + for (int i = 0; i < nFeatures; i++) + { + string name = (featureNames != null && i < featureNames.Length) + ? featureNames[i] + : @"unknown"; + sw.Write('\t'); sw.Write(name); + } + sw.WriteLine(); + + int n = entries.Count; + int[] order = new int[n]; + for (int i = 0; i < n; i++) order[i] = i; + Array.Sort(order, (a, b) => // Array.Sort OK: tie-break on native_position (the input index a/b) makes the comparator total + { + int c = entries[a].EntryId.CompareTo(entries[b].EntryId); + return c != 0 ? c : a.CompareTo(b); + }); + + foreach (int idx in order) + { + var e = entries[idx]; + sw.Write(idx.ToString(inv)); + sw.Write('\t'); sw.Write(e.EntryId.ToString(inv)); + sw.Write('\t'); sw.Write(e.IsDecoy ? @"true" : @"false"); + for (int i = 0; i < e.Features.Length; i++) + { + sw.Write('\t'); + sw.Write(Diagnostics.FormatF64Roundtrip(e.Features[i])); + } + sw.WriteLine(); + } + } + Console.Error.WriteLine(@"Wrote Stage 5 Percolator input dump: {0} ({1} rows)", path, entries.Count); + } + // ============================================================ // Utility // ============================================================ @@ -1960,5 +2339,36 @@ private static Matrix ExtractRows(Matrix matrix, int[] rowIndices) } return Matrix.WrapNoClone(data, nRows, nCols); } + + /// + /// Variant of that writes into a + /// caller-supplied buffer (must be + /// at least rowIndices.Length * matrix.Cols long) and + /// wraps the prefix as a Matrix. Avoids the ~8 MB LOH allocation + /// per call on HRAM Astral. The trailing unused suffix of + /// is left untouched (Matrix.Rows + /// hides it). + /// + private static Matrix ExtractRowsInto(Matrix matrix, int[] rowIndices, double[] destData) + { + int nCols = matrix.Cols; + int nRows = rowIndices.Length; + int need = nRows * nCols; + if (destData.Length < need) + throw new ArgumentException( + string.Format("destData length {0} < required {1}", destData.Length, need)); + double[] src = matrix.Data; + for (int i = 0; i < nRows; i++) + { + int srcOffset = rowIndices[i] * nCols; + int dstOffset = i * nCols; + Array.Copy(src, srcOffset, destData, dstOffset, nCols); + } + // Pool-friendly wrap: Matrix.WrapPrefixNoClone accepts a + // backing array >= rows*cols. The trailing suffix of + // destData (from prior larger calls) is left untouched and + // never read. + return Matrix.WrapPrefixNoClone(destData, nRows, nCols); + } } } diff --git a/pwiz_tools/OspreySharp/OspreySharp.FDR/ProteinFdr.cs b/pwiz_tools/OspreySharp/OspreySharp.FDR/ProteinFdr.cs index cc3c23b6ff9..554e05f4524 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.FDR/ProteinFdr.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.FDR/ProteinFdr.cs @@ -33,6 +33,7 @@ using System; using System.Collections.Generic; +using System.Linq; using pwiz.OspreySharp.Core; namespace pwiz.OspreySharp.FDR @@ -445,36 +446,50 @@ public static ProteinFdrResult ComputeProteinFdr( } // Step 2: Pair picking. Iterate parsimony.Groups in deterministic - // order. Each group yields one winner: target if t >= d, else decoy. + // order. Each group yields one winner: target if t >= d, else + // decoy. Carry a sorted-accessions string on each winner so the + // Step 3 sort tiebreak can use a cross-impl-deterministic key. + // The numeric GroupId is HashMap-iteration-order-derived in + // BuildProteinParsimony and so picks different positions on + // ties cross-impl once upstream arithmetic is bit-equal. var winners = new List(); foreach (var group in parsimony.Groups) { + // Build sort_key once per group: accessions sorted then + // joined with semicolons (matches the Rust port and the + // Stage 7 diagnostic dump). + var sortedAccs = new List(group.Accessions); + sortedAccs.Sort(StringComparer.Ordinal); + string sortKey = string.Join(";", sortedAccs); + bool hasT = targetScore.TryGetValue(group.Id, out double t); bool hasD = decoyScore.TryGetValue(group.Id, out double d); if (hasT && hasD) { if (t >= d) - winners.Add(new ProteinWinner { GroupId = group.Id, Score = t, IsDecoy = false }); + winners.Add(new ProteinWinner { GroupId = group.Id, SortKey = sortKey, Score = t, IsDecoy = false }); else - winners.Add(new ProteinWinner { GroupId = group.Id, Score = d, IsDecoy = true }); + winners.Add(new ProteinWinner { GroupId = group.Id, SortKey = sortKey, Score = d, IsDecoy = true }); } else if (hasT) { - winners.Add(new ProteinWinner { GroupId = group.Id, Score = t, IsDecoy = false }); + winners.Add(new ProteinWinner { GroupId = group.Id, SortKey = sortKey, Score = t, IsDecoy = false }); } else if (hasD) { - winners.Add(new ProteinWinner { GroupId = group.Id, Score = d, IsDecoy = true }); + winners.Add(new ProteinWinner { GroupId = group.Id, SortKey = sortKey, Score = d, IsDecoy = true }); } } // Step 3: Cumulative FDR. Sort winners by score descending, - // tiebreak by group_id ascending for determinism. - winners.Sort((a, b) => + // tiebreak by sorted accessions ASCENDING — cross-impl- + // deterministic, unlike GroupId which is HashMap-iteration- + // order from BuildProteinParsimony. + winners.Sort((a, b) => // Array.Sort OK: SortKey is the sorted-accessions string from BuildProteinParsimony, which assigns a unique accessions list to each ProteinGroup (identical sets are merged), so the comparator never returns 0 and unstable-sort tie reorder cannot fire { int cmp = b.Score.CompareTo(a.Score); if (cmp != 0) return cmp; - return a.GroupId.CompareTo(b.GroupId); + return string.CompareOrdinal(a.SortKey, b.SortKey); }); var rawQvalues = new double[winners.Count]; @@ -495,11 +510,13 @@ public static ProteinFdrResult ComputeProteinFdr( // to GroupQvalues / GroupScores. var groupQvalues = new Dictionary(); var groupScores = new Dictionary(); + var monotonicQvalues = new double[winners.Count]; double minQ = 1.0; for (int i = winners.Count - 1; i >= 0; i--) { if (rawQvalues[i] < minQ) minQ = rawQvalues[i]; + monotonicQvalues[i] = minQ; var w = winners[i]; if (!w.IsDecoy) { @@ -508,6 +525,18 @@ public static ProteinFdrResult ComputeProteinFdr( } } + // Cross-impl bisection dump (env-var-gated, no-op in production). + // The flag check short-circuits the LINQ projection below; in + // the disabled-dump path this whole block is one field read. + // Dump function lives in FdrDiagnostics so the file I/O stays + // isolated from the protein-FDR algorithm code. + if (FdrDiagnostics.DumpStage7Winners) + { + FdrDiagnostics.WriteStage7WinnersDump( + winners.Select(w => (w.Score, w.IsDecoy)).ToList(), + rawQvalues, monotonicQvalues); + } + // Step 5: Propagate to peptides. Each peptide's q-value is the min // (best) q-value across its groups. Peptides whose only groups lost // the pair stay at q = 1.0. @@ -537,6 +566,7 @@ public static ProteinFdrResult ComputeProteinFdr( private struct ProteinWinner { public uint GroupId; + public string SortKey; // sorted-accessions string for cross-impl tiebreak public double Score; public bool IsDecoy; } @@ -576,6 +606,12 @@ public static Dictionary CollectBestPeptideScores( } } } + + // Cross-impl bisection dump (env-var-gated, no-op in production). + // See FdrDiagnostics.WriteBestPeptideScoresDump for context. + if (FdrDiagnostics.DumpBestPeptideScores) + FdrDiagnostics.WriteBestPeptideScoresDump(best); + return best; } @@ -602,5 +638,51 @@ public static void PropagateProteinQvalues( } } } + + /// + /// First-pass protein FDR: build parsimony from peptides passing peptide-level + /// run FDR, run picked-protein FDR at (1x Savitski + /// gate), and write the resulting q-values into + /// on every stub. Mirrors Rust pipeline.rs::run_analysis first-pass block + /// (around line 4292). Caller is responsible for any logging, dump diagnostics, + /// and downstream consumption. + /// + /// Used by FirstJoinTask for the in-process pipeline (runs after first-pass FDR, + /// before compaction) and by PerFileRescoreTask for the --join-at-pass=2 + /// rehydration path (runs after sidecar load, before compaction) so the protein- + /// rescue branch of compaction has fresh RunProteinQvalue values matching + /// what Rust computes inline. Without it, the rehydrated C# pipeline used only + /// the RunProteinQvalue values stored in the 1st-pass FDR sidecar; for + /// single-file --join-at-pass=2 runs that left 19 peptides outside Rust's + /// post-compaction detected set on Stellar Single, causing a 1-protein delta in + /// Stage 7 picked-protein output. + /// + public static void RunFirstPassProteinFdr( + IList>> perFileEntries, + IList fullLibrary, + OspreyConfig config) + { + // Detected-peptide gate: targets passing peptide-level run FDR. + // Matches Rust pipeline.rs:4301 (e.run_peptide_qvalue <= config.run_fdr). + var detectedPeptides = new HashSet(StringComparer.Ordinal); + foreach (var kvp in perFileEntries) + { + foreach (var entry in kvp.Value) + { + if (!entry.IsDecoy && entry.RunPeptideQvalue <= config.RunFdr) + detectedPeptides.Add(entry.ModifiedSequence); + } + } + + var parsimony = BuildProteinParsimony( + fullLibrary, config.SharedPeptides, detectedPeptides); + var bestScores = CollectBestPeptideScores(perFileEntries); + var proteinFdr = ComputeProteinFdr(parsimony, bestScores, config.RunFdr); + + // Set RunProteinQvalue ONLY. ExperimentProteinQvalue is set by the + // post-output Stage 7 second-pass protein FDR (Rust's second-pass). + PropagateProteinQvalues(perFileEntries, proteinFdr, + setRun: true, setExperiment: false); + } } } diff --git a/pwiz_tools/OspreySharp/OspreySharp.FDR/Reconciliation/ConsensusRts.cs b/pwiz_tools/OspreySharp/OspreySharp.FDR/Reconciliation/ConsensusRts.cs index 0e7e4cbdbb3..914abafd5a9 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.FDR/Reconciliation/ConsensusRts.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.FDR/Reconciliation/ConsensusRts.cs @@ -36,8 +36,6 @@ namespace pwiz.OspreySharp.FDR.Reconciliation /// public static class ConsensusRts { - private const string DECOY_PREFIX = @"DECOY_"; - /// /// For each target peptide passing at the /// run-precursor level (hard gate), and its paired decoy, computes a @@ -83,30 +81,38 @@ public static IReadOnlyList Compute( // 1. Collect target peptides passing the run-level hard gate // (or rescued by protein FDR for peptide-level borderline cases). + // We also record the set of qualifying target *base_ids* so that + // paired decoys can be identified by base_id linkage + // (entry_id & 0x7FFFFFFF) rather than by stripping a "DECOY_" + // prefix from modified_sequence. The prefix-strip approach only + // works for Osprey-generated decoys; it silently misses library- + // supplied decoys (Carafe etc.) whose modified sequence carries + // no prefix. Pairing was already established by the FDRBench + // manifest or composition fallback during library load. Mirrors + // Rust reconciliation.rs::compute_consensus_rts. var targetPeptides = new HashSet(StringComparer.Ordinal); + var targetBaseIds = new HashSet(); foreach (var kvp in perFileEntries) { foreach (var entry in kvp.Value) { if (Qualifies(entry, consensusFdr, proteinFdrThreshold)) + { targetPeptides.Add(entry.ModifiedSequence); + targetBaseIds.Add(entry.EntryId & 0x7FFFFFFFu); + } } } if (targetPeptides.Count == 0) return Array.Empty(); - // 2. Collect paired decoy peptides (DECOY_). + // 2. Collect paired decoy peptides via base_id linkage. var decoyPeptides = new HashSet(StringComparer.Ordinal); foreach (var kvp in perFileEntries) { foreach (var entry in kvp.Value) { - if (!entry.IsDecoy) - continue; - var targetSeq = entry.ModifiedSequence.StartsWith(DECOY_PREFIX, StringComparison.Ordinal) - ? entry.ModifiedSequence.Substring(DECOY_PREFIX.Length) - : entry.ModifiedSequence; - if (targetPeptides.Contains(targetSeq)) + if (entry.IsDecoy && targetBaseIds.Contains(entry.EntryId & 0x7FFFFFFFu)) decoyPeptides.Add(entry.ModifiedSequence); } } diff --git a/pwiz_tools/OspreySharp/OspreySharp.IO/DiannTsvLoader.cs b/pwiz_tools/OspreySharp/OspreySharp.IO/DiannTsvLoader.cs index e43e47e7aa5..bbca3351aa2 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.IO/DiannTsvLoader.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.IO/DiannTsvLoader.cs @@ -26,6 +26,7 @@ using System.Globalization; using System.IO; using System.Text; +using System.Xml; using pwiz.OspreySharp.Core; namespace pwiz.OspreySharp.IO @@ -508,18 +509,43 @@ private static string GetFieldOrNull(string[] fields, int index) private static double ParseDouble(string s, string name, int rowNum) { - double value; - if (!double.TryParse(s, NumberStyles.Float, CultureInfo.InvariantCulture, out value)) + // XmlConvert.ToDouble is IEEE-754 correct (XML schema spec requires + // correct rounding) whereas .NET Framework 4.7.2's double.TryParse + // can be off by a few ULPs on 16-digit scientific values. For TSV + // library values like "400.1277954005887", .NET Framework parses + // to a different f64 than the IEEE-correct round-to-nearest-even, + // producing cross-impl drift in mz_min/mz_max and downstream bin + // widths. Rust's `str::parse::` is IEEE-correct, so using + // XmlConvert brings the two parsers into bit-for-bit agreement. + try + { + return XmlConvert.ToDouble(s); + } + catch (FormatException) + { throw new InvalidDataException(string.Format("Invalid {0} '{1}' at row {2}", name, s, rowNum)); - return value; + } + catch (OverflowException) + { + throw new InvalidDataException(string.Format("Invalid {0} '{1}' at row {2}", name, s, rowNum)); + } } private static float ParseFloat(string s, string name, int rowNum) { - float value; - if (!float.TryParse(s, NumberStyles.Float, CultureInfo.InvariantCulture, out value)) + // XmlConvert for IEEE-754 correct parsing - see ParseDouble note. + try + { + return XmlConvert.ToSingle(s); + } + catch (FormatException) + { throw new InvalidDataException(string.Format("Invalid {0} '{1}' at row {2}", name, s, rowNum)); - return value; + } + catch (OverflowException) + { + throw new InvalidDataException(string.Format("Invalid {0} '{1}' at row {2}", name, s, rowNum)); + } } private static byte ParseByte(string s, string name, int rowNum) @@ -534,10 +560,19 @@ private static double ParseDoubleOrDefault(string s, double defaultValue) { if (string.IsNullOrEmpty(s)) return defaultValue; - double value; - if (double.TryParse(s, NumberStyles.Float, CultureInfo.InvariantCulture, out value)) - return value; - return defaultValue; + // XmlConvert for IEEE-754 correct parsing - see ParseDouble note. + try + { + return XmlConvert.ToDouble(s); + } + catch (FormatException) + { + return defaultValue; + } + catch (OverflowException) + { + return defaultValue; + } } #endregion diff --git a/pwiz_tools/OspreySharp/OspreySharp.IO/FileSaver.cs b/pwiz_tools/OspreySharp/OspreySharp.IO/FileSaver.cs index 356673454dc..ed6d6e83666 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.IO/FileSaver.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.IO/FileSaver.cs @@ -61,7 +61,7 @@ public sealed class FileSaver : IDisposable /// Resolves against the current /// working directory so a bare-filename argument (no directory /// component) lands the temp file alongside its destination - /// rather than failing inside + /// rather than failing inside Path.GetDirectoryName /// downstream. Diverges from the SharedBatch original which /// happens to never see relative paths. /// diff --git a/pwiz_tools/OspreySharp/OspreySharp.IO/MzmlReader.cs b/pwiz_tools/OspreySharp/OspreySharp.IO/MzmlReader.cs index 99444d5cc00..1a18a55bbf7 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.IO/MzmlReader.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.IO/MzmlReader.cs @@ -60,8 +60,6 @@ public class MzmlReader private const string CV_ZLIB_COMPRESSION = "MS:1000574"; private const string CV_NO_COMPRESSION = "MS:1000576"; - private const double DEFAULT_ISOLATION_HALF_WIDTH = 12.5; - #endregion /// @@ -238,33 +236,27 @@ private static void ParseSpectrumElement(XmlReader reader, uint spectrumIndex, else if (currentContext == "selectedIon") { if (accession == CV_SELECTED_ION_MZ && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out precursorMz); + TryParseXmlDouble(value, out precursorMz); else if (accession == CV_PEAK_INTENSITY && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out precursorIntensity); + TryParseXmlDouble(value, out precursorIntensity); } else if (currentContext == "isolationWindow") { if (accession == CV_ISOLATION_WINDOW_TARGET && value != null) { - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out isoTarget); + TryParseXmlDouble(value, out isoTarget); hasIsolationWindow = true; } else if (accession == CV_ISOLATION_WINDOW_LOWER && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out isoLower); + TryParseXmlDouble(value, out isoLower); else if (accession == CV_ISOLATION_WINDOW_UPPER && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out isoUpper); + TryParseXmlDouble(value, out isoUpper); } else if (currentContext == "scan") { if (accession == CV_SCAN_START_TIME_MINUTES && value != null) { - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out retentionTime); + TryParseXmlDouble(value, out retentionTime); // Check unitName for seconds string unitName = subtree.GetAttribute("unitName"); if (unitName != null && unitName.IndexOf("second", StringComparison.OrdinalIgnoreCase) >= 0) @@ -273,8 +265,7 @@ private static void ParseSpectrumElement(XmlReader reader, uint spectrumIndex, else if (accession == CV_RETENTION_TIME_SECONDS && value != null) { double seconds; - if (double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out seconds)) + if (TryParseXmlDouble(value, out seconds)) retentionTime = seconds / 60.0; } else if (accession == CV_MS_LEVEL && value != null) @@ -362,10 +353,28 @@ private static void ParseSpectrumElement(XmlReader reader, uint spectrumIndex, if (center <= 0) return; // Skip spectra without precursor info - double lowerOffset = isoLower > 0 ? isoLower : DEFAULT_ISOLATION_HALF_WIDTH; - double upperOffset = isoUpper > 0 ? isoUpper : DEFAULT_ISOLATION_HALF_WIDTH; + // Fail fast on missing offsets rather than substituting a + // 12.5 hardcoded default. DIA processing cannot proceed + // without true isolation windows, and a silent default + // produces bogus results that are very hard to diagnose + // downstream. The error names the spectrum index and which + // cvParam is missing. Mirrors the equivalent fail-fast + // change in osprey/crates/osprey-io/src/mzml/parser.rs + // (PR #39 on maccoss/osprey). + if (isoLower <= 0) + throw new InvalidDataException(string.Format( + "spectrum index {0}: no valid isolation-window lower offset " + + "(cvParam MS:1000828 missing or non-positive); cannot process DIA data " + + "without true isolation windows.", + spectrumIndex)); + if (isoUpper <= 0) + throw new InvalidDataException(string.Format( + "spectrum index {0}: no valid isolation-window upper offset " + + "(cvParam MS:1000829 missing or non-positive); cannot process DIA data " + + "without true isolation windows.", + spectrumIndex)); - var isoWindow = new IsolationWindow(center, lowerOffset, upperOffset); + var isoWindow = new IsolationWindow(center, isoLower, isoUpper); ms2List.Add(new Spectrum { @@ -520,30 +529,25 @@ private static RawSpectrumData ParseSpectrumRaw(XmlReader reader, uint spectrumI else if (currentContext == "selectedIon") { if (accession == CV_SELECTED_ION_MZ && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out raw.PrecursorMz); + TryParseXmlDouble(value, out raw.PrecursorMz); } else if (currentContext == "isolationWindow") { if (accession == CV_ISOLATION_WINDOW_TARGET && value != null) { - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out raw.IsoTarget); + TryParseXmlDouble(value, out raw.IsoTarget); raw.HasIsolationWindow = true; } else if (accession == CV_ISOLATION_WINDOW_LOWER && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out raw.IsoLower); + TryParseXmlDouble(value, out raw.IsoLower); else if (accession == CV_ISOLATION_WINDOW_UPPER && value != null) - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out raw.IsoUpper); + TryParseXmlDouble(value, out raw.IsoUpper); } else if (currentContext == "scan") { if (accession == CV_SCAN_START_TIME_MINUTES && value != null) { - double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out raw.RetentionTime); + TryParseXmlDouble(value, out raw.RetentionTime); string unitName = subtree.GetAttribute("unitName"); if (unitName != null && unitName.IndexOf("second", StringComparison.OrdinalIgnoreCase) >= 0) raw.RetentionTime /= 60.0; @@ -551,8 +555,7 @@ private static RawSpectrumData ParseSpectrumRaw(XmlReader reader, uint spectrumI else if (accession == CV_RETENTION_TIME_SECONDS && value != null) { double seconds; - if (double.TryParse(value, System.Globalization.NumberStyles.Float, - System.Globalization.CultureInfo.InvariantCulture, out seconds)) + if (TryParseXmlDouble(value, out seconds)) raw.RetentionTime = seconds / 60.0; } else if (accession == CV_MS_LEVEL && value != null) @@ -774,15 +777,27 @@ public Spectrum ToMs2Spectrum() if (center <= 0) return null; - double lowerOffset = IsoLower > 0 ? IsoLower : DEFAULT_ISOLATION_HALF_WIDTH; - double upperOffset = IsoUpper > 0 ? IsoUpper : DEFAULT_ISOLATION_HALF_WIDTH; + // Fail fast on missing offsets (see the equivalent block + // in the linear convert path above for the rationale). + if (IsoLower <= 0) + throw new InvalidDataException(string.Format( + "spectrum index {0}: no valid isolation-window lower offset " + + "(cvParam MS:1000828 missing or non-positive); cannot process DIA data " + + "without true isolation windows.", + Index)); + if (IsoUpper <= 0) + throw new InvalidDataException(string.Format( + "spectrum index {0}: no valid isolation-window upper offset " + + "(cvParam MS:1000829 missing or non-positive); cannot process DIA data " + + "without true isolation windows.", + Index)); return new Spectrum { ScanNumber = Index, RetentionTime = RetentionTime, PrecursorMz = PrecursorMz > 0 ? PrecursorMz : center, - IsolationWindow = new IsolationWindow(center, lowerOffset, upperOffset), + IsolationWindow = new IsolationWindow(center, IsoLower, IsoUpper), Mzs = MzArray, Intensities = IntensityArray, }; @@ -790,6 +805,29 @@ public Spectrum ToMs2Spectrum() } #endregion + + /// + /// IEEE-754 correct double parser for mzML cvParam values. + /// .NET Framework 4.7.2's double.TryParse can be off by 1-2 ULPs + /// on 16-digit scientific values; XmlConvert.ToDouble is required + /// by XML schema spec to be IEEE-correct, matching Rust mzdata's + /// parser. Without this, cvParams like scan_start_time + /// "17.0286203070330" parse to slightly different f64 cross-impl, + /// producing 2-ULP apex_rt drift downstream. + /// + private static bool TryParseXmlDouble(string s, out double result) + { + try + { + result = XmlConvert.ToDouble(s); + return true; + } + catch + { + result = 0.0; + return false; + } + } } /// diff --git a/pwiz_tools/OspreySharp/OspreySharp.IO/ParquetScoreCache.cs b/pwiz_tools/OspreySharp/OspreySharp.IO/ParquetScoreCache.cs index 87c18e1be29..103bf731a60 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.IO/ParquetScoreCache.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.IO/ParquetScoreCache.cs @@ -217,9 +217,21 @@ public static void WriteScoresParquet(string path, List en for (int f = 0; f < NUM_PIN_FEATURES; f++) featureArrays[f] = new double[n]; + // Iterate in canonical sorted order (entry_id, charge, scan_number) + // so per-side parquets have identical physical row layout across the + // 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. Mirrors Rust + // pipeline.rs::write_scores_parquet_with_metadata. + var sortedIndices = Enumerable.Range(0, n) + .OrderBy(idx => entries[idx].EntryId) + .ThenBy(idx => entries[idx].Charge) + .ThenBy(idx => entries[idx].ScanNumber) + .ToArray(); + for (int i = 0; i < n; i++) { - var entry = entries[i]; + var entry = entries[sortedIndices[i]]; entryIds[i] = entry.EntryId; isDecoys[i] = entry.IsDecoy; sequences[i] = entry.Sequence ?? string.Empty; @@ -340,9 +352,41 @@ public static void WriteScoresParquet(string path, List entries, for (int f = 0; f < NUM_PIN_FEATURES; f++) featureArrays[f] = new double[n]; + // 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. Mirrors Rust + // pipeline.rs::write_scores_parquet_with_metadata. ParquetIndex is + // assigned to the post-sort destination row below. + var sortedIndices = Enumerable.Range(0, n) + .OrderBy(idx => entries[idx].EntryId) + .ThenBy(idx => entries[idx].Charge) + .ThenBy(idx => entries[idx].ScanNumber) + .ToArray(); + for (int i = 0; i < n; i++) { - var entry = entries[i]; + var entry = entries[sortedIndices[i]]; + // Assign ParquetIndex to match the row position we are + // about to write. Mirrors LoadFdrStubsFromParquet, which + // assigns ParquetIndex = row on read. Without this + // assignment, in-memory entries reach Stage 5 + // ReconciliationPlanner with ParquetIndex = 0 (FdrEntry + // default), and every entry's per-file CWT lookup + // (fileCwt[entry.ParquetIndex]) grabs the first row's + // CwtCandidate list instead of its own -- the planner + // then force-integrates almost every entry because the + // wrong CWT list has no candidate near the expected RT. + // The HPC chain path was unaffected because its entries + // are reloaded via LoadFdrStubsFromParquet, which sets + // ParquetIndex correctly. Found by C# in-memory vs + // C# HPC-chain strict-rehydration bisection on Stellar + // (Stage 5 boundary check: .1st-pass.fdr_scores.bin + // byte-identical but reconciliation.json action shape + // diverged -- 35K use_cwt actions on HPC side, 814 on + // in-memory side, total identical). + entry.ParquetIndex = (uint)i; entryIds[i] = entry.EntryId; isDecoys[i] = entry.IsDecoy; charges[i] = entry.Charge; diff --git a/pwiz_tools/OspreySharp/OspreySharp.IO/ReconciliationFile.cs b/pwiz_tools/OspreySharp/OspreySharp.IO/ReconciliationFile.cs index 6e42a630c30..638b29ee3fe 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.IO/ReconciliationFile.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.IO/ReconciliationFile.cs @@ -50,28 +50,43 @@ namespace pwiz.OspreySharp.IO /// public class ReconciliationFile { - /// Current schema version. Bump on incompatible changes. - public const int CurrentFormatVersion = 1; + /// + /// Current schema version. Bump on incompatible changes. + /// + /// v1: initial format. + /// v2: added file_stems so per-file Stage 6 rescore workers + /// can compute the reconciliation parameter hash that the + /// downstream --join-at-pass=2 merge node expects (the + /// hash is computed over all files in the join, not the + /// worker's single parquet). Old v1 files deserialize with an + /// empty list; the worker falls back + /// to its OspreyConfig.InputFiles stems in that case, + /// preserving v1 behavior. + /// + public const int CurrentFormatVersion = 2; + + [JsonProperty("file_stems", Order = 0)] + public List FileStems { get; set; } - [JsonProperty("forced_integration_actions", Order = 0)] + [JsonProperty("forced_integration_actions", Order = 1)] public List ForcedIntegrationActions { get; set; } - [JsonProperty("format_version", Order = 1)] + [JsonProperty("format_version", Order = 2)] public int FormatVersion { get; set; } - [JsonProperty("gap_fill_targets", Order = 2)] + [JsonProperty("gap_fill_targets", Order = 3)] public List GapFillTargets { get; set; } - [JsonProperty("library_hash", Order = 3)] + [JsonProperty("library_hash", Order = 4)] public string LibraryHash { get; set; } - [JsonProperty("refined_rt_calibration", Order = 4, NullValueHandling = NullValueHandling.Include)] + [JsonProperty("refined_rt_calibration", Order = 5, NullValueHandling = NullValueHandling.Include)] public RefinedRtCalibrationJson RefinedRtCalibration { get; set; } - [JsonProperty("search_hash", Order = 5)] + [JsonProperty("search_hash", Order = 6)] public string SearchHash { get; set; } - [JsonProperty("use_cwt_peak_actions", Order = 6)] + [JsonProperty("use_cwt_peak_actions", Order = 7)] public List UseCwtPeakActions { get; set; } /// @@ -96,6 +111,22 @@ public static ReconciliationFile Load(string path) "Reconciliation file {0} has unsupported format_version {1} (expected {2})", path, parsed.FormatVersion, CurrentFormatVersion)); } + // v2 envelopes must carry the planner's full join file_stems set; + // a deserialized v2 file with file_stems missing or empty would + // silently flow through RescoreHydration with joinFileStems = [] + // and cause downstream --join-at-pass=2 to compute a single-file + // ReconciliationParameterHash for what was meant to be a + // multi-file join. JsonProperty does not enforce required, so + // assert it here. Matches the Rust serde behavior (file_stems is + // a required field, not defaulted). + if (parsed.FileStems == null || parsed.FileStems.Count == 0) + { + throw new InvalidDataException(string.Format( + "Reconciliation file {0} has format_version {1} but file_stems is missing " + + "or empty; v{1} envelopes are required to carry the planner's full join " + + "file set.", + path, CurrentFormatVersion)); + } return parsed; } diff --git a/pwiz_tools/OspreySharp/OspreySharp.ML/LinearSvmClassifier.cs b/pwiz_tools/OspreySharp/OspreySharp.ML/LinearSvmClassifier.cs index 33c173d0186..2689cbaf5b1 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.ML/LinearSvmClassifier.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.ML/LinearSvmClassifier.cs @@ -34,9 +34,230 @@ // Licensed under the MIT License using System; +using System.Collections.Concurrent; +using System.Threading; namespace pwiz.OspreySharp.ML { + /// + /// Minimal explicit-thread parallel-for, modeled after + /// pwiz.Common.SystemUtil.ParallelEx (Skyline). The key difference + /// vs System.Threading.Tasks.Parallel.For is direct + /// thread allocation: a fixed count of dedicated Threads pull + /// indices from a shared atomic counter and call the body. No + /// TaskReplicator heuristics, no ThreadPool hill-climbing, no + /// scheduler throttling. Used by the Percolator hot path where + /// TPL Parallel.For was failing to scale beyond ~2.5x on HRAM + /// Astral (vs Rust rayon's ~9x on the same workload). + /// + public static class OspreyParallel + { + public static void For(int fromInclusive, int toExclusive, int threadCount, Action body) + { + int count = toExclusive - fromInclusive; + if (count <= 0) return; + if (count == 1 || threadCount <= 1) + { + for (int i = fromInclusive; i < toExclusive; i++) + body(i); + return; + } + + int effectiveThreads = Math.Min(threadCount, count); + int next = fromInclusive - 1; // pre-decrement; Interlocked.Increment yields fromInclusive on first call + Exception firstException = null; + object exLock = new object(); + + var threads = new Thread[effectiveThreads]; + for (int t = 0; t < effectiveThreads; t++) + { + threads[t] = new Thread(() => + { + while (true) + { + int i = Interlocked.Increment(ref next); + if (i >= toExclusive) return; + try + { + body(i); + } + catch (Exception ex) + { + lock (exLock) + { + if (firstException == null) firstException = ex; + } + return; + } + } + }); + threads[t].IsBackground = true; + threads[t].Name = "OspreyParallel"; + threads[t].Start(); + } + foreach (var th in threads) th.Join(); + if (firstException != null) + throw new AggregateException("Exception in OspreyParallel.For", firstException); + } + } + + + /// + /// Reusable per-call buffers for LinearSvmClassifier.Train. + /// On large training sets each Train call allocated five fresh arrays + /// (y, diag, alpha, w, indices) totalling ~2 MB at HRAM-Astral scale + /// (n ~ 51K). With ~570 Train calls running 8-way parallel under + /// Percolator's grid search, that allocation pressure showed up as + /// per-call slowdown vs Rust (which uses stack-frame-local Vec but + /// no GC pressure). Pool one scratch per parallel worker; rent/return + /// around each Train call. + /// + public sealed class SvmTrainScratch + { + // n-sized (resized as the largest seen n grows; never shrinks) + public double[] Y; + public double[] Diag; + public double[] Alpha; + public int[] Indices; + // (p+1)-sized; p (feature count) is constant in a Percolator run + public double[] W; + + // Pooled row-major data buffers for ExtractRows results used in + // Percolator's grid search. Each one wraps an (rows * p)-size + // double[] that the caller fills via ExtractRowsInto and then + // hands as a Matrix to LinearSvm.Train / DecisionFunction. For + // HRAM Astral these would otherwise be ~8 MB LOH allocations + // ~540x per file. TrainData and TestData are paired so a single + // grid-search iteration can hold both simultaneously. + public double[] TrainData; + public double[] TestData; + + // Pooled buffers for PercolatorFdr.CountPassing's two per-call + // arrays (allIndices: 0..n-1; qValues: per-winner). Sized to + // initialN at scratch construction; EnsureCountPassingCapacity + // grows on rare oversize requests. + public int[] CountPassingIndices; + public double[] CountPassingQvalues; + + // Pooled output buffers for the hot-path CompeteFromIndicesInto + // helper. Sized to initialN. The active prefix length is + // returned by the helper; callers read only [0..count). + public int[] CompetitionWinnerIndices; + public double[] CompetitionWinnerScores; + public bool[] CompetitionWinnerIsDecoy; + + public SvmTrainScratch(int initialN, int p) + { + Y = new double[initialN]; + Diag = new double[initialN]; + Alpha = new double[initialN]; + Indices = new int[initialN]; + W = new double[p + 1]; + // Pre-allocate the ExtractRows buffers up front -- they + // dominate per-call allocation pressure (8+ MB each for HRAM + // Astral). The pool constructor knows the largest expected + // subset size (subN); sizing here avoids the first-iteration + // LOH stampede when ~20 parallel scratches each lazily + // allocate 17 MB simultaneously (showed up as 10s OwnTime + // in EnsureExtractCapacity in dotTrace). + int extractCap = initialN * p; + TrainData = new double[extractCap]; + TestData = new double[extractCap]; + CountPassingIndices = new int[initialN]; + CountPassingQvalues = new double[initialN]; + CompetitionWinnerIndices = new int[initialN]; + CompetitionWinnerScores = new double[initialN]; + CompetitionWinnerIsDecoy = new bool[initialN]; + } + + [System.Runtime.CompilerServices.MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] + public void EnsureCountPassingCapacity(int n) + { + if (CountPassingIndices.Length < n) + CountPassingIndices = new int[n]; + if (CountPassingQvalues.Length < n) + CountPassingQvalues = new double[n]; + if (CompetitionWinnerIndices.Length < n) + CompetitionWinnerIndices = new int[n]; + if (CompetitionWinnerScores.Length < n) + CompetitionWinnerScores = new double[n]; + if (CompetitionWinnerIsDecoy.Length < n) + CompetitionWinnerIsDecoy = new bool[n]; + } + + public void EnsureCapacity(int n, int p) + { + if (Y.Length < n) + { + Y = new double[n]; + Diag = new double[n]; + Alpha = new double[n]; + Indices = new int[n]; + } + if (W.Length < p + 1) + W = new double[p + 1]; + } + + /// + /// Ensure / each + /// have capacity for * + /// doubles. The constructor pre-sizes both to the expected max, + /// so this is a no-op in the steady state; the branch covers + /// rare cases where a caller's actual subset exceeds the + /// pool's initialN. + /// + [System.Runtime.CompilerServices.MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] + public void EnsureExtractCapacity(int rows, int p) + { + int need = rows * p; + if (TrainData.Length < need) + TrainData = new double[need]; + if (TestData.Length < need) + TestData = new double[need]; + } + } + + /// + /// Concurrent pool of sets. Same pattern + /// as XcorrScratchPool: organic growth up to the parallel-worker + /// high-water mark; arrays live in gen-2 LOH for the lifetime of the + /// Percolator run; no LOH allocation in steady state. + /// + public sealed class SvmTrainScratchPool + { + private readonly ConcurrentBag _bag = new ConcurrentBag(); + private readonly int _initialN; + private readonly int _p; + private int _allocCount; + + public SvmTrainScratchPool(int initialN, int p) + { + _initialN = initialN; + _p = p; + } + + public int AllocCount { get { return _allocCount; } } + + public SvmTrainScratch Rent() + { + SvmTrainScratch s; + if (_bag.TryTake(out s)) + return s; + Interlocked.Increment(ref _allocCount); + return new SvmTrainScratch(_initialN, _p); + } + + public void Return(SvmTrainScratch s) + { + if (s == null) + return; + // No zeroing on return: Train re-initializes Y, Alpha, W, + // Diag, and Indices from scratch each call. + _bag.Add(s); + } + } + + /// /// Deterministic xorshift64 PRNG for reproducible shuffling. /// Matches the Rust implementation: x ^= x << 13; x ^= x >> 7; x ^= x << 17. @@ -201,6 +422,20 @@ public LinearSvmClassifier(double[] weights, double bias) /// Random seed for reproducible shuffling /// Trained LinearSvmClassifier model public static LinearSvmClassifier Train(Matrix features, bool[] labels, double c, ulong seed) + { + return Train(features, labels, c, seed, null); + } + + /// + /// Overload that reuses pre-allocated working buffers from + /// . Pass a per-worker + /// rented from + /// to avoid the five n-sized + /// array allocations per call. Pass null to allocate fresh (the + /// pre-pool behavior, retained for tests and ad-hoc callers). + /// + public static LinearSvmClassifier Train(Matrix features, bool[] labels, double c, ulong seed, + SvmTrainScratch scratch) { if (features.Rows != labels.Length) throw new ArgumentException("Feature rows must match label count"); @@ -217,15 +452,41 @@ public static LinearSvmClassifier Train(Matrix features, bool[] labels, double c double[] data = features.Data; int cols = p; + // Working buffers: either from the scratch pool (no allocation) + // or fresh per-call (the legacy path). Pool mode reuses arrays + // sized by the largest-seen n; we re-initialize the prefix we + // actually use, so leftover bytes past n are harmless. + double[] y, diag, alpha, w; + int[] indices; + if (scratch != null) + { + scratch.EnsureCapacity(n, p); + y = scratch.Y; + diag = scratch.Diag; + alpha = scratch.Alpha; + w = scratch.W; + indices = scratch.Indices; + // alpha and w start at zero each call; the loop reads them + // before writing so we must clear the prefix we use. + Array.Clear(alpha, 0, n); + Array.Clear(w, 0, p + 1); + } + else + { + y = new double[n]; + diag = new double[n]; + alpha = new double[n]; + w = new double[p + 1]; + indices = new int[n]; + } + // Convert labels: target (false) -> +1, decoy (true) -> -1 - var y = new double[n]; for (int i = 0; i < n; i++) y[i] = labels[i] ? -1.0 : 1.0; double inv2c = 1.0 / (2.0 * c); // Precompute diagonal: D_ii = ||x_i||^2 + 1.0 (bias feature) + 1/(2C) - var diag = new double[n]; for (int i = 0; i < n; i++) { double normSq = 0.0; @@ -238,14 +499,8 @@ public static LinearSvmClassifier Train(Matrix features, bool[] labels, double c diag[i] = normSq + 1.0 + inv2c; } - // Initialize dual variables and primal weight vector - // w has p+1 elements: w[0..p] = feature weights, w[p] = bias - var alpha = new double[n]; - var w = new double[p + 1]; - // RNG for index permutation var rng = new XorShift64(seed); - var indices = new int[n]; for (int i = 0; i < n; i++) indices[i] = i; @@ -253,7 +508,7 @@ public static LinearSvmClassifier Train(Matrix features, bool[] labels, double c for (int iter = 0; iter < MAX_ITER; iter++) { - FisherYatesShuffle(indices, rng); + FisherYatesShuffle(indices, n, rng); double maxPgViolation = 0.0; @@ -349,10 +604,13 @@ public bool[] Predict(Matrix features) /// /// Fisher-Yates shuffle using XorShift64 PRNG. /// Matches the Rust implementation exactly. + /// Shuffles slice[0..length] -- callers that pass a pooled + /// over-sized buffer must supply the active prefix length so + /// shuffling stops there. /// - private static void FisherYatesShuffle(int[] slice, XorShift64 rng) + private static void FisherYatesShuffle(int[] slice, int length, XorShift64 rng) { - for (int i = slice.Length - 1; i >= 1; i--) + for (int i = length - 1; i >= 1; i--) { int j = (int)(rng.Next() % (ulong)(i + 1)); int tmp = slice[i]; diff --git a/pwiz_tools/OspreySharp/OspreySharp.ML/Matrix.cs b/pwiz_tools/OspreySharp/OspreySharp.ML/Matrix.cs index 48a2f97739c..ca8352e8bf3 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.ML/Matrix.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.ML/Matrix.cs @@ -73,6 +73,24 @@ internal static Matrix WrapNoClone(double[] data, int rows, int cols) return new Matrix(data, rows, cols, takeOwnership: true); } + /// + /// Pool-friendly wrap: accepts a backing array that is at least + /// rows * cols elements long. Only the first + /// rows * cols cells are read; the trailing suffix + /// (from a larger pooled buffer) is ignored. Use this when + /// renting an over-sized scratch double[] and only filling a + /// prefix this call. + /// + internal static Matrix WrapPrefixNoClone(double[] data, int rows, int cols) + { + int need = rows * cols; + if (data.Length < need) + throw new ArgumentException( + string.Format("data length {0} < required prefix {1} for shape ({2}, {3})", + data.Length, need, rows, cols)); + return new Matrix(data, rows, cols, takeOwnership: true); + } + private Matrix(double[] data, int rows, int cols, bool takeOwnership) { _data = data; @@ -356,7 +374,11 @@ public void AddInPlace(Matrix rhs) { if (_rows != rhs._rows || _cols != rhs._cols) throw new ArgumentException("Matrices must have equal shape to add"); - for (int i = 0; i < _data.Length; i++) + // Bound the loop by active rows*cols, not _data.Length, so + // matrices wrapped via WrapPrefixNoClone (pool-backed, + // possibly oversized) don't touch the suffix. + int active = _rows * _cols; + for (int i = 0; i < active; i++) _data[i] += rhs._data[i]; } @@ -365,8 +387,11 @@ public void AddInPlace(Matrix rhs) /// public Matrix Divide(double divisor) { - var result = new double[_data.Length]; - for (int i = 0; i < _data.Length; i++) + // Allocate result sized exactly to rows*cols, not _data.Length + // (which may be larger if this matrix wraps a pooled buffer). + int active = _rows * _cols; + var result = new double[active]; + for (int i = 0; i < active; i++) result[i] = _data[i] / divisor; return new Matrix(result, _rows, _cols); } diff --git a/pwiz_tools/OspreySharp/OspreySharp.Scoring/SpectralScorer.cs b/pwiz_tools/OspreySharp/OspreySharp.Scoring/SpectralScorer.cs index effc80ded48..099a7de6383 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Scoring/SpectralScorer.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Scoring/SpectralScorer.cs @@ -110,17 +110,14 @@ public double[] PreprocessSpectrumForXcorr(Spectrum spectrum) /// /// Pool-aware preprocessing that writes the final sliding-window - /// result into the caller-supplied buffer. - /// All intermediate steps run in pure f32 to match Rust upstream - /// maccoss/osprey's native f32 XCorr path, avoiding the ~4e-6 drift - /// that f64-compute-then-narrow produces on HRAM data. HRAM main - /// search uses this variant so its per-window cache is f32 - /// (halving the 800 KB -> 400 KB per-spectrum cost vs f64). - /// - /// The parameter is currently unused; - /// a future change can replace it with a pooled f32 scratch set - /// if LOH pressure on the per-spectrum float[] allocations becomes - /// a hotspot. + /// result into the caller-supplied f32 + /// buffer (the per-spectrum HRAM cache). Internal binning, + /// windowing, and sliding-window math run in f64 on the pooled + /// double[] scratch fields; only the final cache store narrows to + /// f32. The cache stays f32 to preserve the 400 KB-per-spectrum + /// HRAM budget, while values now carry single-cast precision + /// rather than f32-cascade noise. Bit-equal cross-impl with Rust + /// preprocess_spectrum_for_xcorr_into. /// public void PreprocessSpectrumForXcorrInto( Spectrum spectrum, XcorrScratch scratch, float[] output) @@ -137,8 +134,28 @@ public void PreprocessSpectrumForXcorrInto( return; } - float[] pre = PreprocessSpectrumForXcorrF32(spectrum); - Array.Copy(pre, output, n); + // True scratch-pool path: f64 scratch (Binned / Windowed / + // Prefix double[]) for the windowing cascade; final sliding- + // window result is narrowed to f32 directly into `output`. + // No per-call allocation, no intermediate copy. + // + // Fallback: if no scratch was passed in, allocate as before. + // The fallback exists only for safety and code that hasn't + // been migrated; the HRAM per-window path (the hot one) + // always supplies scratch. + if (scratch != null && scratch.Binned.Length >= n) + { + // Binned accumulates via +=, so zero it per spectrum. + Array.Clear(scratch.Binned, 0, n); + PreprocessSpectrumForXcorrF32IntoBuffers( + spectrum, n, scratch.Binned, scratch.Windowed, + scratch.Prefix, output); + } + else + { + float[] pre = PreprocessSpectrumForXcorrF32(spectrum); + Array.Copy(pre, output, n); + } } /// @@ -478,77 +495,59 @@ private static void ApplySlidingWindowD(double[] spectrum, double[] prefix, doub } /// - /// Pure-f32 preprocess: bin + windowing normalization + sliding - /// window subtraction, all in f32. Mirrors Rust upstream - /// preprocess_spectrum_for_xcorr in - /// osprey-scoring/src/lib.rs. Used by the calibration XCorr path - /// to stay bit-equivalent with Rust's native f32 arithmetic (vs - /// the ~4e-6 drift that f64-compute-then-narrow produces). + /// Allocating wrapper around the f64-internal / f32-storage + /// preprocess pipeline. Bin + windowing normalization + sliding + /// window subtraction run in f64; final cache store narrows to + /// f32. Mirrors Rust preprocess_spectrum_for_xcorr on the + /// HRAM cache write boundary. Used by code paths that don't + /// supply a pooled . /// public float[] PreprocessSpectrumForXcorrF32(Spectrum spectrum) { int n = _binConfig.NBins; - float[] binned = new float[n]; - float[] windowed = new float[n]; - float[] prefix = new float[n + 1]; + double[] binned = new double[n]; + double[] windowed = new double[n]; + double[] prefix = new double[n + 1]; float[] preprocessed = new float[n]; - - for (int i = 0; i < spectrum.Mzs.Length; i++) - { - int bin = _binConfig.MzToBin(spectrum.Mzs[i]); - if (bin >= 0 && bin < n) - binned[bin] += (float)Math.Sqrt(spectrum.Intensities[i]); - } - ApplyWindowingNormalizationF(binned, windowed); - ApplySlidingWindowF(windowed, prefix, preprocessed); + PreprocessSpectrumForXcorrF32IntoBuffers( + spectrum, n, binned, windowed, prefix, preprocessed); return preprocessed; } - private static void ApplyWindowingNormalizationF(float[] spectrum, float[] result) + // Body of the f64-internal / f32-storage preprocessing pipeline, + // parameterized over pre-allocated work buffers. Caller is + // responsible for zeroing `binned` (this method accumulates into + // it via +=). The other buffers are fully overwritten. Bin sqrt + // is `Math.Sqrt(float)` (implicit float->double widening before + // sqrt) for bit-for-bit parity with Rust `(intensity as f64).sqrt()`. + private void PreprocessSpectrumForXcorrF32IntoBuffers( + Spectrum spectrum, int n, + double[] binned, double[] windowed, double[] prefix, float[] preprocessed) { - int n = spectrum.Length; - const int numWindows = 10; - int windowSize = (n / numWindows) + 1; - - float globalMax = 0.0f; - for (int i = 0; i < n; i++) - if (spectrum[i] > globalMax) - globalMax = spectrum[i]; - float threshold = globalMax * 0.05f; - - Array.Clear(result, 0, n); - - for (int w = 0; w < numWindows; w++) + for (int i = 0; i < spectrum.Mzs.Length; i++) { - int start = w * windowSize; - int end = Math.Min((w + 1) * windowSize, n); - if (start >= end) - break; - - float windowMax = 0.0f; - for (int i = start; i < end; i++) - if (spectrum[i] > windowMax) - windowMax = spectrum[i]; - - if (windowMax > 0.0f) - { - float normFactor = 50.0f / windowMax; - for (int i = start; i < end; i++) - { - if (spectrum[i] > threshold) - result[i] = spectrum[i] * normFactor; - } - } + int bin = _binConfig.MzToBin(spectrum.Mzs[i]); + if (bin >= 0 && bin < n) + binned[bin] += Math.Sqrt(spectrum.Intensities[i]); } + ApplyWindowingNormalizationD(binned, windowed); + ApplySlidingWindowDIntoF32(windowed, prefix, preprocessed); } - private static void ApplySlidingWindowF(float[] spectrum, float[] prefix, float[] result) + /// + /// Comet-style sliding-window subtraction in f64 with a final + /// narrowing store into an f32 cache buffer. Twin of + /// + /// where the result array is f32. Used by the f64-internal / + /// f32-storage cache build path. + /// + private static void ApplySlidingWindowDIntoF32(double[] spectrum, double[] prefix, float[] result) { int n = spectrum.Length; const int offset = XCORR_WINDOW_OFFSET; - float normFactor = 1.0f / (2 * offset); + double normFactor = 1.0 / (2 * offset); - prefix[0] = 0.0f; + prefix[0] = 0.0; for (int i = 0; i < n; i++) prefix[i + 1] = prefix[i] + spectrum[i]; @@ -556,9 +555,13 @@ private static void ApplySlidingWindowF(float[] spectrum, float[] prefix, float[ { int left = Math.Max(0, i - offset); int right = Math.Min(n, i + offset + 1); - float windowSum = prefix[right] - prefix[left]; - float sumExcludingCenter = windowSum - spectrum[i]; - result[i] = spectrum[i] - sumExcludingCenter * normFactor; + double windowSum = prefix[right] - prefix[left]; + double sumExcludingCenter = windowSum - spectrum[i]; + double centered = spectrum[i] - sumExcludingCenter * normFactor; + // f64 -> f32 narrowing at the final store: single + // deterministic rounding, identical bits cross-impl when + // the f64 inputs agree. + result[i] = (float)centered; } } diff --git a/pwiz_tools/OspreySharp/OspreySharp.Scoring/TukeyMedianPolish.cs b/pwiz_tools/OspreySharp/OspreySharp.Scoring/TukeyMedianPolish.cs index 248265ddd52..0d47a237cf4 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Scoring/TukeyMedianPolish.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Scoring/TukeyMedianPolish.cs @@ -536,25 +536,32 @@ private static double ComputeR2(List predicted, List observed) private static double PearsonCorrelationRaw(List x, List y) { + // Single-pass moment form: matches Rust pearson_correlation_raw in + // osprey-scoring/src/lib.rs and the existing PearsonCorrelation.Pearson + // helper. Aligning here brings median_polish_residual_correlation to + // cross-impl bit-equality. int n = Math.Min(x.Count, y.Count); if (n < 2) return 0.0; - double mx = 0.0, my = 0.0; - for (int i = 0; i < n; i++) { mx += x[i]; my += y[i]; } - mx /= n; my /= n; - double sxy = 0.0, sxx = 0.0, syy = 0.0; + double dn = n; + double sx = 0.0, sy = 0.0; + double sx2 = 0.0, sy2 = 0.0, sxy = 0.0; for (int i = 0; i < n; i++) { - double dx = x[i] - mx; - double dy = y[i] - my; - sxy += dx * dy; - sxx += dx * dx; - syy += dy * dy; + double xi = x[i]; + double yi = y[i]; + sx += xi; + sy += yi; + sx2 += xi * xi; + sy2 += yi * yi; + sxy += xi * yi; } - if (sxx < 1e-30 || syy < 1e-30) + + double denom = (dn * sx2 - sx * sx) * (dn * sy2 - sy * sy); + if (denom < 1e-30) return 0.0; - return sxy / Math.Sqrt(sxx * syy); + return (dn * sxy - sx * sy) / Math.Sqrt(denom); } } } diff --git a/pwiz_tools/OspreySharp/OspreySharp.Scoring/XcorrScratchPool.cs b/pwiz_tools/OspreySharp/OspreySharp.Scoring/XcorrScratchPool.cs index ed4ffed714c..f13c75e3832 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Scoring/XcorrScratchPool.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Scoring/XcorrScratchPool.cs @@ -31,10 +31,12 @@ namespace pwiz.OspreySharp.Scoring /// A single set of scratch buffers reused across XCorr calls to avoid /// per-call LOH allocation on HRAM (NBins ~100K, ~800 KB per array). /// - /// Each field is a plain Large Object once pooled; the pool hands out - /// and takes back sets of four arrays, never the arrays individually, - /// so the binned/windowed/prefix/preprocessed buffers stay co-located - /// across the full preprocess -> score pipeline. + /// All preprocessing math runs in f64 on these buffers; the per- + /// spectrum cache that feeds the dot product is f32 and is owned by + /// the caller (see SpectralScorer.PreprocessSpectrumForXcorrInto). + /// The f64 scratch is shared between the calibration path + /// (XcorrAtScan) and the HRAM main-search cache build path, so a + /// single rented scratch services both. /// public sealed class XcorrScratch { @@ -108,6 +110,8 @@ public void Return(XcorrScratch s) return; Array.Clear(s.Binned, 0, s.Binned.Length); Array.Clear(s.VisitedBins, 0, s.VisitedBins.Length); + // BinnedF is zeroed by callers on each preprocess call so it does + // not need post-window zeroing here (would just duplicate work). _scratchBag.Add(s); } diff --git a/pwiz_tools/OspreySharp/OspreySharp.Test/CalibrationTest.cs b/pwiz_tools/OspreySharp/OspreySharp.Test/CalibrationTest.cs index 51499d5f6f9..f41bea2fede 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Test/CalibrationTest.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Test/CalibrationTest.cs @@ -406,8 +406,14 @@ public void TestCalibrationFilename() [TestMethod] public void TestCalibrationFilenameForInput() { + // Use Path.Combine for the directory-prefixed case so the test + // runs on both Windows (\) and Linux (/). Path.GetFileNameWithoutExtension + // only recognizes its host OS's separator, so a hard-coded + // Windows literal would fail under Linux/WSL even though the + // production code handles either OS correctly given an + // OS-native input path. Assert.AreEqual("sample.calibration.json", - CalibrationIO.CalibrationFilenameForInput(@"C:\data\sample.mzML")); + CalibrationIO.CalibrationFilenameForInput(Path.Combine("data", "sample.mzML"))); Assert.AreEqual("test.dia.calibration.json", CalibrationIO.CalibrationFilenameForInput("test.dia.mzML")); Assert.AreEqual("experiment.calibration.json", diff --git a/pwiz_tools/OspreySharp/OspreySharp.Test/CwtCandidateCodecTest.cs b/pwiz_tools/OspreySharp/OspreySharp.Test/CwtCandidateCodecTest.cs index 223d2e6eacd..f245c303b37 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Test/CwtCandidateCodecTest.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Test/CwtCandidateCodecTest.cs @@ -357,7 +357,7 @@ public void TestCsScoringPopulatesCwtCandidates() public void TestLoadCwtCandidatesFromRustParquet() { string baseDir = System.Environment.GetEnvironmentVariable(@"OSPREY_TEST_BASE_DIR") - ?? @"D:\test\osprey-runs"; + ?? DefaultTestBaseDir(); string path = System.IO.Path.Combine(baseDir, @"astral", @"_stage6_planning", @"Astral", @"Ast-2024-12-05_HeLa_3mzDIA_6mIIT_400-900_49.scores.parquet"); @@ -408,8 +408,8 @@ public void TestLoadCwtCandidatesFromRustParquet() /// Resolve the Stellar test data directory via the /// OSPREY_TEST_BASE_DIR environment variable used by /// ai/scripts/OspreySharp/Dataset-Config.ps1, falling back - /// to D:\test\osprey-runs for the default Windows - /// developer setup. Tests that read parquets call + /// to for the default developer + /// setup. Tests that read parquets call /// Assert.Inconclusive when the resolved file is missing, /// keeping the suite portable to environments without the test /// dataset. @@ -417,10 +417,24 @@ public void TestLoadCwtCandidatesFromRustParquet() private static string StellarBaseDir() { string baseDir = System.Environment.GetEnvironmentVariable(@"OSPREY_TEST_BASE_DIR") - ?? @"D:\test\osprey-runs"; + ?? DefaultTestBaseDir(); return System.IO.Path.Combine(baseDir, @"stellar"); } + /// + /// Default test-data root when OSPREY_TEST_BASE_DIR is not + /// set. Maps the Windows D:\test\osprey-runs developer layout + /// to its WSL/drvfs equivalent /mnt/d/test/osprey-runs on + /// Linux, so the parquet-dependent tests run out of the box under + /// WSL without requiring the env var to be set manually. + /// + private static string DefaultTestBaseDir() + { + return System.IO.Path.DirectorySeparatorChar == '/' + ? @"/mnt/d/test/osprey-runs" + : @"D:\test\osprey-runs"; + } + private static void AssertBitEqual(double expected, double actual, string label) { long expBits = System.BitConverter.DoubleToInt64Bits(expected); diff --git a/pwiz_tools/OspreySharp/OspreySharp.Test/IOTest.cs b/pwiz_tools/OspreySharp/OspreySharp.Test/IOTest.cs index b43cd86fcd1..90ab715710e 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.Test/IOTest.cs +++ b/pwiz_tools/OspreySharp/OspreySharp.Test/IOTest.cs @@ -1859,6 +1859,7 @@ private static ReconciliationFile MakeSampleReconciliationFile() new ForcedIntegrationEntry { EntryId = 201, ExpectedRt = 18.5, HalfWidth = 0.05 }, }, FormatVersion = ReconciliationFile.CurrentFormatVersion, + FileStems = new List { "round_trip" }, GapFillTargets = new List { new GapFillEntry @@ -2145,6 +2146,7 @@ public void TestRescoreHydrationRoundTrip() var reconFile = new ReconciliationFile { FormatVersion = ReconciliationFile.CurrentFormatVersion, + FileStems = new List { stem }, SearchHash = "abc123", LibraryHash = "lib-h", UseCwtPeakActions = new List @@ -2280,6 +2282,7 @@ public void TestRescoreHydrationRejectsActionEntryIdNotInStubs() var reconFile = new ReconciliationFile { FormatVersion = ReconciliationFile.CurrentFormatVersion, + FileStems = new List { stem }, SearchHash = "x", LibraryHash = "y", UseCwtPeakActions = new List @@ -2323,23 +2326,38 @@ public void TestRescoreHydrationRejectsActionEntryIdNotInStubs() /// does between first-pass FDR and Stage 6 — see /// AnalysisPipeline's "First-pass compaction" block. /// + /// Compaction predicate is the UNION of (a) the local-FDR + /// predicate (peptide_q OR protein_q pass) AND (b) every + /// base_id that has a reconciliation action emitted by the + /// planner. The planner runs cross-file consensus rescue + /// (ConsensusRts.Compute), so an entry whose own file fails + /// local first-pass FDR can still be a reconciliation target + /// when its peptide passes FDR in a sibling file. Without the + /// union, the local-FDR-only predicate drops those entries and + /// their planner actions are silently dropped too -- the + /// reconciled .scores.parquet ends up with stale Stage 4 apex_rt + /// / bounds for the rescued entries and the blib output diverges + /// from the in-memory straight-through pipeline. + /// /// Test layout (single file, five entries): /// idx 0: target id=1, peptide_q=0.005 (PASS peptide) /// idx 1: decoy id=0x80000001, base=1 (retained via target's base_id) - /// idx 2: target id=2, peptide_q=0.5 (FAIL — non-passing) - /// idx 3: decoy id=0x80000002, base=2 (dropped because base=2 not in pass set) + /// idx 2: target id=2, peptide_q=0.5 (FAIL local; PLANNER ACTION at (f,2)) + /// idx 3: decoy id=0x80000002, base=2 (retained via target's base_id from action) /// idx 4: target id=3, peptide_q=0.5, protein_q=0.005 (PASS via protein-rescue) /// - /// With ProteinFdr=0.01, base_ids {1, 3} pass. After compaction: - /// idx 0: id=1, idx 1: id=0x80000001, idx 2: id=3. + /// With ProteinFdr=0.01 AND the planner-action union: + /// local-FDR pass: base_ids {1, 3} + /// planner action targets: base_ids {1, 2, 3} + /// UNION: base_ids {1, 2, 3} — all 5 entries survive compaction. /// /// Reconciliation actions are seeded at: - /// (f, 0) on id=1 → should remain at (f, 0) - /// (f, 4) on id=3 → should move to (f, 2) - /// (f, 2) on id=2 → should be dropped (entry compacted away) + /// (f, 0) on id=1 → remains at (f, 0) (no compaction shift) + /// (f, 4) on id=3 → remains at (f, 4) + /// (f, 2) on id=2 → remains at (f, 2) (entry survives via union) /// [TestMethod] - public void TestRescoreCompactionRekeysActionsAndDropsNonpassing() + public void TestRescoreCompactionUnionsActionsWithLocalFdrPredicate() { const string fileName = "f1"; var perFile = new List>> @@ -2375,30 +2393,30 @@ public void TestRescoreCompactionRekeysActionsAndDropsNonpassing() ProteinFdr = 0.01, }); - // Compaction stats. + // Compaction stats: planner-action union keeps base 2 alive. Assert.AreEqual(5, stats.EntriesBefore); - Assert.AreEqual(3, stats.EntriesAfter); - Assert.AreEqual(2, stats.FirstPassBaseIds); // base_ids {1, 3} - Assert.AreEqual(1, stats.DroppedActions); // the (f, 2) action + Assert.AreEqual(5, stats.EntriesAfter); + Assert.AreEqual(3, stats.FirstPassBaseIds); // base_ids {1, 2, 3} + Assert.AreEqual(0, stats.DroppedActions); // no actions dropped - // Per-file list compacted. + // Per-file list unchanged (all base_ids survived). Assert.AreEqual(1, inputs.PerFileEntries.Count); var got = inputs.PerFileEntries[0].Value; - Assert.AreEqual(3, got.Count); + Assert.AreEqual(5, got.Count); Assert.AreEqual(1u, got[0].EntryId); Assert.AreEqual(0x80000001u, got[1].EntryId); - Assert.AreEqual(3u, got[2].EntryId); + Assert.AreEqual(2u, got[2].EntryId); + Assert.AreEqual(0x80000002u, got[3].EntryId); + Assert.AreEqual(3u, got[4].EntryId); - // Reconciliation actions re-keyed. - Assert.AreEqual(2, inputs.ReconciliationActions.Count); + // All 3 reconciliation actions preserved. + Assert.AreEqual(3, inputs.ReconciliationActions.Count); Assert.IsInstanceOfType( inputs.ReconciliationActions[(fileName, 0)], typeof(ReconcileAction.UseCwtPeak)); Assert.IsInstanceOfType( - inputs.ReconciliationActions[(fileName, 2)], - typeof(ReconcileAction.ForcedIntegration)); - // The action at the dropped entry is gone, NOT silently re-keyed - // to a different surviving entry. - Assert.IsFalse(inputs.ReconciliationActions.ContainsKey((fileName, 1))); + inputs.ReconciliationActions[(fileName, 2)], typeof(ReconcileAction.UseCwtPeak)); + Assert.IsInstanceOfType( + inputs.ReconciliationActions[(fileName, 4)], typeof(ReconcileAction.ForcedIntegration)); } /// diff --git a/pwiz_tools/OspreySharp/OspreySharp.sln b/pwiz_tools/OspreySharp/OspreySharp.sln index 77fc194f6e0..da1d36be9b7 100644 --- a/pwiz_tools/OspreySharp/OspreySharp.sln +++ b/pwiz_tools/OspreySharp/OspreySharp.sln @@ -1,4 +1,4 @@ - + Microsoft Visual Studio Solution File, Format Version 12.00 # Visual Studio Version 17 VisualStudioVersion = 17.0.31903.59 diff --git a/pwiz_tools/OspreySharp/OspreySharp/AnalysisPipeline.cs b/pwiz_tools/OspreySharp/OspreySharp/AnalysisPipeline.cs index 8f99b7d40ee..5713fca70bf 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/AnalysisPipeline.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/AnalysisPipeline.cs @@ -212,6 +212,23 @@ private static bool RunTask(OspreyTask task, PipelineContext ctx) ctx.LogInfo(string.Format(@"[task] {0}: done ({1:F1}s)", task.Name, sw.Elapsed.TotalSeconds)); + // [STAGE-WALL] one line per task->stage with parseable format + // for Measure-Pipeline.ps1 / Osprey-workflow.html perf tables. + // MergeNodeTask emits its own stage7 + blib lines internally + // (one task -> two pipeline stages). + string stageName = task.Name switch + { + "PerFileScoring" => "stage1to4", + "FirstJoin" => "stage5", + "PerFileRescore" => "stage6", + _ => null, + }; + if (stageName != null) + { + ctx.LogInfo(string.Format(@"[STAGE-WALL] {0}: {1:F1}s", + stageName, sw.Elapsed.TotalSeconds)); + } + // Write sidecars whenever the task ran without setting a // non-zero exit code. Several tasks intentionally return // false on success to stop the pipeline at a configured diff --git a/pwiz_tools/OspreySharp/OspreySharp/OspreyDiagnostics.cs b/pwiz_tools/OspreySharp/OspreySharp/OspreyDiagnostics.cs index 9135557107a..8818cfd518b 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/OspreyDiagnostics.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/OspreyDiagnostics.cs @@ -102,6 +102,18 @@ public static class OspreyDiagnostics /// OSPREY_CAL_MATCH_ONLY: exit after cs_cal_match.txt dump. public static readonly bool CalMatchOnly = IsOne(@"OSPREY_CAL_MATCH_ONLY"); + /// + /// OSPREY_DUMP_MS2_CAL_ERRORS: dump the per-fragment MS2 mass errors + /// that feed the MS2 mass calibration mean/sd computation, to + /// cs_ms2_cal_errors.txt. Used to bisect a 1-ULP MS2 calibration + /// mean divergence cross-impl by checking whether the input + /// fragment errors are identical in value AND order across impls. + /// + public static readonly bool DumpMs2CalErrors = IsOne(@"OSPREY_DUMP_MS2_CAL_ERRORS"); + + /// OSPREY_MS2_CAL_ERRORS_ONLY: exit after cs_ms2_cal_errors.txt dump. + public static readonly bool Ms2CalErrorsOnly = IsOne(@"OSPREY_MS2_CAL_ERRORS_ONLY"); + /// /// OSPREY_DUMP_LDA_SCORES: dump per-entry LDA discriminant + q-value /// after LDA scoring (cs_lda_scores.txt). @@ -293,6 +305,15 @@ public static class OspreyDiagnostics /// OSPREY_STAGE7_PROTEIN_FDR_ONLY: exit after cs_stage7_protein_fdr.tsv dump. public static readonly bool Stage7ProteinFdrOnly = IsOne(@"OSPREY_STAGE7_PROTEIN_FDR_ONLY"); + /// + /// OSPREY_DUMP_DETECTED_PEPTIDES: dump the sorted set of detected + /// target peptides handed to BuildProteinParsimony for Stage 7 to + /// cs_stage7_detected_peptides.txt. Mirrors Rust + /// dump_stage7_detected_peptides; lets the protein-FDR input set + /// be diffed cross-impl before debugging downstream divergence. + /// + public static readonly bool DumpDetectedPeptides = IsOne(@"OSPREY_DUMP_DETECTED_PEPTIDES"); + /// /// OSPREY_DUMP_LOESS_FIT: dump the per-point LOESS fit state of the /// Stage 6 refit RTCalibration to cs_stage6_loess_fit.tsv. Used to @@ -418,14 +439,14 @@ public static void WriteCalScalarsAndGridDump( w.WriteLine(@"n_targets" + "\t" + targets.Count); w.WriteLine(@"n_decoys" + "\t" + decoys.Count); w.WriteLine(@"bins_per_axis" + "\t" + binsPerAxis); - w.WriteLine(@"rt_min" + "\t" + rtMin.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"rt_max" + "\t" + rtMax.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"mz_min" + "\t" + mzMin.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"mz_max" + "\t" + mzMax.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"rt_range" + "\t" + rtRange.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"mz_range" + "\t" + mzRange.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"rt_bin_width" + "\t" + rtBinWidth.ToString(@"F17", CultureInfo.InvariantCulture)); - w.WriteLine(@"mz_bin_width" + "\t" + mzBinWidth.ToString(@"F17", CultureInfo.InvariantCulture)); + w.WriteLine(@"rt_min" + "\t" + rtMin.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"rt_max" + "\t" + rtMax.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"mz_min" + "\t" + mzMin.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"mz_max" + "\t" + mzMax.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"rt_range" + "\t" + rtRange.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"mz_range" + "\t" + mzRange.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"rt_bin_width" + "\t" + rtBinWidth.ToString(@"G17", CultureInfo.InvariantCulture)); + w.WriteLine(@"mz_bin_width" + "\t" + mzBinWidth.ToString(@"G17", CultureInfo.InvariantCulture)); w.WriteLine(@"n_occupied" + "\t" + nOccupied); w.WriteLine(@"per_cell" + "\t" + perCell); w.WriteLine(@"seed" + "\t" + seed); @@ -552,8 +573,18 @@ public static void WriteCalMatchDump(int passNumber, double snr; if (!snrByEntryId.TryGetValue(entry.Id, out snr)) snr = 0.0; + // G17 (17 significant digits) for round-trip-safe f64. + // .NET Framework 4.7.2's F17 truncates output at ~15 + // significant digits and pads with zeros, so a + // string -> parse round-trip yields a different f64 + // than the original. G17 prints enough digits for the + // result to round-trip exactly back to the same f64. + // The cross-impl comparator parses both sides as + // numbers, so variable-width G17 output on C# vs + // fixed {:.17} on Rust is fine - both round-trip to + // the same f64 when the underlying value matches. w.WriteLine(string.Format(inv, - "{0}\t{1}\t{2}\t1\t{3}\t{4:F10}\t{5:F10}\t{6:F10}\t{7}\t{8:F10}\t{9:F10}", + "{0}\t{1}\t{2}\t1\t{3}\t{4:G17}\t{5:G17}\t{6:G17}\t{7}\t{8:G17}\t{9:G17}", entry.Id, entry.IsDecoy ? 1 : 0, entry.Charge, @@ -587,6 +618,57 @@ public static void WriteCalMatchDump(int passNumber, /// Dump per-entry LDA discriminant + q-value (cs_lda_scores.txt), /// sorted by entry_id, F10-formatted. /// + /// + /// Dump per-fragment MS2 mass errors that feed the MS2 calibration + /// (cs_ms2_cal_errors.txt). One row per + /// (entry_id, fragment_order_index, error) triple, sorted by + /// (entry_id, fragment_order_index). Mirrors Rust + /// dump_ms2_cal_errors. Use to bisect 1-ULP MS2 calibration mean + /// divergence cross-impl: if dumps match exactly, divergence is in + /// the mean/sd computation (e.g. naive sum vs Welford); if not, + /// divergence is in per-match top-N fragment selection or per- + /// fragment error arithmetic. + /// + /// Pass only the matches contributing to calibration (i.e. post- + /// LDA passing targets with q_value <= calibration_fdr and + /// snr >= MIN_SNR_FOR_RT_CAL) so this dump aligns with what + /// MzCalibration.CalculateSingleCalibration actually sees. + /// + public static void WriteMs2CalErrorsDump(IEnumerable contributingMatches) + { + var sortedByEntry = contributingMatches.OrderBy(m => m.EntryId).ToArray(); + var inv = CultureInfo.InvariantCulture; + int nRows = 0; + int nMatches = 0; + using (var w = new StreamWriter(@"cs_ms2_cal_errors.txt")) + { + w.WriteLine("entry_id\tfrag_order_idx\terror"); + foreach (var m in sortedByEntry) + { + nMatches++; + if (m.Ms2MassErrors == null) continue; + for (int i = 0; i < m.Ms2MassErrors.Length; i++) + { + double err = m.Ms2MassErrors[i]; + // G17: round-trip-safe canonical form. C# G17 and + // Rust ryu (the `{}` default) format the same f64 + // value with different strings (e.g. "1.23E-05" + // vs "0.0000123"), so this dump is intended to be + // compared numerically via a Python diff script, + // not via SHA byte-equality. Inside one impl, the + // G17 strings round-trip the f64 bits exactly. + w.WriteLine(string.Format(inv, + "{0}\t{1}\t{2:G17}", + m.EntryId, i, err)); + nRows++; + } + } + } + LogAction(string.Format(inv, + @"[COUNT] Wrote MS2 cal errors dump: cs_ms2_cal_errors.txt ({0} rows across {1} matches)", + nRows, nMatches)); + } + public static void WriteLdaScoresDump(int passNumber, IEnumerable matchArray) { var sortedByEntry = matchArray.OrderBy(m => m.EntryId).ToArray(); @@ -596,8 +678,9 @@ public static void WriteLdaScoresDump(int passNumber, IEnumerable + /// Write the sorted set of detected target peptides handed to + /// BuildProteinParsimony for Stage 7 to cs_stage7_detected_peptides.txt. + /// Mirrors Rust dump_stage7_detected_peptides. Gated by + /// . + /// + public static void WriteStage7DetectedPeptidesDump(HashSet detectedPeptides) + { + const string path = @"cs_stage7_detected_peptides.txt"; + var sorted = new List(detectedPeptides); + sorted.Sort(StringComparer.Ordinal); + File.WriteAllLines(path, sorted); + LogAction(string.Format(@"[DIAG] Wrote {0} ({1} entries)", path, sorted.Count)); + } + } } diff --git a/pwiz_tools/OspreySharp/OspreySharp/Program.cs b/pwiz_tools/OspreySharp/OspreySharp/Program.cs index fde179c9a01..8b141e896fc 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/Program.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/Program.cs @@ -38,7 +38,18 @@ static class Program // Tracks the Rust Osprey upstream version this OspreySharp port // is aligned with. Used in parquet footer metadata; the Phase 3 // validator requires same major.minor across cross-impl handoff. - internal const string VERSION = "26.6.0"; + // TODO: 26.6.1 bumped the version string but the algorithmic + // payload of v26.6.1 (reconciliation pairing library-supplied + // decoys by base_id instead of stripping a DECOY_ prefix in + // compute_consensus_rts + plan_reconciliation) is NOT yet + // ported to this side. It does not affect reverse-decoy mode + // (Stellar, DecoysInLibrary=false), but it WILL affect any + // dataset run with --decoys-in-library. See osprey + // release-notes/RELEASE_NOTES_v26.6.1.md and the + // test_consensus_rts_pairs_library_decoy_by_base_id + + // test_plan_reconciliation_includes_library_decoy_via_base_id + // regression tests on the Rust side. + internal const string VERSION = "26.6.1"; internal const string VERSION_STRING = VERSION; static int Main(string[] args) diff --git a/pwiz_tools/OspreySharp/OspreySharp/RescoreCompaction.cs b/pwiz_tools/OspreySharp/OspreySharp/RescoreCompaction.cs index d7cc545ebe5..793a09ac04a 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/RescoreCompaction.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/RescoreCompaction.cs @@ -147,6 +147,26 @@ public static Stats Apply(RescoreInputs inputs, OspreyConfig config) // vec_idx values. We have to do this BEFORE per_file_entries // shrinks, because the saved (file, vec_idx) keys are about // to become stale. + // + // ALSO: union firstPassBaseIds with the base_ids of every + // entry the planner emits a reconciliation action for. The + // planner runs compute_consensus_rts (cross-file consensus + // rescue), so an entry whose own file fails local first-pass + // FDR can still be a reconciliation target when its peptide + // passes FDR in a sibling file. The local-FDR-only predicate + // above drops those entries, then this loop silently drops + // their planner actions (counted as "DroppedActions"), and + // the rescore engine never applies them -- the reconciled + // .scores.parquet ends up with stale Stage 4 apex_rt / + // bounds for ~200 rows per Stellar file (0.04%) and the + // blib output diverges from the in-memory straight-through + // pipeline. Mirrors the Rust Option B fix in + // rescore.rs::run_rescore: first_pass_base_ids = UNION of + // local-FDR predicate AND entry_ids in + // reconciliation_actions_pre. Bisected via the C# in-memory + // vs C# HPC-chain strict-rehydration test (Stage 6 boundary + // revealed apex_rt / bounds drift on cross-file-rescued + // entries that were dropped by worker compaction). var actionsById = new Dictionary<(string FileName, uint EntryId), ReconcileAction>( inputs.ReconciliationActions.Count); // Build a lookup map so the per-action entry_id resolution is @@ -163,7 +183,13 @@ public static Stats Apply(RescoreInputs inputs, OspreyConfig config) continue; if (vecIdx < 0 || vecIdx >= entries.Count) continue; - actionsById[(fileName, entries[vecIdx].EntryId)] = kvp.Value; + uint entryId = entries[vecIdx].EntryId; + actionsById[(fileName, entryId)] = kvp.Value; + // Extend firstPassBaseIds so the entry survives the local-FDR + // compaction predicate above. Adding the masked base_id + // (decoy bit stripped) keeps both the target and its paired + // decoy alive, preserving the target-decoy invariant. + firstPassBaseIds.Add(entryId & BASE_ID_MASK); } // 3. Compact each per-file entry list in place. diff --git a/pwiz_tools/OspreySharp/OspreySharp/RescoreHydration.cs b/pwiz_tools/OspreySharp/OspreySharp/RescoreHydration.cs index 296a6351805..0dcf9f4a634 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/RescoreHydration.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/RescoreHydration.cs @@ -85,6 +85,19 @@ public class RescoreInputs /// public Dictionary> PerFileConsensusTargets { get; set; } + /// + /// Full set of file stems participating in the planner's join, as + /// read from reconciliation.json's file_stems field + /// (v2+). Per-file Stage 6 rescore workers carry this through so + /// they can compute the join-wide reconciliation parameter hash — + /// the worker's OspreyConfig.InputFiles only has its single + /// parquet, but the hash that the downstream + /// --join-at-pass=2 merge node validates is computed over + /// all files. Empty list when reading a v1 envelope (the worker + /// falls back to its InputFiles stems in that case). + /// + public List JoinFileStems { get; set; } + /// Total non-Keep reconciliation actions across all files. public int TotalActions => ReconciliationActions.Count; @@ -215,6 +228,11 @@ public static RescoreInputs HydrateReconciliationOverlay( var refinedCalibrations = new Dictionary(); var perFileGapFill = new Dictionary>(); var reconciliationActions = new Dictionary<(string, int), ReconcileAction>(); + // Captured from the first envelope's file_stems field (v2+); + // every subsequent envelope must carry the same set so the + // worker's join-wide reconciliation hash matches the planner's. + // Mirrors the consistency check in Rust hydrate_for_rescore. + List joinFileStems = null; for (int i = 0; i < perFileEntries.Count; i++) { @@ -249,6 +267,32 @@ public static RescoreInputs HydrateReconciliationOverlay( reconPath, ex.Message), ex); } + // Capture / validate file_stems across all envelopes. + // ReconciliationFile.Load already rejects any envelope whose + // format_version != CurrentFormatVersion (currently 2), so by + // the time we reach here `envelope.FileStems` must be the + // planner's full join file set -- a non-empty list, identical + // across every envelope produced by a single planner step. + // Any disagreement (including unexpected empty stems) means + // the on-disk envelopes were produced by different planner + // steps and indicates a corrupted hand-off. Mirrors the + // consistency check in Rust's hydrate_for_rescore. + var envelopeStems = NormalizeStems(envelope.FileStems); + if (joinFileStems == null) + { + joinFileStems = envelopeStems; + } + else if (!StemsEqual(joinFileStems, envelopeStems)) + { + throw new InvalidDataException(string.Format( + "HydrateReconciliationOverlay: reconciliation.json {0} carries a " + + "different file_stems set than its siblings (planner inconsistency). " + + "Expected: [{1}]; got: [{2}]", + reconPath, + string.Join(", ", joinFileStems), + string.Join(", ", envelopeStems))); + } + // Build entry_id -> vec_idx map from the loaded stubs so the // planner's entry_id-keyed actions can be rehomed onto // (file_name, vec_idx) keys the rescore engine consumes. @@ -322,9 +366,50 @@ public static RescoreInputs HydrateReconciliationOverlay( RefinedCalibrations = refinedCalibrations, PerFileGapFill = perFileGapFill, PerFileConsensusTargets = null, + JoinFileStems = joinFileStems ?? new List(), }; } + /// + /// Sort + dedup a list of file stems (Ordinal). Returns a new list; + /// the input is not mutated. Empty / null input becomes an empty + /// list. Used to canonicalize the file_stems field from + /// each reconciliation.json envelope before consistency + /// checks across siblings. + /// + private static List NormalizeStems(IList stems) + { + if (stems == null || stems.Count == 0) + return new List(); + var result = new List(stems.Count); + foreach (var s in stems) + { + if (!string.IsNullOrEmpty(s)) + result.Add(s); + } + result.Sort(StringComparer.Ordinal); + for (int i = result.Count - 1; i > 0; i--) + { + if (string.Equals(result[i], result[i - 1], StringComparison.Ordinal)) + result.RemoveAt(i); + } + return result; + } + + /// + /// Ordinal element-wise equality on two pre-normalized stem lists. + /// + private static bool StemsEqual(IList a, IList b) + { + if (a.Count != b.Count) return false; + for (int i = 0; i < a.Count; i++) + { + if (!string.Equals(a[i], b[i], StringComparison.Ordinal)) + return false; + } + return true; + } + /// /// Inverse of scores_path_for_input: given /// /data/sample1.scores.parquet, produce a synthetic input diff --git a/pwiz_tools/OspreySharp/OspreySharp/RescoreWorker.cs b/pwiz_tools/OspreySharp/OspreySharp/RescoreWorker.cs index df97636d13d..5a14057171a 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/RescoreWorker.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/RescoreWorker.cs @@ -22,7 +22,6 @@ */ using pwiz.OspreySharp.Core; -using pwiz.OspreySharp.Tasks; namespace pwiz.OspreySharp { diff --git a/pwiz_tools/OspreySharp/OspreySharp/Tasks/AbstractScoringTask.cs b/pwiz_tools/OspreySharp/OspreySharp/Tasks/AbstractScoringTask.cs index 079e5f2ef01..cae12a8adba 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/Tasks/AbstractScoringTask.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/Tasks/AbstractScoringTask.cs @@ -2975,6 +2975,18 @@ protected List DeduplicatePairs(List entries) deduped.Add(pair.Value); } + // Sort by EntryId for deterministic order regardless of Dictionary + // iteration. Mirrors Rust deduplicate_pairs's final sort_by_key + // (pipeline.rs:6123) and its comment: "Without this, the random + // HashMap order propagates to SVM feature matrix row ordering, + // causing non-deterministic gradient updates and model weights." + // Cross-impl: the straight-through path feeds these entries + // directly to Percolator (no parquet round-trip to mask the + // un-sorted order), so an unsorted dedup output cascades into + // SVM working-set divergence and ~190-precursor / ~270-peptide + // first-pass FDR drift on Stellar Single. + deduped.Sort((a, b) => a.EntryId.CompareTo(b.EntryId)); + int removed = entries.Count - deduped.Count; if (removed > 0) { diff --git a/pwiz_tools/OspreySharp/OspreySharp/Tasks/FirstJoinTask.cs b/pwiz_tools/OspreySharp/OspreySharp/Tasks/FirstJoinTask.cs index f815672f253..5784193359c 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/Tasks/FirstJoinTask.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/Tasks/FirstJoinTask.cs @@ -270,11 +270,12 @@ public override bool Run(PipelineContext ctx) // before compaction drops any rows, so the cross-impl diff // sees both targets and decoys. if (OspreyDiagnostics.DumpPercolator) - { OspreyDiagnostics.WriteStage5PercolatorDump(perFileEntries); - if (OspreyDiagnostics.PercolatorOnly) - OspreyDiagnostics.ExitAfterDump(@"OSPREY_PERCOLATOR_ONLY"); - } + // 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. + if (OspreyDiagnostics.PercolatorOnly) + OspreyDiagnostics.ExitAfterDump(@"OSPREY_PERCOLATOR_ONLY"); // First-pass protein FDR: runs on the full pre-compaction // peptide pool so target and decoy proteins compete on a @@ -992,6 +993,24 @@ private int WriteReconciliationFiles( gapFillByFileOut[kvp.Key] = copy; } + // The multi-file stems set goes into every per-file + // reconciliation.json so a worker rescoring its single + // parquet can compute the join-wide reconciliation hash that + // --join-at-pass=2 will validate against. Sort + dedup once; + // BuildReconciliationFile copies the list into the wire form. + var joinFileStems = new List(perFileEntries.Count); + foreach (var fEntry in perFileEntries) + { + if (!string.IsNullOrEmpty(fEntry.Key)) + joinFileStems.Add(fEntry.Key); + } + joinFileStems.Sort(StringComparer.Ordinal); + for (int i = joinFileStems.Count - 1; i > 0; i--) + { + if (string.Equals(joinFileStems[i], joinFileStems[i - 1], StringComparison.Ordinal)) + joinFileStems.RemoveAt(i); + } + int failures = 0; foreach (var kvp in perFileEntries) { @@ -1014,7 +1033,7 @@ private int WriteReconciliationFiles( var reconFile = BuildReconciliationFile( fileEntries, fileActions, fileGapFill, refinedCalibrations.TryGetValue(fileName, out var fileCal) ? fileCal : null, - searchHash, libraryHash); + searchHash, libraryHash, joinFileStems); try { ReconciliationFile.Save(reconPath, reconFile); @@ -1093,7 +1112,8 @@ private static ReconciliationFile BuildReconciliationFile( IReadOnlyList gapFillTargets, RTCalibration refinedCalibration, string searchHash, - string libraryHash) + string libraryHash, + IReadOnlyList joinFileStems) { var useCwt = new List(); var forced = new List(); @@ -1164,8 +1184,15 @@ private static ReconciliationFile BuildReconciliationFile( } } + // Defensive copy so a later caller-side mutation of + // joinFileStems doesn't leak into the serialized envelope. + var fileStems = joinFileStems != null + ? new List(joinFileStems) + : new List(); + return new ReconciliationFile { + FileStems = fileStems, ForcedIntegrationActions = forced, FormatVersion = ReconciliationFile.CurrentFormatVersion, GapFillTargets = gap, @@ -1187,7 +1214,7 @@ private void RunFdr( switch (config.FdrMethod) { case FdrMethod.Percolator: - RunPercolatorFdr(perFileEntries, fullLibrary, config); + RunPercolatorFdr(perFileEntries, fullLibrary, config, _ctx); break; case FdrMethod.Simple: @@ -1206,12 +1233,42 @@ private void RunFdr( /// /// Run Percolator-based FDR control. /// Builds PercolatorEntry objects from FdrEntry stubs and runs Percolator. + /// Static + internal so can call it for + /// the 2nd-pass run after Stage 6 reconciliation (the HPC distribution + /// case where workers wrote reconciled .scores.parquet but no + /// .2nd-pass.fdr_scores.bin sidecars; mirrors Rust pipeline.rs:4394-4468). /// - private void RunPercolatorFdr( + internal static void RunPercolatorFdr( List>> perFileEntries, List fullLibrary, - OspreyConfig config) + OspreyConfig config, + PipelineContext ctx, + string passLabel = "First-pass") { + // Sort each file's entries by EntryId 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 + // DeduplicatePairs (AbstractScoringTask.cs), 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. Mirrors Rust pipeline.rs::run_percolator_fdr. + foreach (var kvp in perFileEntries) + { + kvp.Value.Sort((a, b) => + { + int c = a.EntryId.CompareTo(b.EntryId); + if (c != 0) return c; + c = a.Charge.CompareTo(b.Charge); + if (c != 0) return c; + c = a.ScanNumber.CompareTo(b.ScanNumber); + if (c != 0) return c; + return a.ParquetIndex.CompareTo(b.ParquetIndex); + }); + } + // Build PercolatorEntry list from all files var percEntries = new List(); @@ -1250,9 +1307,23 @@ private void RunPercolatorFdr( nInputDecoys++; else nInputTargets++; + // PSM Id must uniquely identify each observation so the + // result -> FdrEntry write-back can score every row + // independently. EntryId alone is NOT unique within a + // file: a single base_id with multiple scan-time + // observations (different scan numbers, same charge, + // same modified_sequence) shares one EntryId. Using + // "{fileName}_{EntryId}" collided those rows in + // resultMap, leaving the last-inserted score + // overwriting every same-EntryId observation's + // FdrEntry.Score and producing 176-185 score + // divergences per file vs. Rust's 4-component psm_id. + // Mirrors osprey-fdr/src/percolator.rs:5978-5980. percEntries.Add(new PercolatorEntry { - Id = string.Format("{0}_{1}", fileName, fdrEntry.EntryId), + Id = string.Format("{0}_{1}_{2}_{3}", + fileName, fdrEntry.ModifiedSequence, + fdrEntry.Charge, fdrEntry.ScanNumber), FileName = fileName, Peptide = fdrEntry.ModifiedSequence, Charge = fdrEntry.Charge, @@ -1263,14 +1334,15 @@ private void RunPercolatorFdr( } } - _ctx.LogInfo(string.Format( - "[COUNT] Percolator input: {0} entries ({1} targets, {2} decoys, {3} features)", - percEntries.Count, nInputTargets, nInputDecoys, NUM_PIN_FEATURES)); - _ctx.LogInfo(string.Format( - "[COUNT] Percolator features computed: {0} entries with PIN features, {1} fallback", - nWithFeatures, nWithoutFeatures)); + ctx.LogInfo(string.Format( + "[COUNT] {0} Percolator input: {1} entries ({2} targets, {3} decoys, {4} features)", + passLabel, percEntries.Count, nInputTargets, nInputDecoys, NUM_PIN_FEATURES)); + ctx.LogInfo(string.Format( + "[COUNT] {0} Percolator features computed: {1} entries with PIN features, {2} fallback", + passLabel, nWithFeatures, nWithoutFeatures)); - _ctx.LogInfo(string.Format("Running Percolator on {0} entries...", percEntries.Count)); + ctx.LogInfo(string.Format("Running {0} Percolator on {1} entries...", + passLabel, percEntries.Count)); var percConfig = new PercolatorConfig { @@ -1294,7 +1366,7 @@ private void RunPercolatorFdr( if (percConfig.MaxTrainSize > 0 && percEntries.Count > percConfig.MaxTrainSize * 2) { - results = RunPercolatorStreaming(percEntries, percConfig); + results = RunPercolatorStreaming(percEntries, percConfig, ctx, passLabel); } else { @@ -1311,7 +1383,12 @@ private void RunPercolatorFdr( string fileName = kvp.Key; foreach (var fdrEntry in kvp.Value) { - string id = string.Format("{0}_{1}", fileName, fdrEntry.EntryId); + // 4-component psm_id matches the construction in + // the loop above so each FdrEntry pulls back its + // own PercolatorResult. Mirrors Rust direct path. + string id = string.Format("{0}_{1}_{2}_{3}", + fileName, fdrEntry.ModifiedSequence, + fdrEntry.Charge, fdrEntry.ScanNumber); PercolatorResult result; if (resultMap.TryGetValue(id, out result)) { @@ -1324,7 +1401,6 @@ private void RunPercolatorFdr( } } } - // Log FDR results int nTargetPassing = 0; int nDecoyPassing = 0; @@ -1342,19 +1418,19 @@ private void RunPercolatorFdr( fileTargets++; } } - _ctx.LogInfo(string.Format( - "[COUNT] Percolator pass [{0}]: {1} targets, {2} decoys at {3:P0} FDR", - kvp.Key, fileTargets, fileDecoys, config.RunFdr)); + ctx.LogInfo(string.Format( + "[COUNT] {0} Percolator pass [{1}]: {2} targets, {3} decoys at {4:P0} FDR", + passLabel, kvp.Key, fileTargets, fileDecoys, config.RunFdr)); nTargetPassing += fileTargets; nDecoyPassing += fileDecoys; } - _ctx.LogInfo(string.Format( - "Percolator results: {0} targets, {1} decoys pass {2:P1} FDR", - nTargetPassing, nDecoyPassing, config.RunFdr)); - _ctx.LogInfo(string.Format( - "[COUNT] First-pass total across files: {0}", - nTargetPassing)); + ctx.LogInfo(string.Format( + "{0} Percolator results: {1} targets, {2} decoys pass {3:P1} FDR", + passLabel, nTargetPassing, nDecoyPassing, config.RunFdr)); + ctx.LogInfo(string.Format( + "[COUNT] {0} total across files: {1}", + passLabel, nTargetPassing)); // Compute unique precursors across files (best q-value per modseq+charge) var bestQByPrecursor = new Dictionary(StringComparer.Ordinal); @@ -1373,9 +1449,9 @@ private void RunPercolatorFdr( bestQByPrecursor[pkey] = q; } } - _ctx.LogInfo(string.Format( - "[COUNT] First-pass unique precursors (best q across files): {0}", - bestQByPrecursor.Count)); + ctx.LogInfo(string.Format( + "[COUNT] {0} unique precursors (best q across files): {1}", + passLabel, bestQByPrecursor.Count)); } /// @@ -1385,7 +1461,7 @@ private void RunPercolatorFdr( /// 21-feature vector is computed during coelution scoring in /// and stored on the entry. /// - private double[] BuildBasicFeatures( + private static double[] BuildBasicFeatures( FdrEntry entry, Dictionary libraryById) { double[] features = new double[NUM_PIN_FEATURES]; @@ -1460,9 +1536,11 @@ private double[] BuildBasicFeatures( /// same helpers the direct path calls internally, so both paths /// select identical 300K subsets when given identical input. /// - private PercolatorResults RunPercolatorStreaming( + private static PercolatorResults RunPercolatorStreaming( List percEntries, - PercolatorConfig percConfig) + PercolatorConfig percConfig, + PipelineContext ctx, + string passLabel) { int n = percEntries.Count; int maxTrain = percConfig.MaxTrainSize; @@ -1487,9 +1565,9 @@ private PercolatorResults RunPercolatorStreaming( if (labels[bestIdx[i]]) dedupDecoys++; else dedupTargets++; } - _ctx.LogInfo(string.Format( - "[COUNT] Percolator streaming best-per-precursor: {0} entries ({1} targets, {2} decoys) from {3} total", - bestIdx.Length, dedupTargets, dedupDecoys, n)); + ctx.LogInfo(string.Format( + "[COUNT] {0} Percolator streaming best-per-precursor: {1} entries ({2} targets, {3} decoys) from {4} total", + passLabel, bestIdx.Length, dedupTargets, dedupDecoys, n)); // 2. Peptide-grouped subsample if dedup count still exceeds MaxTrainSize. int[] trainSubsetGlobalIdx; @@ -1522,9 +1600,9 @@ private PercolatorResults RunPercolatorStreaming( if (labels[trainSubsetGlobalIdx[i]]) subDecoys++; else subTargets++; } - _ctx.LogInfo(string.Format( - "[COUNT] Percolator streaming subsample: {0} entries ({1} targets, {2} decoys)", - trainSubsetGlobalIdx.Length, subTargets, subDecoys)); + ctx.LogInfo(string.Format( + "[COUNT] {0} Percolator streaming subsample: {1} entries ({2} targets, {3} decoys)", + passLabel, trainSubsetGlobalIdx.Length, subTargets, subDecoys)); // 3. Build subset entry list + train. var subsetEntries = new List(trainSubsetGlobalIdx.Length); @@ -1605,31 +1683,35 @@ private void RunFirstPassProteinFdr( List fullLibrary, OspreyConfig config) { - // Detected-peptide gate: targets passing peptide-level run FDR. - // Matches Rust pipeline.rs:3048 (e.run_peptide_qvalue <= run_fdr). - var detectedPeptides = new HashSet(StringComparer.Ordinal); + // Detected-peptide count for logging only; the core computation + + // propagation lives in ProteinFdr.RunFirstPassProteinFdr so the + // join-at-pass=2 rehydration path (PerFileRescoreTask) can run the + // same logic without duplicating it. + int detectedCount = 0; + var detectedTracker = new HashSet(StringComparer.Ordinal); foreach (var kvp in perFileEntries) { foreach (var entry in kvp.Value) { if (!entry.IsDecoy && entry.RunPeptideQvalue <= config.RunFdr) - detectedPeptides.Add(entry.ModifiedSequence); + detectedTracker.Add(entry.ModifiedSequence); } } + detectedCount = detectedTracker.Count; _ctx.LogInfo(string.Format( "[COUNT] First-pass detected peptides for protein FDR: {0} unique", - detectedPeptides.Count)); + detectedCount)); - var parsimony = ProteinFdr.BuildProteinParsimony( - fullLibrary, config.SharedPeptides, detectedPeptides); + ProteinFdr.RunFirstPassProteinFdr(perFileEntries, fullLibrary, config); - // Best peptide score across all files for picked-protein TDC. + // Recompute summary counters for log parity with the prior inline + // implementation. The static helper has already mutated entries + // via PropagateProteinQvalues; the parsimony / FDR objects below + // are rebuilt for logging + the diagnostic dump only. + var parsimony = ProteinFdr.BuildProteinParsimony( + fullLibrary, config.SharedPeptides, detectedTracker); var bestScores = ProteinFdr.CollectBestPeptideScores(perFileEntries); - - // First-pass gate is config.RunFdr exactly — matches Rust - // pipeline.rs:3062 (compute_protein_fdr at config.run_fdr). var proteinFdr = ProteinFdr.ComputeProteinFdr(parsimony, bestScores, config.RunFdr); - int nAtRunFdr = 0; foreach (var qv in proteinFdr.GroupQvalues.Values) { @@ -1647,11 +1729,6 @@ private void RunFirstPassProteinFdr( if (OspreyDiagnostics.ProteinFdrOnly) OspreyDiagnostics.ExitAfterDump(@"OSPREY_PROTEIN_FDR_ONLY"); } - - // Set RunProteinQvalue ONLY. Experiment-protein-q is set by the - // post-output Stage 8 protein FDR pass (Rust calls it second-pass). - ProteinFdr.PropagateProteinQvalues(perFileEntries, proteinFdr, - setRun: true, setExperiment: false); } } } diff --git a/pwiz_tools/OspreySharp/OspreySharp/Tasks/MergeNodeTask.cs b/pwiz_tools/OspreySharp/OspreySharp/Tasks/MergeNodeTask.cs index cf024c7214f..2762858a46b 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/Tasks/MergeNodeTask.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/Tasks/MergeNodeTask.cs @@ -102,6 +102,167 @@ public override bool Run(PipelineContext ctx) // Stage 8: Protein FDR (optional) if (config.ProteinFdr.HasValue) { + // Run 2nd-pass Percolator on the post-reconciliation + // entries when any 2nd-pass FDR sidecar is missing. + // Mirrors Rust pipeline.rs:4394-4468. After Stage 6 + // reconciliation, the entries' Features have been + // overwritten with rescored values, but their Scores + // are still the 1st-pass Percolator output (from + // FirstJoinTask). Without this 2nd-pass run, protein + // FDR (Stage 8) and the blib output would use stale + // 1st-pass scores; in the HPC distribution case the + // straight-through pipeline would silently lose ~25% + // of the precursors it produces -- the missing + // 2nd-pass step was the root cause behind the C# + // Stage 7 algorithmic divergence (issue: "Bug C"). + if (perFileParquetPaths.Count > 0 && config.InputFiles != null) + { + var inputByFileName = new Dictionary(StringComparer.Ordinal); + foreach (var inputFile in config.InputFiles) + inputByFileName[Path.GetFileNameWithoutExtension(inputFile)] = inputFile; + + // Surface any perFileEntries key that has no matching + // entry in config.InputFiles -- a silent skip here would + // hide a name-drift bug that the standard cross-impl gate + // (where keys always match) cannot catch. + var unmatchedKeys = perFileEntries + .Where(kvp => !inputByFileName.ContainsKey(kvp.Key)) + .Select(kvp => kvp.Key) + .ToList(); + if (unmatchedKeys.Count > 0) + { + ctx.LogWarning(string.Format( + "--join-at-pass=2: {0} perFileEntries key(s) have no matching " + + "config.InputFiles entry and will be skipped: [{1}]. This usually " + + "indicates an input-file rename or path drift between Stage 5 and " + + "Stage 7; the skipped files will not get a 2nd-pass sidecar.", + unmatchedKeys.Count, string.Join(", ", unmatchedKeys))); + } + + int missingPass2 = 0; + int totalFiles = 0; + foreach (var kvp in perFileEntries) + { + totalFiles++; + if (!inputByFileName.TryGetValue(kvp.Key, out string probeInput)) + continue; + if (!File.Exists(FdrScoresSidecar.Pass2Path(probeInput))) + missingPass2++; + } + if (missingPass2 > 0) + { + ctx.LogInfo(string.Format( + "--join-at-pass=2: {0}/{1} file(s) lack a 2nd-pass sidecar -- running " + + "second-pass FDR to compute scores from reconciled features " + + "(HPC distribution path; mirrors Rust pipeline.rs:4394-4468).", + missingPass2, totalFiles)); + ctx.LogInfo(string.Empty); + ctx.LogInfo("Second-pass FDR"); + // Reload PIN features from the reconciled parquets. + // PerFileScoringTask's bundle-hydration path + // explicitly nulls Features after stub load (see + // PerFileScoringTask.cs ~line 710) to keep + // PerFileRescoreTask.WriteReconciledParquet's + // "Features != null means this entry was rescored" + // criterion. That assumption was safe when Stage 7 + // didn't run Percolator -- with the Bug C 2nd-pass + // wired in below, we now need the 21-PIN features + // for SVM training, so pull them back from the + // post-Stage-6 reconciled parquet. The features + // there are the rescored values that Stage 6 wrote + // back, so they are the correct input for 2nd-pass + // Percolator. Mirrors Rust pipeline.rs:4209-4218 + // (run_search loads PIN features from parquet + // before second-pass FDR via the cache path). + var swReloadFeats = Stopwatch.StartNew(); + int nReloaded = 0; + foreach (var kvp in perFileEntries) + { + if (!perFileParquetPaths.TryGetValue(kvp.Key, out string parquetPath)) + { + // No first-join parquet was produced (or mapped) for this + // file. The {0} entries below will go into the second-pass + // Percolator with stale / null Features, which silently + // regresses 2nd-pass FDR -- log so the operator can detect + // an incomplete first-join hand-off. + ctx.LogWarning(string.Format( + "Second-pass FDR: no parquet path mapped for file '{0}' " + + "({1} entries will run with stale/null features). " + + "Check first-join output completeness.", + kvp.Key, kvp.Value.Count)); + continue; + } + List featRows; + try + { + featRows = ParquetScoreCache.LoadPinFeaturesFromParquet(parquetPath); + } + catch (Exception ex) + { + ctx.LogWarning(string.Format( + "Second-pass FDR: failed to reload PIN features from {0}: {1}", + parquetPath, ex.Message)); + continue; + } + int nMapped = 0; + foreach (var entry in kvp.Value) + { + int idx = (int)entry.ParquetIndex; + if (idx >= 0 && idx < featRows.Count) + { + entry.Features = featRows[idx]; + nMapped++; + } + } + // An entry whose ParquetIndex lies past the loaded row count + // is a stub/parquet mismatch (e.g., the first-join parquet + // was regenerated with fewer rows than the in-memory FDR + // stubs reference). Such entries silently keep their stale + // Features and corrupt 2nd-pass FDR; warn so the mismatch + // is visible. + if (nMapped < kvp.Value.Count) + { + ctx.LogWarning(string.Format( + "Second-pass FDR: file '{0}' parquet has {1} feature rows " + + "but {2} FDR entries reference it; {3} entries will run with " + + "stale/null features. Stub/parquet mismatch -- check first-join " + + "output integrity.", + kvp.Key, featRows.Count, kvp.Value.Count, kvp.Value.Count - nMapped)); + } + nReloaded += nMapped; + } + swReloadFeats.Stop(); + ctx.LogInfo(string.Format( + "[TIMING] Reloaded PIN features for {0} entries: {1:F1}s", + nReloaded, swReloadFeats.Elapsed.TotalSeconds)); + + var swPass2 = Stopwatch.StartNew(); + switch (config.FdrMethod) + { + case FdrMethod.Percolator: + FirstJoinTask.RunPercolatorFdr( + perFileEntries, fullLibrary, config, ctx, "Second-pass"); + break; + // Simple / Mokapot 2nd-pass paths intentionally + // not implemented yet -- the in-process pipeline's + // FirstJoinTask.RunFdr already covers Simple, and + // Mokapot is not used in OspreySharp's current + // scope. If those become relevant for an HPC chain, + // mirror the Rust dispatch in pipeline.rs:4424-4448. + default: + ctx.LogWarning(string.Format( + "Second-pass FDR: {0} is not supported in MergeNodeTask; " + + "skipping (protein FDR will run on first-pass scores)", + config.FdrMethod)); + break; + } + swPass2.Stop(); + ctx.LogInfo(string.Format( + "[STAGE-WALL] second-pass-fdr: {0:F1}s", + swPass2.Elapsed.TotalSeconds)); + } + } + // Persist post-Stage-6 per-file 2nd-pass FDR scores // BEFORE RunProteinFdr. The sidecar holds Score + // run/experiment precursor/peptide q-values + Pep + @@ -128,6 +289,22 @@ public override bool Run(PipelineContext ctx) foreach (var inputFile in config.InputFiles) inputByFileName[Path.GetFileNameWithoutExtension(inputFile)] = inputFile; + // Surface any perFileEntries key not in config.InputFiles + // -- a silent skip below would mean that file gets no + // .2nd-pass sidecar written and the next resume re-runs + // its second-pass FDR unnecessarily. + var unmatchedSidecarKeys = perFileEntries + .Where(kvp => !inputByFileName.ContainsKey(kvp.Key)) + .Select(kvp => kvp.Key) + .ToList(); + if (unmatchedSidecarKeys.Count > 0) + { + ctx.LogWarning(string.Format( + "2nd-pass sidecar write: {0} perFileEntries key(s) have no matching " + + "config.InputFiles entry and will be skipped: [{1}].", + unmatchedSidecarKeys.Count, string.Join(", ", unmatchedSidecarKeys))); + } + // Compute the task validity key once so each per-file // .MergeNode.osprey.task sidecar carries an identical // key. AnalysisPipeline.WriteTaskSidecars also writes @@ -198,13 +375,68 @@ public override bool Run(PipelineContext ctx) } } + // Re-load 2nd-pass FDR sidecar onto the post-compaction stub list. + // After the post-Stage-6 rehydration path, every stub still carries + // the 1st-pass q-values from RescoreHydration's 1st-pass sidecar + // overlay (PerFileScoringTask). The 2nd-pass q-values produced by + // Stage 6's reconciliation-aware rescore live in the + // .2nd-pass.fdr_scores.bin sidecar (or were just computed above and + // written to it). RunProteinFdr's detected_peptides gate filters on + // ExperimentPrecursorQvalue, which has to be the 2nd-pass value to + // match Rust pipeline.rs:4480-4494's reload-then-second-pass-FDR + // sequence. Without this reload, single-file --join-at-pass=2 runs + // include ~19 borderline peptides whose 1st-pass q-value passes + // <=1% but 2nd-pass q-value does not, producing a 1-protein delta + // in the Stage 7 picked-protein output cross-impl. + if (perFileParquetPaths.Count > 0 && config.InputFiles != null) + { + var inputByName = new Dictionary(StringComparer.Ordinal); + foreach (var inputFile in config.InputFiles) + inputByName[Path.GetFileNameWithoutExtension(inputFile)] = inputFile; + int filesReloaded = 0; + int filesMissing = 0; + foreach (var kvp in perFileEntries) + { + if (!inputByName.TryGetValue(kvp.Key, out string inputFile4)) + continue; + string pass2Path = FdrScoresSidecar.Pass2Path(inputFile4); + if (!File.Exists(pass2Path)) + { + filesMissing++; + continue; + } + var byEntryId = new Dictionary(kvp.Value.Count); + foreach (var e in kvp.Value) + byEntryId[e.EntryId] = e; + if (FdrScoresSidecar.TryReadOverlay( + pass2Path, byEntryId, FdrScoresSidecar.Pass.SecondPass)) + { + filesReloaded++; + } + else + { + filesMissing++; + ctx.LogWarning(string.Format( + "Failed to reload 2nd-pass FDR sidecar for {0} ({1}); " + + "protein FDR will use stale 1st-pass q-values", + kvp.Key, pass2Path)); + } + } + if (filesReloaded > 0) + { + ctx.LogInfo(string.Format( + "--join-at-pass=2: reloaded 2nd-pass FDR sidecar for {0}/{1} file(s) post-compaction", + filesReloaded, filesReloaded + filesMissing)); + } + } + ctx.LogInfo(string.Empty); ctx.LogInfo(string.Format(@"Running protein-level FDR at {0:P1}...", config.ProteinFdr.Value)); var swProtein = Stopwatch.StartNew(); RunProteinFdr(perFileEntries, fullLibrary, config); swProtein.Stop(); - ctx.LogInfo(string.Format(@"[TIMING] Protein FDR: {0:F1}s", + ctx.LogInfo(string.Format(@"[STAGE-WALL] stage7: {0:F1}s", swProtein.Elapsed.TotalSeconds)); } @@ -214,7 +446,7 @@ public override bool Run(PipelineContext ctx) var swBlib = Stopwatch.StartNew(); WriteBlibOutput(perFileEntries, fullLibrary, libraryById, config); swBlib.Stop(); - ctx.LogInfo(string.Format(@"[TIMING] Blib output: {0:F1}s", + ctx.LogInfo(string.Format(@"[STAGE-WALL] blib: {0:F1}s", swBlib.Elapsed.TotalSeconds)); return true; } @@ -269,6 +501,10 @@ private void RunProteinFdr( "[COUNT] Detected peptides for protein FDR: {0} unique", detectedPeptides.Count)); + // Cross-impl bisection dump (env-var-gated, no-op in production). + if (OspreyDiagnostics.DumpDetectedPeptides) + OspreyDiagnostics.WriteStage7DetectedPeptidesDump(detectedPeptides); + // Build protein parsimony var parsimony = ProteinFdr.BuildProteinParsimony( fullLibrary, config.SharedPeptides, detectedPeptides); diff --git a/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileRescoreTask.cs b/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileRescoreTask.cs index f5c45bcb805..f5baa57b9a3 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileRescoreTask.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileRescoreTask.cs @@ -27,6 +27,7 @@ using System.IO; using pwiz.OspreySharp.Chromatography; using pwiz.OspreySharp.Core; +using pwiz.OspreySharp.FDR; using pwiz.OspreySharp.FDR.Reconciliation; using pwiz.OspreySharp.IO; using pwiz.OspreySharp.Scoring; @@ -165,6 +166,76 @@ public override bool Run(PipelineContext ctx) var perFileScoring = ctx.GetTask(); _perFileEntries = perFileScoring.GetPerFileEntries(ctx); + // Hard short-circuit for --join-at-pass=2: every input parquet + // already has osprey.reconciled = "true" (asserted by + // ParquetScoreCache.CheckParquetMetadata when ExpectReconciledInput + // is set), so Stage 5 first-pass Percolator AND Stage 6 + // planning / rescore have ALREADY been performed upstream by + // the worker nodes that wrote those parquets. We must NOT + // touch FirstJoinTask here -- doing so transitively triggers + // FirstJoinTask.Run via EnsureHydrated, which re-runs Stage 5 + // first-pass Percolator from scratch on the reconciled parquets + // (producing wildly different action counts than the planner + // saw on the raw Stage 4 inputs) and then attempts a Stage 6 + // rescore that needs mzML files the merge node does not have + // (in production HPC, the merge node ships only sidecars + + // reconciled parquets, no mzMLs). MergeNodeTask is responsible + // for 2nd-pass Percolator (Bug C) and protein FDR + blib + // output starting from this hydrated, reconciled state. + // Mirrors Rust pipeline.rs:3313-3344 which gates the entire + // Stage 5+6 block on `!config.expect_reconciled_input`. + // + // Compaction still needs to run though: PerFileScoringTask's + // bundle-hydration path loads ALL entries from the parquet, + // including ones that failed first-pass FDR. FirstJoinTask's + // normal flow would run this compaction inline after first-pass + // Percolator (and we're skipping FirstJoinTask entirely here). + // Without it, MergeNodeTask's 2nd-pass Percolator would train + // on ~3x too many entries -- specifically the non-passing + // first-pass entries whose 1st-pass q-values are 1.0 -- and + // the SVM would learn a much worse decision boundary than the + // in-memory pipeline's, producing different per-precursor + // scores and different protein-FDR results. The compaction + // reads first-pass q-values that are already overlaid onto + // each entry from the .1st-pass.fdr_scores.bin sidecar by + // PerFileScoringTask's bundle hydration; no fresh FDR + // computation needed. + if (ctx.Config.ExpectReconciledInput) + { + var bundle = perFileScoring.GetRescoreInputs(ctx); + if (bundle != null) + { + // First-pass protein FDR BEFORE compaction. The 1st-pass FDR + // sidecar v3 already carries RunProteinQvalue from the original + // straight-through pipeline, but Rust pipeline.rs:4292 (gated by + // `!can_skip_fdr || config.expect_reconciled_input`) recomputes + // it inline in the --join-at-pass=2 path. The recompute uses the + // post-rehydration detected_peptides set + best_peptide_scores + // (which differ from the original write-time inputs whenever any + // upstream rebuild has nudged peptide q-values or score values + // even at the ULP level). Without this matching recompute on the + // C# side, the protein-rescue branch of compaction below sees + // slightly stale RunProteinQvalue values and the post-compaction + // detected_peptides set diverges from Rust by ~19 peptides on + // Stellar Single (1 protein delta at Stage 7). Only runs when + // protein FDR is enabled — the recompute is the protein-rescue + // input and is meaningless otherwise. Mirrors Rust pipeline.rs: + // 4292-4358. + if (ctx.Config.ProteinFdr.HasValue && bundle.PerFileEntries.Count > 0) + { + var fullLibrary = perFileScoring.GetFullLibrary(ctx); + ProteinFdr.RunFirstPassProteinFdr( + bundle.PerFileEntries, fullLibrary, ctx.Config); + } + var stats = RescoreCompaction.Apply(bundle, ctx.Config); + ctx.LogInfo(string.Format( + @"--join-at-pass=2 compaction: {0} -> {1} entries ({2} passing base_ids; {3} action(s) dropped)", + stats.EntriesBefore, stats.EntriesAfter, + stats.FirstPassBaseIds, stats.DroppedActions)); + } + return true; + } + // Self-gate: rescore + reconciliation only run when there is // planning state to act on AND the rescore hasn't already been // done upstream. State comes from either FirstJoinTask's @@ -178,9 +249,11 @@ public override bool Run(PipelineContext ctx) // the no-op alongside the no-state case. Probe-the-disk on // 2nd-pass sidecar presence replaces the prior // ExpectReconciledInput gate (Phase C: mechanism-driven, not - // flag-driven). Downstream MergeNodeTask still gets - // _perFileEntries via our accessor (falls through to the - // upstream reference). + // flag-driven) for the worker self-gate cases below; + // ExpectReconciledInput keeps the hard short-circuit above for + // the strict --join-at-pass=2 merge path. Downstream + // MergeNodeTask still gets _perFileEntries via our accessor + // (falls through to the upstream reference). var firstJoin = ctx.GetTask(); bool didPlan = firstJoin.DidPlan(ctx); var rescoreBundle = perFileScoring.GetRescoreInputs(ctx); @@ -204,6 +277,30 @@ public override bool Run(PipelineContext ctx) // preserve the valid sidecars for already-rescored files and // only invalidate the file(s) about to be re-rescored. + // Join file stems for the reconciled parquet metadata hash. + // In the in-process pipeline _perFileEntries has every file in + // the run; in worker mode (--join-at-pass=1 --no-join) it has + // a single file and the planner's full set comes from + // RescoreInputs.JoinFileStems (read from reconciliation.json + // v2+). Pass _perFileEntries keys when there's more than one; + // else fall through to the bundle's JoinFileStems. Null / + // empty means "let ExecuteRescore fall back to the + // InputFiles-derived hash" (preserves v1 behavior). + IReadOnlyList joinFileStems = null; + if (_perFileEntries != null && _perFileEntries.Count > 1) + { + var stems = new List(_perFileEntries.Count); + foreach (var kv in _perFileEntries) + stems.Add(kv.Key); + joinFileStems = stems; + } + else if (rescoreBundle != null + && rescoreBundle.JoinFileStems != null + && rescoreBundle.JoinFileStems.Count > 0) + { + joinFileStems = rescoreBundle.JoinFileStems; + } + var rescoreStats = ExecuteRescore( _perFileEntries, firstJoin.GetPerFileConsensusTargets(ctx), @@ -213,7 +310,8 @@ public override bool Run(PipelineContext ctx) firstJoin.GetPerFileGapFillForRescore(ctx), perFileScoring.GetPerFileParquetPaths(ctx), perFileScoring.GetFullLibrary(ctx), - ctx.Config); + ctx.Config, + joinFileStems); ctx.LogInfo(string.Format( @"Stage 6 rescore: {0} entries re-scored ({1} reconciliation actions executed)", rescoreStats.TotalRescored, rescoreStats.TotalReconciliation)); @@ -267,7 +365,8 @@ public override bool Run(PipelineContext ctx) /// Mirrors Rust pipeline.rs:3050-3110. /// private void WriteReconciledParquet(string parquetPath, List fdrEntries, - string fileName, List fullLibrary, OspreyConfig config) + string fileName, List fullLibrary, OspreyConfig config, + IReadOnlyList joinFileStems) { // 1. Reload the original parquet's per-row state. List fullEntries; @@ -336,13 +435,27 @@ private void WriteReconciledParquet(string parquetPath, List fdrEntrie // 5. Reconciliation metadata (mirrors Rust // build_reconciled_metadata). osprey.version is what the // next reload's CacheValidity check compares against. + // The reconciliation_hash must be the JOIN-wide hash + // (over every file in the planner step), not the worker's + // single-file InputFiles hash; without that, a worker + // rescoring a single parquet stamps a single-file hash + // that the downstream --join-at-pass=2 merge node rejects + // on hash mismatch. The join file stems come from the + // planner's reconciliation.json (v2+) via + // RescoreInputs.JoinFileStems; fall back to config-derived + // stems when the caller didn't pass any (in-process + // pipeline where config.InputFiles already has all files, + // or v1 backward compat). + string reconciliationHash = (joinFileStems != null && joinFileStems.Count > 0) + ? config.ReconciliationParameterHashForStems(joinFileStems) + : config.ReconciliationParameterHash(); var metadata = new Dictionary { { @"osprey.version", Program.VERSION }, { @"osprey.search_hash", config.SearchParameterHash() }, { @"osprey.library_hash", config.LibraryIdentityHash() }, { @"osprey.reconciled", @"true" }, - { @"osprey.reconciliation_hash", config.ReconciliationParameterHash() }, + { @"osprey.reconciliation_hash", reconciliationHash }, }; try @@ -395,7 +508,8 @@ internal RescoreStats ExecuteRescore( IReadOnlyDictionary> perFileGapFill, IReadOnlyDictionary perFileParquetPaths, List fullLibrary, - OspreyConfig config) + OspreyConfig config, + IReadOnlyList joinFileStems = null) { // Pre-group reconciliation actions by file. Mirrors the Rust // pre-grouping at pipeline.rs:2719-2744 -- a single pass over @@ -442,9 +556,12 @@ internal RescoreStats ExecuteRescore( // --input-scores parquet stems by Program.Main; for in-process, // it's the user's -i mzML list. Either way the stem matches // the file_name keys in perFileEntries. - var fileNameToIdx = new Dictionary(); + var fileNameToIdx = new Dictionary(StringComparer.Ordinal); for (int i = 0; i < config.InputFiles.Count; i++) - fileNameToIdx[Path.GetFileNameWithoutExtension(config.InputFiles[i])] = i; + { + string stem = Path.GetFileNameWithoutExtension(config.InputFiles[i]) ?? string.Empty; + fileNameToIdx[stem] = i; + } int totalRescored = 0; int totalGapCwt = 0; @@ -721,12 +838,25 @@ internal RescoreStats ExecuteRescore( int nGapForced = 0; if (gapFillTargets.Count > 0) { - // Build target+decoy id set from gap_fill_targets. + // Build gap-fill library subset (targets only). + // + // Decoys are intentionally excluded from gap-fill: forcing a + // random decoy sequence to be scored at the target's + // consensus RT has no biological basis (decoys are not + // expected to co-elute with their paired target), and the + // 1st-pass parquet already has a score for every decoy at + // its own natural-but-best peak. Gap-filling decoys also + // re-scored them at consensus RT and APPENDED a second + // parquet row alongside the existing 1st-pass row, + // producing exact-duplicate rows in the reconciled parquet. + // Those duplicates cascaded into different max-per-modseq + // aggregations cross-impl and a 1.1e-4 group_qvalue drift + // on Astral 3-file. Targets are still gap-filled because + // they were missing from this file by definition. var gapFillIds = new HashSet(); foreach (var gf in gapFillTargets) { gapFillIds.Add(gf.TargetEntryId); - gapFillIds.Add(gf.DecoyEntryId); } var gapFillLibrary = new List(gapFillIds.Count); foreach (var libEntry in fullLibrary) @@ -787,10 +917,9 @@ internal RescoreStats ExecuteRescore( cwtHitIds = new HashSet(); } - // Pass 2: Forced integration for entries CWT missed. - // For each gap-fill target, check both the target_id - // and decoy_id; either or both may have missed the CWT - // pass. + // Pass 2: Forced integration for targets CWT missed. + // Decoys are intentionally excluded from gap-fill (see + // gapFillIds build above). var forcedOverrides = new Dictionary(); var forcedIds = new HashSet(); foreach (var gf in gapFillTargets) @@ -802,11 +931,6 @@ internal RescoreStats ExecuteRescore( forcedOverrides[gf.TargetEntryId] = (gf.ExpectedRt, start, end); forcedIds.Add(gf.TargetEntryId); } - if (!cwtHitIds.Contains(gf.DecoyEntryId)) - { - forcedOverrides[gf.DecoyEntryId] = (gf.ExpectedRt, start, end); - forcedIds.Add(gf.DecoyEntryId); - } } if (forcedOverrides.Count > 0) @@ -861,7 +985,7 @@ internal RescoreStats ExecuteRescore( File.Exists(parquetPath)) { WriteReconciledParquet(parquetPath, fdrEntries, fileName, - fullLibrary, config); + fullLibrary, config, joinFileStems); // Per-file resume sidecar: write next to the // reconciled parquet so a subsequent invocation with @@ -937,8 +1061,17 @@ private void LoadSpectraForRescore(string inputFile, string fileName, /// /// Load MS2 + MS1 mass calibrations and the original Stage-4 RT /// calibration MAD from the sibling .calibration.json that - /// Stage 2 wrote. Returns uncalibrated results / null MAD if the - /// file is missing or the relevant section is absent. + /// Stage 2 wrote. Throws if the + /// calibration sidecar is missing or unreadable -- Stage 6 + /// requires the Stage 1-4 calibration to rescore, and silently + /// falling back to uncalibrated would mask a real configuration + /// error (the worker's output would diverge from the + /// straight-through pipeline's output). Mirrors the hard-error + /// behavior in Rust run_rescore at + /// osprey/crates/osprey/src/rescore.rs. Individual calibration + /// sections (Ms1Calibration / Ms2Calibration / RtMad) may still + /// be absent within the file; those leave the corresponding + /// out-param at its uncalibrated / null default. /// private void LoadMassCalibrations(string inputFile, out MzCalibrationResult ms2Cal, out MzCalibrationResult ms1Cal, @@ -950,10 +1083,20 @@ private void LoadMassCalibrations(string inputFile, string parent = Path.GetDirectoryName(Path.GetFullPath(inputFile)); if (string.IsNullOrEmpty(parent)) - return; + { + throw new InvalidDataException(string.Format( + "LoadMassCalibrations: cannot derive sidecar directory from input path `{0}`. " + + "Stage 6 needs to read the Stage 1-4 calibration sidecar; without it the " + + "worker would silently produce uncalibrated rescore output.", inputFile)); + } string calPath = CalibrationIO.CalibrationPathForInput(inputFile, parent); if (!File.Exists(calPath)) - return; + { + throw new InvalidDataException(string.Format( + "LoadMassCalibrations: required calibration JSON not found at `{0}` " + + "(input file: `{1}`). Stage 6 needs the Stage 1-4 calibration sidecar to " + + "rescore. Run Stages 1-4 first or fix the path.", calPath, inputFile)); + } CalibrationParams calParams; try @@ -962,9 +1105,10 @@ private void LoadMassCalibrations(string inputFile, } catch (Exception ex) { - _ctx.LogWarning(string.Format( - "Failed to load calibration JSON {0}: {1}", calPath, ex.Message)); - return; + throw new InvalidDataException(string.Format( + "LoadMassCalibrations: failed to read calibration JSON `{0}`: {1}. The file " + + "exists but could not be parsed -- check that it was written by a matching " + + "OspreySharp version.", calPath, ex.Message), ex); } if (calParams.Ms2Calibration != null && calParams.Ms2Calibration.Calibrated) diff --git a/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileScoringTask.cs b/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileScoringTask.cs index cc2f168d44f..edc8bd5a4e3 100644 --- a/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileScoringTask.cs +++ b/pwiz_tools/OspreySharp/OspreySharp/Tasks/PerFileScoringTask.cs @@ -83,6 +83,18 @@ private class CalibrationPassResult public RTCalibrationStats Stats; public MzCalibrationResult Ms1Calibration; public MzCalibrationResult Ms2Calibration; + // Total matches scored in this pass (before any q-value or S/N + // filtering). Plumbed into CalibrationMetadata.NumSampledPrecursors + // for parity with Rust's accumulated_matches.len(). + public int MatchCount; + // (lib_rt, measured_rt) pairs that were actually fed to the LOESS + // fit for this pass. Exposed so the caller can emit the + // OSPREY_DUMP_LOESS_INPUT diagnostic only for the pass whose + // calibration is actually used (pass 1 always; pass 2 only on + // acceptance) -- mirroring Rust pipeline.rs's "dump reflects + // the calibration actually used" semantics. + public double[] LibRts; + public double[] MeasuredRts; } public override string Name => @"PerFileScoring"; @@ -969,6 +981,11 @@ private List ProcessFile( RTCalibration rtCalibration = null; MzCalibrationResult ms2Cal = MzCalibrationResult.Uncalibrated(); MzCalibrationResult ms1Cal = MzCalibrationResult.Uncalibrated(); + // Total matches scored during pass 1 of calibration; threaded + // into CalibrationMetadata.NumSampledPrecursors to match Rust's + // accumulated_matches.len() (Stellar Single: 192289). Stays 0 + // when calibration is loaded from a cached JSON. + int numSampledPrecursorsForMetadata = 0; // BISECT: load Rust's calibration JSON instead of computing our own. // This eliminates calibration noise from the feature comparison. @@ -1023,7 +1040,7 @@ private List ProcessFile( var swCal = Stopwatch.StartNew(); rtCalibration = RunCalibration( fullLibrary, spectra, ms1Spectra, context, - out ms1Cal, out ms2Cal); + out ms1Cal, out ms2Cal, out numSampledPrecursorsForMetadata); swCal.Stop(); int nPoints = rtCalibration != null ? rtCalibration.Stats().NPoints : 0; _ctx.LogInfo(string.Format( @@ -1048,8 +1065,16 @@ private List ProcessFile( Metadata = new CalibrationMetadata { CalibrationSuccessful = rtCalibration != null, - NumConfidentPeptides = 0, - NumSampledPrecursors = 0, + // Match Rust's CalibrationMetadata field semantics + // (osprey/src/pipeline.rs:1144-1145): + // num_confident_peptides = LOESS points used by the + // fit actually saved (post-S/N filter) + // num_sampled_precursors = total matches scored + // during sampling (pre any q-value / S/N filter) + NumConfidentPeptides = rtCalibration != null + ? rtCalibration.Stats().NPoints + : 0, + NumSampledPrecursors = numSampledPrecursorsForMetadata, Timestamp = DateTime.UtcNow.ToString("o") }, Ms1Calibration = MzCalibrationJson.FromResult(ms1Cal), @@ -1305,9 +1330,13 @@ private RTCalibration RunCalibration( List ms1Spectra, ScoringContext context, out MzCalibrationResult ms1Calibration, - out MzCalibrationResult ms2Calibration) + out MzCalibrationResult ms2Calibration, + out int numSampledPrecursors) { var config = context.Config; + // Default to 0 so early returns / exception paths leave the + // metadata caller in a known state. Overwritten on success. + numSampledPrecursors = 0; _ctx.LogInfo("Running RT calibration..."); // Calculate library and mzML RT ranges @@ -1432,6 +1461,24 @@ private RTCalibration RunCalibration( return null; } + // Pass 1 succeeded -- emit the OSPREY_DUMP_LOESS_INPUT diagnostic + // for the pass-1 fit unconditionally. If pass 2 is later accepted + // it will overwrite this with the pass-2 pairs; if pass 2 is + // rejected (or never runs) the pass-1 dump stays, matching the + // calibration actually used. Mirrors Rust pipeline.rs. + if (OspreyDiagnostics.DumpLoessInput) + { + OspreyDiagnostics.WriteLoessInputDump(1, pass1.LibRts, pass1.MeasuredRts); + if (OspreyDiagnostics.LoessInputOnly) + OspreyDiagnostics.ExitAfterDump("OSPREY_LOESS_INPUT_ONLY"); + } + + // Match Rust accumulated_matches.len() semantics: report pass 1's + // total scored matches (before any q-value / S/N filtering). + // Pass 2 is a refinement using narrowed RT tolerance and does not + // change the sampled-precursor count. + numSampledPrecursors = pass1.MatchCount; + // === Iterative calibration refinement (2-pass) === // Mirrors Rust pipeline.rs:714-839. // MAD * 1.4826 ~ SD for a normal distribution; 3* that covers ~99.7%. @@ -1484,6 +1531,15 @@ private RTCalibration RunCalibration( // by more than 1% (matches Rust pipeline.rs:811). if (pass2.Stats.RSquared >= pass1.Stats.RSquared * 0.99) { + // Overwrite the OSPREY_DUMP_LOESS_INPUT dump with + // pass 2's points so the diagnostic reflects the + // calibration actually being used. Mirrors Rust + // pipeline.rs; only fires on acceptance. + if (OspreyDiagnostics.DumpLoessInput) + { + OspreyDiagnostics.WriteLoessInputDump( + 2, pass2.LibRts, pass2.MeasuredRts); + } ms1Calibration = pass2.Ms1Calibration; ms2Calibration = pass2.Ms2Calibration; return pass2.Calibration; @@ -1726,6 +1782,10 @@ private CalibrationPassResult RunCalibrationScoringPass( // filter). This is the same set used for RT calibration points. var allMs1Errors = new List(); var allMs2Errors = new List(); + // Track contributing matches in a stable list so a cross-impl + // bisection dump can replay them in the same order the + // calibration accumulator sees their errors. + var contributingMatches = new List(); foreach (var m in matchArray) { if (m.Ms2MassErrors == null || m.IsDecoy || m.QValue > CAL_FDR_THRESHOLD) @@ -1736,6 +1796,13 @@ private CalibrationPassResult RunCalibrationScoringPass( if (m.Ms1Error.HasValue) allMs1Errors.Add(m.Ms1Error.Value); allMs2Errors.AddRange(m.Ms2MassErrors); + contributingMatches.Add(m); + } + if (OspreyDiagnostics.DumpMs2CalErrors) + { + OspreyDiagnostics.WriteMs2CalErrorsDump(contributingMatches); + if (OspreyDiagnostics.Ms2CalErrorsOnly) + OspreyDiagnostics.ExitAfterDump("OSPREY_MS2_CAL_ERRORS_ONLY"); } string unitStr = config.FragmentTolerance.Unit == ToleranceUnit.Ppm ? "ppm" : "Th"; var ms1Cal = MzCalibration.CalculateSingleLevel(allMs1Errors.ToArray(), unitStr); @@ -1764,15 +1831,12 @@ private CalibrationPassResult RunCalibrationScoringPass( ClassicalRobustIterations = OspreyEnvironment.LoessClassicalRobust }; - // Cross-implementation diagnostic: dump the (lib_rt, measured_rt) pairs - // fed to LOESS. Used to verify Rust and C# see identical inputs - // before LOESS fitting. - if (OspreyDiagnostics.DumpLoessInput) - { - OspreyDiagnostics.WriteLoessInputDump(passNumber, libRts, measuredRts); - if (OspreyDiagnostics.LoessInputOnly) - OspreyDiagnostics.ExitAfterDump("OSPREY_LOESS_INPUT_ONLY"); - } + // NOTE: the OSPREY_DUMP_LOESS_INPUT diagnostic is emitted by + // the caller (RunCalibration) so that pass 2 only dumps when + // it is actually accepted -- mirrors Rust pipeline.rs (and + // matches the PR #42 semantics of "dump reflects the + // calibration actually used"). The pairs are returned via + // CalibrationPassResult.LibRts/MeasuredRts. var calibrator = new RTCalibrator(calibratorConfig); var rtCal = calibrator.Fit(libRts, measuredRts); @@ -1791,6 +1855,9 @@ private CalibrationPassResult RunCalibrationScoringPass( Stats = stats, Ms1Calibration = ms1Cal, Ms2Calibration = ms2Cal, + MatchCount = matchArray.Length, + LibRts = libRts, + MeasuredRts = measuredRts, }; } catch (Exception ex)