OspreySharp: close C# Stage 7 algorithmic gap + port Rust v2 reconciliation envelope + parity-sprint perf work

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<T> 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];

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

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];

pwiz_tools/OspreySharp/OspreySharp.Core/OspreyConfig.cs

@@ -387,6 +387,33 @@ internal static string EscapeForRustDebug(string s)
         /// hash invariant to invocation order).
         /// </summary>
         public string ReconciliationParameterHash()
+        {
+            var stems = new List<string>(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);
+        }
+
+        /// <summary>
+        /// Compute the reconciliation parameter hash for an explicit set of
+        /// file stems. Used by per-file Stage 6 rescore workers, whose
+        /// <see cref="InputFiles"/> only carries this worker's single
+        /// parquet — the hash that the downstream <c>--join-at-pass=2</c>
+        /// merge node expects is computed over ALL files in the join, so
+        /// the worker must read the full set from the planner's
+        /// <c>reconciliation.json</c> envelope and pass it in here. The
+        /// stems are sorted + deduped internally so the hash is invariant
+        /// to caller ordering. Mirrors Rust
+        /// <c>OspreyConfig::reconciliation_parameter_hash_for_stems</c>.
+        /// </summary>
+        public string ReconciliationParameterHashForStems(IReadOnlyList<string> 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<String> yields ["a", "b"] with the
                 // brackets and double-quoted, comma-space-separated values.
-                var stems = new List<string>(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<string>(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++)
                 {

pwiz_tools/OspreySharp/OspreySharp.FDR/FdrDiagnostics.cs

@@ -0,0 +1,125 @@
+/*
+ * Original author: Brendan MacLean <brendanx .at. uw.edu>,
+ *                  MacCoss Lab, Department of Genome Sciences, UW
+ * AI assistance: Claude Code (Claude Opus 4.7) <noreply .at. anthropic.com>
+ *
+ * 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";
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        public static readonly bool DumpStage7Winners = IsOne(@"OSPREY_DUMP_STAGE7_WINNERS");
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        public static readonly bool DumpBestPeptideScores = IsOne(@"OSPREY_DUMP_BEST_PEPTIDE_SCORES");
+
+        /// <summary>
+        /// Write cs_stage7_winners.tsv. Caller passes a (score, is_decoy)
+        /// tuple list in sort order plus the parallel q-value arrays. Check
+        /// <see cref="DumpStage7Winners"/> first to skip the LINQ projection
+        /// on the disabled-dump path.
+        /// </summary>
+        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])));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Write cs_best_peptide_scores.tsv. Rows sorted by
+        /// modified_sequence for stable cross-impl diff.
+        /// </summary>
+        public static void WriteBestPeptideScoresDump(Dictionary<string, PeptideScore> best)
+        {
+            const string path = @"cs_best_peptide_scores.tsv";
+            var inv = CultureInfo.InvariantCulture;
+            var keys = new List<string>(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)));
+                }
+            }
+        }
+    }
+}