From 3bcf4c38e11721d6073bf123ca431890963aea47 Mon Sep 17 00:00:00 2001 From: trishorts Date: Wed, 10 Jun 2026 07:34:07 -0500 Subject: [PATCH 1/4] feat(parallelsearch): spectral-library (.msl) peptide source with learned mass+RT pre-filter Use mzLib .msl spectral libraries as the transient-search peptide source instead of re-digesting FASTAs. Adds a lean (fragment-less) reader that fragments on the fly, a candidate pre-filter that learns the precursor-mass tolerance and an RT (iRT->observed) calibration from the confident base PSMs and indexes the library by mass and retention time, a single merged "index of indexes" so thousands of databases load as one file, a load-ahead producer/consumer pipeline, and a CPU fragment-matching abstraction (ISpectralScorer / CpuSpectralScorer / struct-of-arrays ScoringBatch). On the SARS spike-in (12,973 virus DBs, 3 raws): ~40% faster wall-clock end-to-end (42 vs 73 min), ~5x less search-engine CPU, the merged load collapses 12,973 file-opens (~558 s) to a single ~23 s read, and zero specificity cost (every confident MSL peptide is also found by the exhaustive FASTA search). Bumps the mzLib dependency to 9.9.23 for the .msl reader types. Tests: 7 unit tests for the candidate filter (mass/RT window, ppm + 0.01 Da margin, iRT==0 mass-only), the merged "db|accession" parse, and the lower-bound helper; both read loops call the shared predicate. Co-Authored-By: Claude Opus 4.8 (1M context) --- MetaMorpheus/CMD/CMD.csproj | 4 +- MetaMorpheus/EngineLayer/EngineLayer.csproj | 4 +- .../ParallelSearch/MslPeptideReader.cs | 324 +++++++++++ .../Scoring/CpuSpectralScorer.cs | 80 +++ .../ParallelSearch/Scoring/ISpectralScorer.cs | 95 +++ .../Scoring/SpectralScorerProvider.cs | 26 + .../Scoring/SpectralScoringData.cs | 147 +++++ .../TransientClassicSearchEngine.cs | 322 ++++++++--- MetaMorpheus/GUI/GUI.csproj | 4 +- MetaMorpheus/GuiFunctions/GuiFunctions.csproj | 4 +- .../IO/StatisticalTestResultFile.cs | 10 + .../IO/TestSummaryResultFile.cs | 9 + .../ParallelSearch/ParallelSearchTask.cs | 546 ++++++++++++++++-- .../AnomalyDetectionFeatureBuilder.cs | 11 +- MetaMorpheus/TaskLayer/TaskLayer.csproj | 4 +- .../MslCandidateFilterTests.cs | 111 ++++ MetaMorpheus/Test/Test.csproj | 4 +- 17 files changed, 1563 insertions(+), 142 deletions(-) create mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs create mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs create mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs create mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs create mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs create mode 100644 MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs diff --git a/MetaMorpheus/CMD/CMD.csproj b/MetaMorpheus/CMD/CMD.csproj index a60989660..94cfdf487 100644 --- a/MetaMorpheus/CMD/CMD.csproj +++ b/MetaMorpheus/CMD/CMD.csproj @@ -1,4 +1,4 @@ - + Exe @@ -24,7 +24,7 @@ - + diff --git a/MetaMorpheus/EngineLayer/EngineLayer.csproj b/MetaMorpheus/EngineLayer/EngineLayer.csproj index 06e043759..1baac2c68 100644 --- a/MetaMorpheus/EngineLayer/EngineLayer.csproj +++ b/MetaMorpheus/EngineLayer/EngineLayer.csproj @@ -1,4 +1,4 @@ - + net8.0 @@ -29,7 +29,7 @@ - + diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs b/MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs new file mode 100644 index 000000000..01e2f83e0 --- /dev/null +++ b/MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs @@ -0,0 +1,324 @@ +using System; +using System.Collections.Generic; +using System.Text; +using Chemistry; +using Omics; +using Omics.Fragmentation; +using Omics.Fragmentation.Peptide; +using Omics.SpectralMatch.MslSpectralLibrary; +using Proteomics; +using Proteomics.ProteolyticDigestion; +using Readers.SpectralLibrary; + +namespace EngineLayer.ParallelSearch +{ + /// + /// Reads a MetaMorpheus .msl spectral library as a PEPTIDE SOURCE for the parallel search, + /// using the library the way the binary format was designed for: + /// keeps fragments on disk; the precursor index (mass / charge / iRT metadata, NO fragment I/O) is + /// filtered on TWO orthogonal axes — precursor mass AND retention time — and only the surviving + /// CANDIDATES have their fragments fetched on demand. + /// + /// Retention time is the axis that makes the filter selective: precursor mass alone is degenerate when + /// tens of thousands of dense experimental precursors saturate the mass window, but a peptide's predicted + /// elution (stored as Chronologer iRT in the .msl, calibrated to this run from the base search) collides + /// with far fewer scans. Entries with iRT==0 (predictor could not place them) fall back to mass-only so + /// they are never lost. + /// + /// Each candidate's stored float fragments are turned into s directly (no + /// re-fragmentation). Accession handling: one shared per accession (concatenation + /// of its candidate peptides with per-peptide offsets) so parsimony/coverage stay correct and in-bounds. + /// + public static class MslPeptideReader + { + private const double MassPreFilterMarginDa = 0.01; + + /// RT↔iRT calibration and learned precursor tolerance applied to the candidate pre-filter. + public readonly struct CandidatePriors + { + public readonly double[] SortedScanMasses; // experimental precursor masses, ASCENDING + public readonly double[] ScanRetentionTimes; // RT of each scan, SAME order as SortedScanMasses + public readonly double PrecursorTolPpm; // learned precursor tolerance (ppm) + public readonly double RtSlope, RtIntercept; // observedRT = RtSlope*iRT + RtIntercept + public readonly double RtWindowMin; // +/- window in observed-RT minutes (k * residualSD) + + public CandidatePriors(double[] sortedScanMasses, double[] scanRetentionTimes, + double precursorTolPpm, double rtSlope, double rtIntercept, double rtWindowMin) + { + SortedScanMasses = sortedScanMasses; + ScanRetentionTimes = scanRetentionTimes; + PrecursorTolPpm = precursorTolPpm; + RtSlope = rtSlope; + RtIntercept = rtIntercept; + RtWindowMin = rtWindowMin; + } + } + + private readonly struct PendingPeptide + { + public readonly string FullSequence; + public readonly int Start; + public readonly int Length; + public readonly List Fragments; + public PendingPeptide(string fullSequence, int start, int length, List fragments) + { FullSequence = fullSequence; Start = start; Length = length; Fragments = fragments; } + } + + private sealed class AccessionGroup + { + public readonly StringBuilder Sequence = new(); + public readonly List Peptides = new(); + public bool IsDecoy; + } + + public static List<(IBioPolymerWithSetMods Peptide, List Fragments)> ReadPeptides( + string mslPath, string databaseName, CandidatePriors priors) + { + var mods = GlobalVariables.AllModsKnownDictionary; + var digestionParams = new DigestionParams("trypsin", maxMissedCleavages: 0); + string fallbackAccession = $"{databaseName}_UNKNOWN"; + var groups = new Dictionary(); + + int totalEntries = 0, massCandidates = 0, massRtCandidates = 0; + var sw = System.Diagnostics.Stopwatch.StartNew(); + + long loadMs, getEntryTicks = 0, buildTicks = 0; + var swPhase = System.Diagnostics.Stopwatch.StartNew(); + var candidateIdx = new List(); + using (var library = MslLibrary.LoadIndexOnly(mslPath)) + { + loadMs = swPhase.ElapsedMilliseconds; // cost of LoadIndexOnly (metadata read + sort) + + // PASS 1 — filter on mass + RT using ONLY the in-memory index (no fragment I/O). + foreach (var idx in library.QueryMzWindow(float.MinValue, float.MaxValue)) + { + totalEntries++; + if (IsCandidate(idx.PrecursorMz, idx.Charge, idx.Irt, priors, out bool massMatched)) + { + massRtCandidates++; + candidateIdx.Add(idx.PrecursorIdx); + } + if (massMatched) + massCandidates++; + } + + // PASS 2 — fetch fragments in ON-DISK (PrecursorIdx/write) order, not the m/z order the + // window query yielded. Fragment blocks are stored in PrecursorIdx order, so a sorted walk + // reads the file front-to-back (sequential, OS read-ahead) instead of one random seek per + // candidate — the dominant reader cost. Result order doesn't affect IDs/scores; the shared + // protein concatenation just follows this order. + candidateIdx.Sort(); + foreach (int pid in candidateIdx) + { + long t0 = swPhase.ElapsedTicks; + MslLibraryEntry entry = library.GetEntry(pid); + getEntryTicks += swPhase.ElapsedTicks - t0; + if (entry is null) continue; + + string fullSequence = entry.FullSequence; + string baseSequence = IBioPolymerWithSetMods.GetBaseSequenceFromFullSequence(fullSequence); + string accession = string.IsNullOrEmpty(entry.ProteinAccession) ? fallbackAccession : entry.ProteinAccession; + + if (!groups.TryGetValue(accession, out var group)) + { + group = new AccessionGroup { IsDecoy = entry.IsDecoy }; + groups[accession] = group; + } + + int start = group.Sequence.Length + 1; + group.Sequence.Append(baseSequence); + long t1 = swPhase.ElapsedTicks; + // Lean libraries store no fragments → return null so the engine fragments the peptide on + // the fly (cheap, and matches the search's exact dissociation/precision). Full libraries + // build Products from the stored float ions. + List built = entry.MatchedFragmentIons is { Count: > 0 } + ? BuildProducts(entry.MatchedFragmentIons) + : null; + buildTicks += swPhase.ElapsedTicks - t1; + group.Peptides.Add(new PendingPeptide(fullSequence, start, baseSequence.Length, built)); + } + } + double getEntryMs = getEntryTicks * 1000.0 / System.Diagnostics.Stopwatch.Frequency; + double buildMs = buildTicks * 1000.0 / System.Diagnostics.Stopwatch.Frequency; + + // Per-database candidate-filter diagnostics. Opt-in (MM_PARALLELSEARCH_DIAG=1) so it is + // available when profiling the .msl reader but silent in normal runs. + if (Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_DIAG") == "1") + Console.WriteLine($"MSL_RT {databaseName}: entries={totalEntries} massCand={massCandidates} " + + $"({Pct(massCandidates, totalEntries)}%) massRtCand={massRtCandidates} ({Pct(massRtCandidates, totalEntries)}%) " + + $"load={loadMs}ms getEntry={getEntryMs:F0}ms build={buildMs:F0}ms " + + $"tolPpm={priors.PrecursorTolPpm:F1} rtWin=+/-{priors.RtWindowMin:F1}min ms={sw.ElapsedMilliseconds}"); + + var result = new List<(IBioPolymerWithSetMods, List)>(); + foreach (var kvp in groups) + { + var protein = new Protein(kvp.Value.Sequence.ToString(), kvp.Key, isDecoy: kvp.Value.IsDecoy); + foreach (var pending in kvp.Value.Peptides) + { + var peptide = new PeptideWithSetModifications( + pending.FullSequence, mods, p: protein, digestionParams: digestionParams, + oneBasedStartResidueInProtein: pending.Start, + oneBasedEndResidueInProtein: pending.Start + pending.Length - 1, + missedCleavages: 0); + result.Add((peptide, pending.Fragments)); + } + } + return result; + } + + /// + /// Reads a MERGED index (many databases in one .msl, each entry's accession stamped "db|accession") + /// with ONE open, runs the mass+RT candidate filter ONCE over all entries, and returns the surviving + /// candidate peptides GROUPED BY their source database. The per-database lists are then searched + /// INDEPENDENTLY (each against the shared base PSMs in its own copy) exactly as in the per-file path — + /// databases never compete; this only collapses 1000s of file opens into one shared in-memory index. + /// + public static Dictionary Fragments)>> + ReadCandidatesGroupedByDatabase(string mergedMslPath, CandidatePriors priors) + { + var mods = GlobalVariables.AllModsKnownDictionary; + var digestionParams = new DigestionParams("trypsin", maxMissedCleavages: 0); + + // dbTag -> (accession-within-db -> shared-protein layout) + var byDb = new Dictionary>(); + var candidateIdx = new List(); + int totalEntries = 0, massRtCandidates = 0; + var swPhase = System.Diagnostics.Stopwatch.StartNew(); + long loadMs; + + using (var library = MslLibrary.LoadIndexOnly(mergedMslPath)) + { + loadMs = swPhase.ElapsedMilliseconds; + // PASS 1 — mass+RT filter over the whole merged index (no fragment I/O). + foreach (var idx in library.QueryMzWindow(float.MinValue, float.MaxValue)) + { + totalEntries++; + if (IsCandidate(idx.PrecursorMz, idx.Charge, idx.Irt, priors, out _)) + { + massRtCandidates++; + candidateIdx.Add(idx.PrecursorIdx); + } + } + + // PASS 2 — fetch candidates in on-disk order (sequential), group by source database. + candidateIdx.Sort(); + foreach (int pid in candidateIdx) + { + MslLibraryEntry entry = library.GetEntry(pid); + if (entry is null) continue; + + var (dbTag, accession) = ParseDbTagAndAccession(entry.ProteinAccession); + + string fullSequence = entry.FullSequence; + string baseSequence = IBioPolymerWithSetMods.GetBaseSequenceFromFullSequence(fullSequence); + + if (!byDb.TryGetValue(dbTag, out var accGroups)) + { + accGroups = new Dictionary(); + byDb[dbTag] = accGroups; + } + if (!accGroups.TryGetValue(accession, out var group)) + { + group = new AccessionGroup { IsDecoy = entry.IsDecoy }; + accGroups[accession] = group; + } + int start = group.Sequence.Length + 1; + group.Sequence.Append(baseSequence); + List built = entry.MatchedFragmentIons is { Count: > 0 } ? BuildProducts(entry.MatchedFragmentIons) : null; + group.Peptides.Add(new PendingPeptide(fullSequence, start, baseSequence.Length, built)); + } + } + + if (Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_DIAG") == "1") + Console.WriteLine($"MSL_MERGED: entries={totalEntries} candidates={massRtCandidates} " + + $"({Pct(massRtCandidates, totalEntries)}%) databases={byDb.Count} load={loadMs}ms total={swPhase.ElapsedMilliseconds}ms"); + + // Materialize per-database peptide lists (each database is searched independently downstream). + var result = new Dictionary)>>(byDb.Count); + foreach (var dbKvp in byDb) + { + var list = new List<(IBioPolymerWithSetMods, List)>(); + foreach (var accKvp in dbKvp.Value) + { + var protein = new Protein(accKvp.Value.Sequence.ToString(), accKvp.Key, isDecoy: accKvp.Value.IsDecoy); + foreach (var pending in accKvp.Value.Peptides) + { + var peptide = new PeptideWithSetModifications( + pending.FullSequence, mods, p: protein, digestionParams: digestionParams, + oneBasedStartResidueInProtein: pending.Start, + oneBasedEndResidueInProtein: pending.Start + pending.Length - 1, + missedCleavages: 0); + list.Add((peptide, pending.Fragments)); + } + } + result[dbKvp.Key] = list; + } + return result; + } + + private static double Pct(int n, int d) => d == 0 ? 0 : Math.Round(100.0 * n / d, 1); + + /// First index i where sorted[i] >= value (sorted ascending). + internal static int LowerBound(double[] sorted, double value) + { + int i = Array.BinarySearch(sorted, value); + return i < 0 ? ~i : i; + } + + /// + /// The .msl candidate pre-filter (extracted for testing): does any experimental scan match this library + /// entry on BOTH precursor mass (within the learned ppm window) AND retention time (within the calibrated + /// window of the entry's predicted RT)? An iRT of 0 (unpredicted) keeps the entry on mass alone. + /// reports whether the mass axis alone matched (for diagnostics). + /// + internal static bool IsCandidate(double precursorMz, int charge, float irt, in CandidatePriors priors, out bool massMatched) + { + massMatched = false; + double neutralMass = precursorMz.ToMass(charge); + double tolDa = neutralMass * priors.PrecursorTolPpm * 1e-6 + MassPreFilterMarginDa; + int a = LowerBound(priors.SortedScanMasses, neutralMass - tolDa); + double hi = neutralMass + tolDa; + if (a >= priors.SortedScanMasses.Length || priors.SortedScanMasses[a] > hi) + return false; // no mass match at all + massMatched = true; + + if (irt == 0f) + return true; // unpredicted RT -> keep on mass alone + double predRt = priors.RtSlope * irt + priors.RtIntercept; + for (int i = a; i < priors.SortedScanMasses.Length && priors.SortedScanMasses[i] <= hi; i++) + if (Math.Abs(priors.ScanRetentionTimes[i] - predRt) <= priors.RtWindowMin) + return true; + return false; + } + + /// + /// Splits a merged-index accession stamped "db|accession" into its (dbTag, accession) parts. An entry + /// with no bar falls back to dbTag "UNKNOWN"; an empty accession becomes "<dbTag>_UNKNOWN". + /// + internal static (string dbTag, string accession) ParseDbTagAndAccession(string tagged) + { + tagged ??= ""; + int bar = tagged.IndexOf('|'); + string dbTag = bar > 0 ? tagged.Substring(0, bar) : "UNKNOWN"; + string accession = bar >= 0 ? tagged.Substring(bar + 1) : tagged; + if (string.IsNullOrEmpty(accession)) + accession = dbTag + "_UNKNOWN"; + return (dbTag, accession); + } + + private static List BuildProducts(List ions) + { + var products = new List(ions.Count); + foreach (var ion in ions) + { + FragmentationTerminus terminus = + TerminusSpecificProductTypes.ProductTypeToFragmentationTerminus.TryGetValue( + ion.ProductType, out var t) ? t : FragmentationTerminus.None; + products.Add(new Product( + ion.ProductType, terminus, ion.Mz.ToMass(ion.Charge), ion.FragmentNumber, + ion.ResiduePosition, ion.NeutralLoss, ion.SecondaryProductType, ion.SecondaryFragmentNumber)); + } + return products; + } + } +} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs new file mode 100644 index 000000000..3c2e6c078 --- /dev/null +++ b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs @@ -0,0 +1,80 @@ +using System; +using Chemistry; +using MzLibUtil; + +namespace EngineLayer.ParallelSearch.Scoring +{ + /// + /// CPU implementation of and the correctness baseline for the + /// GPU path. Reproduces the original TransientClassicSearchEngine matching exactly: for each + /// theoretical fragment, find the closest experimental envelope (binary search), accept it if + /// it is within the product tolerance and its charge does not exceed the precursor charge. + /// + public sealed class CpuSpectralScorer : ISpectralScorer + { + private readonly SpectralScoringData _data; + private readonly Tolerance _productTolerance; + private FragmentMatch[] _matchBuffer = new FragmentMatch[64]; + + public CpuSpectralScorer(SpectralScoringData data, Tolerance productTolerance) + { + _data = data ?? throw new ArgumentNullException(nameof(data)); + _productTolerance = productTolerance ?? throw new ArgumentNullException(nameof(productTolerance)); + } + + public string BackendDescription => "CPU (binary search per fragment)"; + + public void ScoreBatch(ScoringBatch batch, IScoringSink sink) + { + if (batch == null) throw new ArgumentNullException(nameof(batch)); + if (sink == null) throw new ArgumentNullException(nameof(sink)); + + for (int w = 0; w < batch.WorkItemCount; w++) + { + int slot = batch.WorkPeptideSlot[w]; + int scanIndex = batch.WorkScanIndex[w]; + + int fragStart = batch.PeptideFragmentOffsets[slot]; + int fragEnd = batch.PeptideFragmentOffsets[slot + 1]; + int fragCount = fragEnd - fragStart; + + if (_matchBuffer.Length < fragCount) + _matchBuffer = new FragmentMatch[Math.Max(fragCount, _matchBuffer.Length * 2)]; + + int precursorCharge = _data.ScanPrecursorCharges[scanIndex]; + int matchCount = 0; + + for (int local = 0; local < fragCount; local++) + { + double theoreticalMass = batch.FragmentNeutralMasses[fragStart + local]; + + // unknown fragment mass; rare, for sequences with unknown amino acids + if (double.IsNaN(theoreticalMass)) + continue; + + int idx = _data.GetClosestFragmentIndex(scanIndex, theoreticalMass); + if (idx < 0) + continue; + + double experimentalMass = _data.FragmentMonoMasses[idx]; + int experimentalCharge = _data.FragmentCharges[idx]; + + if (_productTolerance.Within(experimentalMass, theoreticalMass) + && Math.Abs(experimentalCharge) <= Math.Abs(precursorCharge)) + { + _matchBuffer[matchCount++] = new FragmentMatch( + local, + experimentalMass.ToMz(experimentalCharge), + _data.FragmentIntensities[idx], + experimentalCharge); + } + } + + if (matchCount > 0) + sink.AcceptWorkItem(w, _matchBuffer, matchCount); + } + } + + public void Dispose() { /* nothing to release on the CPU path */ } + } +} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs new file mode 100644 index 000000000..521a1e5c4 --- /dev/null +++ b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs @@ -0,0 +1,95 @@ +using System; + +namespace EngineLayer.ParallelSearch.Scoring +{ + /// + /// Abstraction over the hot inner step of the transient classic search: for each + /// (peptide, candidate-scan) work item, decide which theoretical fragments match an + /// experimental fragment within tolerance and the charge filter. + /// + /// Deliberately narrow: the scorer does ONLY the matching (the expensive binary-search + + /// tolerance work that is identical across all 30k databases). Score computation, the + /// ScoreCutoff gate, MatchedFragmentIon construction and PSM updates all stay in the + /// engine. + /// + /// Implementation: + /// CpuSpectralScorer — binary search per fragment on the CPU. + /// + /// One instance per thread (not safe for concurrent ScoreBatch calls on the same instance). + /// + public interface ISpectralScorer : IDisposable + { + /// + /// Match every work item in , reporting each item's matched + /// fragments to . Items with no matches need not be reported. + /// + void ScoreBatch(ScoringBatch batch, IScoringSink sink); + + /// Human-readable backend description (for startup logging). + string BackendDescription { get; } + } + + /// One matched theoretical-vs-experimental fragment pairing. + public readonly struct FragmentMatch + { + /// Index of the theoretical product within its peptide's product list. + public readonly int LocalProductIndex; + /// Experimental m/z = monoisotopic mass converted at the envelope charge. + public readonly double ExperimentalMz; + /// Representative (first-peak) intensity of the matched experimental envelope. + public readonly double ExperimentalIntensity; + /// Charge of the matched experimental envelope. + public readonly int ExperimentalCharge; + + public FragmentMatch(int localProductIndex, double experimentalMz, double experimentalIntensity, int experimentalCharge) + { + LocalProductIndex = localProductIndex; + ExperimentalMz = experimentalMz; + ExperimentalIntensity = experimentalIntensity; + ExperimentalCharge = experimentalCharge; + } + } + + /// + /// Receives match results one work item at a time. The array is a + /// reusable buffer owned by the scorer — the sink must consume entries [0, matchCount) during + /// the call and not retain the array. + /// + public interface IScoringSink + { + void AcceptWorkItem(int workIndex, FragmentMatch[] matches, int matchCount); + } + + /// + /// A batch of (peptide, candidate-scan) work items to match, in struct-of-arrays form. + /// Theoretical fragment masses are concatenated per peptide slot (CSR via + /// ); NaN masses are kept in place so + /// aligns with the peptide's product list (the + /// scorer skips NaN, which never matches). + /// + public sealed class ScoringBatch + { + /// Number of work items. + public int WorkItemCount; + + /// Peptide slot for each work item (index into PeptideFragmentOffsets). + public int[] WorkPeptideSlot = Array.Empty(); + /// Candidate scan index for each work item (index into the SpectralScoringData). + public int[] WorkScanIndex = Array.Empty(); + + /// Number of distinct peptide slots in this batch. + public int PeptideCount; + /// CSR offsets into ; length = PeptideCount + 1. + public int[] PeptideFragmentOffsets = new int[1]; + /// Concatenated theoretical product neutral masses for all peptide slots. + public double[] FragmentNeutralMasses = Array.Empty(); + + /// Reset counters for reuse without reallocating the backing arrays. + public void Clear() + { + WorkItemCount = 0; + PeptideCount = 0; + PeptideFragmentOffsets[0] = 0; + } + } +} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs new file mode 100644 index 000000000..f80f56112 --- /dev/null +++ b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs @@ -0,0 +1,26 @@ +using MzLibUtil; + +namespace EngineLayer.ParallelSearch.Scoring +{ + /// + /// Hands out an to each search partition (thread). + /// Each thread gets its own lightweight scorer (per-thread match buffer, no shared state). + /// + public interface ISpectralScorerProvider : System.IDisposable + { + /// Scorer for the calling partition. Do NOT dispose the returned scorer; the provider owns it. + ISpectralScorer GetScorer(); + string BackendDescription { get; } + } + + /// CPU provider: a fresh (cheap) CpuSpectralScorer per partition. + public sealed class CpuScorerProvider : ISpectralScorerProvider + { + private readonly SpectralScoringData _data; + private readonly Tolerance _tolerance; + public CpuScorerProvider(SpectralScoringData data, Tolerance tolerance) { _data = data; _tolerance = tolerance; } + public ISpectralScorer GetScorer() => new CpuSpectralScorer(_data, _tolerance); + public string BackendDescription => "CPU (binary search per fragment)"; + public void Dispose() { } + } +} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs new file mode 100644 index 000000000..77061dc5c --- /dev/null +++ b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs @@ -0,0 +1,147 @@ +using System; +using System.Linq; +using Chemistry; +using MassSpectrometry; + +namespace EngineLayer.ParallelSearch.Scoring +{ + /// + /// Read-only, flattened (struct-of-arrays) view of the experimental MS2 spectra used by + /// the transient parallel search. Built ONCE from the shared, precursor-mass-sorted + /// array and reused, read-only, across every + /// transient database search. + /// + /// This is the data the GPU keeps resident in VRAM: experimental fragment masses are + /// concatenated into one sorted-per-scan array (CSR layout via ), + /// with parallel charge/intensity arrays and per-scan TIC / precursor charge. The CPU scorer + /// reconstructs s from the same data; the GPU scorer reads only + /// the flat arrays. + /// + /// Per-scan fragment slices are sorted ascending by monoisotopic mass (the + /// constructor guarantees ExperimentalFragments is sorted), + /// which is what makes the per-fragment closest-mass lookup a binary search. + /// + public sealed class SpectralScoringData + { + /// The original scans, kept for CPU-side MatchedFragmentIon reconstruction. + public Ms2ScanWithSpecificMass[] Scans { get; } + + /// Precursor masses, ascending (mirrors ClassicSearchEngine.MyScanPrecursorMasses). + public double[] ScanPrecursorMasses { get; } + + /// CSR offsets into the flat fragment arrays; length = Scans.Length + 1. + public int[] ScanFragmentOffsets { get; } + + /// Experimental fragment monoisotopic masses, concatenated; each scan's slice is sorted ascending. + public double[] FragmentMonoMasses { get; } + + /// Charge of each experimental fragment envelope (parallel to FragmentMonoMasses). + public int[] FragmentCharges { get; } + + /// Representative intensity (first peak) of each experimental fragment envelope. + public double[] FragmentIntensities { get; } + + /// Total ion current per scan (scoring normalization). + public double[] ScanTotalIonCurrents { get; } + + /// Precursor charge per scan (fragment-charge acceptance filter). + public int[] ScanPrecursorCharges { get; } + + public int ScanCount => Scans.Length; + + private SpectralScoringData( + Ms2ScanWithSpecificMass[] scans, double[] scanPrecursorMasses, int[] scanFragmentOffsets, + double[] fragmentMonoMasses, int[] fragmentCharges, double[] fragmentIntensities, + double[] scanTotalIonCurrents, int[] scanPrecursorCharges) + { + Scans = scans; + ScanPrecursorMasses = scanPrecursorMasses; + ScanFragmentOffsets = scanFragmentOffsets; + FragmentMonoMasses = fragmentMonoMasses; + FragmentCharges = fragmentCharges; + FragmentIntensities = fragmentIntensities; + ScanTotalIonCurrents = scanTotalIonCurrents; + ScanPrecursorCharges = scanPrecursorCharges; + } + + /// + /// Flatten the (already precursor-mass-sorted) scans into struct-of-arrays form. + /// One-time cost amortized over all transient database searches. + /// + public static SpectralScoringData Build(Ms2ScanWithSpecificMass[] sortedScans) + { + if (sortedScans == null) throw new ArgumentNullException(nameof(sortedScans)); + + int scanCount = sortedScans.Length; + var precursorMasses = new double[scanCount]; + var offsets = new int[scanCount + 1]; + var tics = new double[scanCount]; + var precursorCharges = new int[scanCount]; + + int totalFragments = 0; + for (int i = 0; i < scanCount; i++) + { + var frags = sortedScans[i].ExperimentalFragments; + totalFragments += frags?.Length ?? 0; + } + + var monoMasses = new double[totalFragments]; + var charges = new int[totalFragments]; + var intensities = new double[totalFragments]; + + int write = 0; + for (int i = 0; i < scanCount; i++) + { + var scan = sortedScans[i]; + offsets[i] = write; + precursorMasses[i] = scan.PrecursorMass; + tics[i] = scan.TotalIonCurrent; + precursorCharges[i] = scan.PrecursorCharge; + + var frags = scan.ExperimentalFragments; + if (frags != null) + { + for (int f = 0; f < frags.Length; f++) + { + var env = frags[f]; + monoMasses[write] = env.MonoisotopicMass; + charges[write] = env.Charge; + // Mirrors the existing engine: closestExperimentalMass.Peaks.First().intensity + intensities[write] = env.Peaks.First().intensity; + write++; + } + } + } + offsets[scanCount] = write; + + return new SpectralScoringData(sortedScans, precursorMasses, offsets, + monoMasses, charges, intensities, tics, precursorCharges); + } + + /// + /// Index of the experimental fragment whose mass is closest to + /// within scan 's slice, or -1 if the scan has no fragments. + /// Reproduces Ms2ScanWithSpecificMass.GetClosestFragmentMass over the flattened slice. + /// + public int GetClosestFragmentIndex(int scanIndex, double mass) + { + int lo = ScanFragmentOffsets[scanIndex]; + int hi = ScanFragmentOffsets[scanIndex + 1]; + int len = hi - lo; + if (len == 0) + return -1; + + int index = Array.BinarySearch(FragmentMonoMasses, lo, len, mass); + if (index >= 0) + return index; + + index = ~index; + + if (index == hi) + return index - 1; + if (index == lo || mass - FragmentMonoMasses[index - 1] > FragmentMonoMasses[index] - mass) + return index; + return index - 1; + } + } +} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs b/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs index 7e8a6bb1a..c0a91b0be 100644 --- a/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs +++ b/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs @@ -2,11 +2,13 @@ using System.Collections.Concurrent; using System.Collections.Generic; using System.Linq; +using System.Runtime.CompilerServices; using System.Threading; using System.Threading.Tasks; using Chemistry; using EngineLayer.ClassicSearch; using EngineLayer.FdrAnalysis; +using EngineLayer.ParallelSearch.Scoring; using EngineLayer.Util; using MzLibUtil; using Omics; @@ -29,15 +31,20 @@ public class TransientClassicSearchEngine : ClassicSearchEngine private bool _singleThreadMode; private readonly ConcurrentDictionary UpdatedIndexes = new(); private readonly bool _copyOnWriteEnabled; + // When supplied (e.g. from a .msl spectral library), the search iterates these peptides + // directly instead of digesting Proteins, and matches each peptide's STORED (float) fragments + // instead of re-fragmenting — skipping both digestion and fragmentation. + private readonly List<(IBioPolymerWithSetMods Peptide, List Fragments)> _precomputedPeptides; public TransientClassicSearchEngine(SpectralMatch[] globalPsms, Ms2ScanWithSpecificMass[] arrayOfSortedMS2Scans, - List variableModifications, List fixedModifications, + List variableModifications, List fixedModifications, List proteinList, MassDiffAcceptor searchMode, CommonParameters commonParameters, List<(string FileName, CommonParameters Parameters)> fileSpecificParameters, List nestedIds, - bool copyOnWriteEnabled = false) + bool copyOnWriteEnabled = false, List<(IBioPolymerWithSetMods Peptide, List Fragments)> precomputedPeptides = null) : base(globalPsms, arrayOfSortedMS2Scans, variableModifications, fixedModifications, null, null, null, proteinList, searchMode, commonParameters, fileSpecificParameters, null, nestedIds, false) { UpdatedIndexes = new ConcurrentDictionary(); _singleThreadMode = CommonParameters.MaxThreadsToUsePerFile <= 1; _copyOnWriteEnabled = copyOnWriteEnabled; + _precomputedPeptides = precomputedPeptides; if (!_singleThreadMode) { @@ -50,6 +57,13 @@ public TransientClassicSearchEngine(SpectralMatch[] globalPsms, Ms2ScanWithSpeci } } + // The flattened experimental spectra depend only on the (shared) scan array, so cache by + // its reference: every transient engine over the same file set reuses one build. + private static readonly ConditionalWeakTable _scoringDataCache = new(); + + private static SpectralScoringData GetOrBuildScoringData(Ms2ScanWithSpecificMass[] sortedScans) + => _scoringDataCache.GetValue(sortedScans, SpectralScoringData.Build); + protected override MetaMorpheusEngineResults RunSpecific() { double proteinsSearched = 0; @@ -58,100 +72,115 @@ protected override MetaMorpheusEngineResults RunSpecific() Status("Performing classic search..."); - if (Proteins.Any()) + bool usePrecomputedPeptides = _precomputedPeptides != null && _precomputedPeptides.Count > 0; + if (Proteins.Any() || usePrecomputedPeptides) { + // Experimental spectra are identical and read-only across every transient database + // search, so flatten them once (struct-of-arrays) and reuse across all peptides/threads. + var scoringData = GetOrBuildScoringData(ArrayOfSortedMS2Scans); + using var scorerProvider = new CpuScorerProvider( + scoringData, CommonParameters.ProductMassTolerance); + Status("ParallelSearch scoring backend: " + scorerProvider.BackendDescription); + + // Shared per-peptide work: fragment (in double), queue candidate scans, flush. The + // scorer does ONLY the fragment matching (the GPU-able hot step); the accumulator + // builds MatchedFragmentIons, applies the score cutoff, scores, and updates PSMs. + // A peptide's work items are queued contiguously so per-scan candidate ordering is + // preserved (bit-identical results). + // precomputedFragments != null (from a .msl library) skips the Fragment() call and matches + // the library's stored fragments directly — the fragmentation-time savings. + void ProcessOnePeptide(IBioPolymerWithSetMods peptide, TransientScoringBatch batch, + ISpectralScorer scorer, List peptideTheorProducts, List precomputedFragments = null) + { + Interlocked.Increment(ref peptideCounter); + + List products; + if (precomputedFragments != null) + { + products = precomputedFragments; + } + else + { + peptideTheorProducts.Clear(); + peptide.Fragment(CommonParameters.DissociationType, CommonParameters.DigestionParams.FragmentationTerminus, peptideTheorProducts, CommonParameters.FragmentationParameters); + products = peptideTheorProducts; + } + + int slot = -1; + foreach (ScanWithIndexAndNotchInfo scan in GetAcceptableScans(peptide.MonoisotopicMass, SearchMode)) + { + if (slot < 0) + slot = batch.BeginPeptide(peptide, products); + batch.AddWorkItem(slot, scan.ScanIndex, scan.Notch); + } + + if (batch.ShouldFlush) + batch.Flush(scorer); + } + + // GPU = one shared thread-safe scorer (resident spectra); CPU = a cheap per-partition + // scorer. The provider owns the scorer's lifetime. Action processProteinRange = (start, end) => { + var scorer = scorerProvider.GetScorer(); + var batch = new TransientScoringBatch(this, scoringData); List peptideTheorProducts = new(); - HashSet matchedFragmentIons = new(); - Tolerance productTolerance = CommonParameters.ProductMassTolerance; + for (int i = start; i < end; i++) { - // Stop loop if canceled - if (GlobalVariables.StopLoops) { return; } + if (GlobalVariables.StopLoops) { batch.Flush(scorer); return; } - // digest each protein into peptides and search for each peptide in all spectra within precursor mass tolerance + // digest each protein into peptides and search each peptide in all spectra within precursor mass tolerance foreach (var specificBioPolymer in Proteins[i].Digest(CommonParameters.DigestionParams, FixedModifications, VariableModifications)) - { - Interlocked.Increment(ref peptideCounter); - - peptideTheorProducts.Clear(); - specificBioPolymer.Fragment(CommonParameters.DissociationType, CommonParameters.DigestionParams.FragmentationTerminus, peptideTheorProducts, CommonParameters.FragmentationParameters); - - // score each scan that has an acceptable precursor mass - foreach (ScanWithIndexAndNotchInfo scan in GetAcceptableScans(specificBioPolymer.MonoisotopicMass, SearchMode)) - { - matchedFragmentIons.Clear(); - Ms2ScanWithSpecificMass theScan = ArrayOfSortedMS2Scans[scan.ScanIndex]; - int precursorCharge = theScan.PrecursorCharge; - - // Match Fragment Ions - foreach (var product in peptideTheorProducts) - { - // unknown fragment mass; this only happens rarely for sequences with unknown amino acids - if (double.IsNaN(product.NeutralMass)) - { - continue; - } - - // get the closest peak in the spectrum to the theoretical peak - var closestExperimentalMass = theScan.GetClosestExperimentalIsotopicEnvelope(product.NeutralMass); - - // is the mass error acceptable? - if (closestExperimentalMass != null - && productTolerance.Within(closestExperimentalMass.MonoisotopicMass, product.NeutralMass) - && Math.Abs(closestExperimentalMass.Charge) <= Math.Abs(precursorCharge))//TODO apply this filter before picking the envelope - { - matchedFragmentIons.Add(new MatchedFragmentIon(product, closestExperimentalMass.MonoisotopicMass.ToMz(closestExperimentalMass.Charge), - closestExperimentalMass.Peaks.First().intensity, closestExperimentalMass.Charge)); - } - } - - if (matchedFragmentIons.Count < CommonParameters.ScoreCutoff) - continue; - - // Score the peptide-spectrum match - double tic = theScan.TotalIonCurrent; - double score = 0; - foreach (var ion in matchedFragmentIons) - { - if (ion.NeutralTheoreticalProduct.ProductType != ProductType.D) - { - score += 1 + ion.Intensity / tic; - } - } - - var matchedIons = matchedFragmentIons.ToList(); // materialize before passing to another thread - AddPeptideCandidateToPsm(scan, score, specificBioPolymer, matchedIons); - } - } + ProcessOnePeptide(specificBioPolymer, batch, scorer, peptideTheorProducts); // report search progress (proteins searched so far out of total proteins in database) proteinsSearched++; var percentProgress = (int)((proteinsSearched / Proteins.Count) * 100); - if (percentProgress > oldPercentProgress) { oldPercentProgress = percentProgress; ReportProgress(new ProgressEventArgs(percentProgress, "Performing classic search... ", NestedIds)); } } + + batch.Flush(scorer); }; + // Library (.msl) peptide source: iterate precomputed peptides directly, no digestion. + Action processPeptideRange = (start, end) => + { + var scorer = scorerProvider.GetScorer(); + var batch = new TransientScoringBatch(this, scoringData); + List peptideTheorProducts = new(); + + for (int i = start; i < end; i++) + { + if (GlobalVariables.StopLoops) { batch.Flush(scorer); return; } + var (peptide, fragments) = _precomputedPeptides[i]; + ProcessOnePeptide(peptide, batch, scorer, peptideTheorProducts, fragments); + } + + batch.Flush(scorer); + }; + + int itemCount = usePrecomputedPeptides ? _precomputedPeptides.Count : Proteins.Count; + Action processRange = usePrecomputedPeptides ? processPeptideRange : processProteinRange; + if (_singleThreadMode) { - processProteinRange(0, Proteins.Count); + processRange(0, itemCount); } else { - var proteinPartioner = Partitioner.Create(0, Proteins.Count); + var partitioner = Partitioner.Create(0, itemCount); Parallel.ForEach( - proteinPartioner, + partitioner, new ParallelOptions { MaxDegreeOfParallelism = CommonParameters.MaxThreadsToUsePerFile }, (range, loopState) => { - processProteinRange(range.Item1, range.Item2); + processRange(range.Item1, range.Item2); }); } } @@ -319,6 +348,167 @@ private int GetFirstScanWithMassOverOrEqual(double minimum) // index of the first element that is larger than value return index; } + + /// + /// Per-thread accumulator that batches (peptide, candidate-scan) work for the + /// and, on flush, turns each item's matched fragments into a + /// PSM update — reproducing the original per-scan logic (HashSet dedup, ScoreCutoff gate, + /// (1 + intensity/TIC) scoring) exactly, just deferred to flush time. Work items are queued + /// in digestion order with each peptide's items contiguous, so PSM update order (and thus + /// tie resolution) matches the unbatched engine. + /// + private sealed class TransientScoringBatch : IScoringSink + { + private const int WorkItemFlushThreshold = 16384; + private const int PeptideFlushThreshold = 4096; + + private readonly TransientClassicSearchEngine _engine; + private readonly SpectralScoringData _data; + private readonly ScoringBatch _batch = new(); + // Reused across work items / peptides to keep the batched path allocation-free in the + // hot loop (the unbatched engine reused one cleared HashSet and never copied products). + private readonly HashSet _matchedIonScratch = new(); + private IBioPolymerWithSetMods[] _slotPeptides = new IBioPolymerWithSetMods[1024]; + private int[] _slotProductOffset = new int[1025]; + private Product[] _productPool = new Product[4096]; + private int[] _workNotch = new int[1024]; + private int _fragmentWrite; + private int _productWrite; + + public TransientScoringBatch(TransientClassicSearchEngine engine, SpectralScoringData data) + { + _engine = engine; + _data = data; + } + + public bool ShouldFlush => + _batch.WorkItemCount >= WorkItemFlushThreshold || _batch.PeptideCount >= PeptideFlushThreshold; + + /// Register a peptide and its theoretical products; returns its batch slot. + public int BeginPeptide(IBioPolymerWithSetMods peptide, List products) + { + int slot = _batch.PeptideCount; + int count = products.Count; + + EnsureSlotCapacity(slot + 1); + _slotPeptides[slot] = peptide; + _slotProductOffset[slot] = _productWrite; + + // Snapshot the products (the caller reuses/clears its list) into a pooled array, + // and mirror their neutral masses into the flat batch buffer for the scorer. + EnsureProductPoolCapacity(_productWrite + count); + EnsureFragmentCapacity(_fragmentWrite + count); + for (int k = 0; k < count; k++) + { + Product p = products[k]; + _productPool[_productWrite + k] = p; + _batch.FragmentNeutralMasses[_fragmentWrite + k] = p.NeutralMass; + } + _productWrite += count; + _fragmentWrite += count; + + EnsureOffsetCapacity(slot + 2); + _batch.PeptideFragmentOffsets[slot + 1] = _fragmentWrite; + _slotProductOffset[slot + 1] = _productWrite; + _batch.PeptideCount = slot + 1; + return slot; + } + + public void AddWorkItem(int slot, int scanIndex, int notch) + { + int w = _batch.WorkItemCount; + EnsureWorkCapacity(w + 1); + _batch.WorkPeptideSlot[w] = slot; + _batch.WorkScanIndex[w] = scanIndex; + _workNotch[w] = notch; + _batch.WorkItemCount = w + 1; + } + + public void Flush(ISpectralScorer scorer) + { + if (_batch.WorkItemCount > 0) + scorer.ScoreBatch(_batch, this); + Reset(); + } + + private void Reset() + { + _batch.Clear(); + _fragmentWrite = 0; + _productWrite = 0; + } + + void IScoringSink.AcceptWorkItem(int workIndex, FragmentMatch[] matches, int matchCount) + { + int slot = _batch.WorkPeptideSlot[workIndex]; + int scanIndex = _batch.WorkScanIndex[workIndex]; + int notch = _workNotch[workIndex]; + int productBase = _slotProductOffset[slot]; + + HashSet matchedFragmentIons = _matchedIonScratch; + matchedFragmentIons.Clear(); + for (int k = 0; k < matchCount; k++) + { + FragmentMatch m = matches[k]; + matchedFragmentIons.Add(new MatchedFragmentIon( + _productPool[productBase + m.LocalProductIndex], m.ExperimentalMz, m.ExperimentalIntensity, m.ExperimentalCharge)); + } + + if (matchedFragmentIons.Count < _engine.CommonParameters.ScoreCutoff) + return; + + double tic = _data.ScanTotalIonCurrents[scanIndex]; + double score = 0; + foreach (var ion in matchedFragmentIons) + { + if (ion.NeutralTheoreticalProduct.ProductType != ProductType.D) + score += 1 + ion.Intensity / tic; + } + + var matchedIons = matchedFragmentIons.ToList(); + _engine.AddPeptideCandidateToPsm( + new ScanWithIndexAndNotchInfo(notch, scanIndex), score, _slotPeptides[slot], matchedIons); + } + + private void EnsureSlotCapacity(int neededSlots) + { + if (_slotPeptides.Length < neededSlots) + { + int n = Math.Max(neededSlots, _slotPeptides.Length * 2); + Array.Resize(ref _slotPeptides, n); + Array.Resize(ref _slotProductOffset, n + 1); + } + } + + private void EnsureProductPoolCapacity(int needed) + { + if (_productPool.Length < needed) + Array.Resize(ref _productPool, Math.Max(needed, Math.Max(4096, _productPool.Length * 2))); + } + + private void EnsureFragmentCapacity(int needed) + { + if (_batch.FragmentNeutralMasses.Length < needed) + Array.Resize(ref _batch.FragmentNeutralMasses, Math.Max(needed, Math.Max(1024, _batch.FragmentNeutralMasses.Length * 2))); + } + + private void EnsureOffsetCapacity(int needed) + { + if (_batch.PeptideFragmentOffsets.Length < needed) + Array.Resize(ref _batch.PeptideFragmentOffsets, Math.Max(needed, _batch.PeptideFragmentOffsets.Length * 2)); + } + + private void EnsureWorkCapacity(int needed) + { + if (_batch.WorkPeptideSlot.Length < needed) + { + int n = Math.Max(needed, Math.Max(1024, _batch.WorkPeptideSlot.Length * 2)); + Array.Resize(ref _batch.WorkPeptideSlot, n); + Array.Resize(ref _batch.WorkScanIndex, n); + Array.Resize(ref _workNotch, n); + } + } + } } public class TransientSearchEngineResults( diff --git a/MetaMorpheus/GUI/GUI.csproj b/MetaMorpheus/GUI/GUI.csproj index f1b354242..658d8ec0e 100644 --- a/MetaMorpheus/GUI/GUI.csproj +++ b/MetaMorpheus/GUI/GUI.csproj @@ -1,4 +1,4 @@ - + WinExe @@ -63,7 +63,7 @@ - + diff --git a/MetaMorpheus/GuiFunctions/GuiFunctions.csproj b/MetaMorpheus/GuiFunctions/GuiFunctions.csproj index 9c2aaa847..7dd533525 100644 --- a/MetaMorpheus/GuiFunctions/GuiFunctions.csproj +++ b/MetaMorpheus/GuiFunctions/GuiFunctions.csproj @@ -1,4 +1,4 @@ - + net8.0-windows @@ -16,7 +16,7 @@ - + diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs b/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs index c86caf07d..3731d0bcd 100644 --- a/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs +++ b/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs @@ -32,6 +32,16 @@ public StatisticalTestResultFile(double alpha = 0.05) : base() public override void LoadResults() { + // Writer-only usage (parameterless ctor with in-memory Results) has no file to read. + // When the in-memory result list is empty the base lazy getter still calls LoadResults; + // guard so that does not throw on an empty/missing FilePath (e.g. a run where no + // statistical results were produced). Set via the setter to avoid getter recursion. + if (string.IsNullOrEmpty(FilePath) || !File.Exists(FilePath)) + { + Results = new List(); + return; + } + using var reader = new StreamReader(FilePath); using var csv = new CsvReader(reader, new CsvConfiguration(CultureInfo.InvariantCulture) { diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs b/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs index f45121d2e..53b7c5736 100644 --- a/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs +++ b/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs @@ -1,4 +1,5 @@ #nullable enable +using System.Collections.Generic; using System.Globalization; using System.IO; using System.Linq; @@ -27,6 +28,14 @@ public TestSummaryResultFile() : base() { } public override void LoadResults() { + // Writer-only usage (parameterless ctor with in-memory Results) has no file to read; + // guard against the base lazy getter calling this on an empty/missing FilePath. + if (string.IsNullOrEmpty(FilePath) || !File.Exists(FilePath)) + { + Results = new List(); + return; + } + using var csv = new CsvReader(new StreamReader(FilePath), CsvConfiguration); Results = csv.GetRecords().ToList(); } diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs index d9ef752aa..18fa371b3 100644 --- a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs +++ b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs @@ -9,12 +9,15 @@ using EngineLayer.Util; using FlashLFQ; using Nett; +using Chromatography.RetentionTimePrediction; using Omics; +using Omics.Fragmentation; using Omics.Modifications; using Omics.SpectrumMatch; using System; using System.Collections.Concurrent; using System.Collections.Generic; +using System.Diagnostics; using System.IO; using System.Linq; using System.Threading; @@ -38,6 +41,12 @@ public class ParallelSearchTask : SearchTask private int _dashboardInitialFinishedCount; private int _dashboardCompletedThisRun; + // Phase timing (P0 of GPU plan): reveals where wall-clock actually goes so GPU effort + // targets the real hot spot. Inner-loop split is accumulated across the parallel + // per-database loop via Interlocked; coarse phases are timed in RunSpecific. + private long _searchEngineTicks; // time inside TransientClassicSearchEngine.Run + private long _postAnalysisTicks; // time inside PerformPostSearchAnalysis + private const string DashboardPhaseInitializing = "Initializing"; private const string DashboardPhaseSearching = "Searching"; private const string DashboardPhaseStatisticalAnalysis = "Statistical analysis"; @@ -106,6 +115,14 @@ public override SearchParameters SearchParameters [TomlIgnore] public Ms2ScanWithSpecificMass[] AllSortedMs2Scans { get; private set; } = []; [TomlIgnore] private SpectralMatch[] BaseSearchPsms = null!; // PSMs from base database search + // Calibration LEARNED from the base search (S1), applied to the .msl candidate pre-filter (S3): + // precursor tolerance, and the iRT→observed-RT regression + window. Null until learned (and stays + // null when no transient database is a .msl library). See LearnMslCalibration. + [TomlIgnore] private MslCandidateCalibration? _mslCalibration; + // Scan precursor masses (ascending) + their RTs, cached once (identical across all .msl databases). + [TomlIgnore] private double[]? _sortedScanMasses; + [TomlIgnore] private double[]? _scanRetentionTimes; + // Optimization caches for FDR alignment and parsimony [TomlIgnore] private readonly PsmSpectralMatchFdrAlignmentService _psmFdrAlignmentService = new(); @@ -156,35 +173,125 @@ protected override MyTaskResults RunSpecific(string outputFolder, goto Finalization; } + // Phase timing (P0 of GPU plan) — coarse wall-clock per phase. + var swInit = Stopwatch.StartNew(); + // Initialize all necessary data structures including base search Initialize(taskId, dbFilenameList, currentRawFileList, fileSettingsList, outputFolder); InitializeCompletedDatabaseWriter(); + swInit.Stop(); Status($"Starting search of {TotalDatabases} transient databases...", taskId); ReportTaskDashboard(taskId, ParallelSearchDashboardUpdateKind.TaskStatus, DashboardPhaseSearching, $"Searching {TotalDatabases} transient databases..."); + // MERGED-INDEX mode: one .msl file holds many databases (entries tagged "db|accession"). The number + // of databases to search comes from inside the file, not the (count==1) file list, so parallelism + // must not be capped by the file count. + bool mergedMode = Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_MERGED") == "1" + && ParallelSearchParameters.TransientDatabases.Count == 1 + && ParallelSearchParameters.TransientDatabases[0].FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase); + // Determine optimal thread allocation int totalAvailableThreads = Environment.ProcessorCount; - int databaseParallelism = Math.Min(ParallelSearchParameters.MaxSearchesInParallel, - ParallelSearchParameters.TransientDatabases.Count); + int databaseParallelism = mergedMode + ? ParallelSearchParameters.MaxSearchesInParallel + : Math.Min(ParallelSearchParameters.MaxSearchesInParallel, + ParallelSearchParameters.TransientDatabases.Count); int threadsPerDatabase = Math.Max(1, totalAvailableThreads / databaseParallelism); CommonParameters.MaxThreadsToUsePerFile = threadsPerDatabase; - // Loop through each transient database - try + // S4: LOAD-AHEAD PRODUCER/CONSUMER. Database loading (now an I/O-bound .msl index query + candidate + // fragment fetch) runs ahead on producer threads and hands PREPARED databases to searcher threads + // through a bounded queue, so loading overlaps with the CPU-bound search instead of blocking it. + // The result WRITE channel was already pipelined; this pipelines the read side too. The bounded + // capacity caps how many loaded databases are held in memory at once. + var swTransientLoop = Stopwatch.StartNew(); + using (var loadedQueue = new System.Collections.Concurrent.BlockingCollection( + boundedCapacity: Math.Max(2, databaseParallelism))) { - Parallel.ForEach(ParallelSearchParameters.TransientDatabases, - new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism }, - transientDbPath => + // MERGED-INDEX mode (mergedMode computed above): a single .msl holds every database (each entry + // tagged "db|accession"). Load it ONCE and run the candidate filter ONCE, then emit one prepared + // database per source db-group. Each is searched INDEPENDENTLY by the consumer below (own + // base-PSM copy) — databases never compete; this only collapses 1000s of file opens into one + // shared in-memory index. + + // Producers: load databases concurrently, blocking on the bounded queue when it is full. + var producer = Task.Run(() => + { + try { - ProcessTransientDatabase(transientDbPath, outputFolder, taskId); - }); - } - finally - { - CompleteCompletedDatabaseWriter(); + if (mergedMode) + { + Status("Loading merged transient index...", taskId); + var grouped = MslPeptideReader.ReadCandidatesGroupedByDatabase( + ParallelSearchParameters.TransientDatabases[0].FilePath, BuildCandidatePriors()); + Status($"Merged index: {grouped.Count} databases with candidates.", taskId); + Parallel.ForEach(grouped, + new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism }, + kvp => + { + if (GlobalVariables.StopLoops) return; + var loaded = BuildLoadedFromCandidates(kvp.Key, kvp.Value, outputFolder, taskId); + if (loaded != null) loadedQueue.Add(loaded); + }); + } + else + { + Parallel.ForEach(ParallelSearchParameters.TransientDatabases, + new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism }, + transientDbPath => + { + if (GlobalVariables.StopLoops) return; + LoadedTransientDatabase? loaded; + try + { + loaded = LoadTransientDatabaseForPipeline(transientDbPath, outputFolder, taskId); + } + catch (Exception ex) + { + WriteTransientProcessError(transientDbPath, outputFolder, ex); + throw; + } + if (loaded != null) + loadedQueue.Add(loaded); + }); + } + } + finally + { + loadedQueue.CompleteAdding(); + } + }); + + try + { + // Consumers: search each loaded database as it becomes available. + Parallel.ForEach(loadedQueue.GetConsumingEnumerable(), + new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism }, + loaded => + { + if (GlobalVariables.StopLoops) return; + try + { + SearchLoadedTransientDatabase(loaded, taskId); + } + catch (Exception ex) + { + WriteTransientProcessError(loaded.TransientDb, outputFolder, ex); + throw; + } + }); + } + finally + { + CompleteCompletedDatabaseWriter(); + } + + producer.GetAwaiter().GetResult(); // surface any producer (load) exception } + swTransientLoop.Stop(); + LogPhaseTimingBreakdown(taskId, outputFolder, swInit.Elapsed, swTransientLoop.Elapsed, databaseParallelism); Finalization: @@ -237,7 +344,16 @@ protected override MyTaskResults RunSpecific(string outputFolder, Status("Writing Final Results...", taskId); ReportTaskDashboard(taskId, ParallelSearchDashboardUpdateKind.TaskStatus, DashboardPhaseWritingFinalResults, "Writing final results..."); - WriteFinalOutputs(outputFolder, taskId, currentRawFileList.Count); + try + { + WriteFinalOutputs(outputFolder, taskId, currentRawFileList.Count); + } + catch (Exception ex) + { + // Dump the full stack to a labeled file — the task runner only surfaces the message. + try { File.WriteAllText(Path.Combine(outputFolder, "WriteFinalOutputs_ERROR.txt"), ex.ToString()); } catch { } + throw; + } ReportTaskDashboard(taskId, ParallelSearchDashboardUpdateKind.TaskCompleted, DashboardPhaseCompleted, "Many search task complete!"); @@ -246,6 +362,44 @@ protected override MyTaskResults RunSpecific(string outputFolder, return MyTaskResults; } + /// + /// Logs where wall-clock time went so GPU effort targets the real hot spot. + /// The per-database search and post-analysis times are summed across all parallel + /// workers (CPU-seconds); dividing by the database parallelism approximates the + /// wall-clock each contributed to the transient loop. + /// + private void LogPhaseTimingBreakdown(string taskId, string outputFolder, TimeSpan init, TimeSpan transientLoop, int databaseParallelism) + { + double searchCpuSec = _searchEngineTicks / (double)Stopwatch.Frequency; + double postCpuSec = _postAnalysisTicks / (double)Stopwatch.Frequency; + int par = Math.Max(1, databaseParallelism); + + string summary = + "Phase timing — " + + $"Initialize (load + base search): {init.TotalSeconds:F1}s | " + + $"Transient loop wall-clock: {transientLoop.TotalSeconds:F1}s | " + + $"search engine: {searchCpuSec:F1} CPU-s (~{searchCpuSec / par:F1}s wall) | " + + $"post-search analysis: {postCpuSec:F1} CPU-s (~{postCpuSec / par:F1}s wall) | " + + $"db parallelism: {par}."; + + Status(summary, taskId); + + // Also persist to a flushed file so the profile survives any later-stage failure + // (the breakdown is the whole point of the run; don't let finalization bugs eat it). + try + { + string detail = + summary + Environment.NewLine + Environment.NewLine + + "Breakdown (transient loop only; Initialize is a one-time cost):" + Environment.NewLine + + $" search engine (TransientClassicSearchEngine): {searchCpuSec,10:F1} CPU-s" + Environment.NewLine + + $" post-search analysis (stats/collectors/FDR): {postCpuSec,10:F1} CPU-s" + Environment.NewLine + + $" search : post ratio = {searchCpuSec / Math.Max(1e-9, postCpuSec):F2} : 1" + Environment.NewLine + + $" db parallelism = {par}" + Environment.NewLine; + File.WriteAllText(Path.Combine(outputFolder, "PhaseTiming.txt"), detail); + } + catch { /* timing file is best-effort */ } + } + #region Initialization private void Initialize(string taskId, List dbFilenameList, @@ -317,6 +471,19 @@ private void Initialize(string taskId, List dbFilenameList, _peptideFdrAlignmentService.BuildBaselineCache(baselinePsms); BuildBaselinePeptideToScanIndexLookup(); + // S1: when any transient database is a .msl library, learn the precursor tolerance and the + // RT (iRT→observed) calibration from the confident base-search PSMs. These become the .msl + // candidate pre-filter priors (S3). Cheap relative to the search and done once. + if (ParallelSearchParameters.TransientDatabases.Any( + d => d.FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase))) + { + _mslCalibration = LearnMslCalibration(baselinePsms, taskId); + // Cache the (mass-sorted) scan mass + RT arrays ONCE — they are identical for every .msl + // database, so rebuilding them per database (1000s of times) would be pure waste. + _sortedScanMasses = Array.ConvertAll(AllSortedMs2Scans, s => s.PrecursorMass); + _scanRetentionTimes = Array.ConvertAll(AllSortedMs2Scans, s => s.RetentionTime); + } + if (SearchParameters.DoParsimony) { Status("Preparing baseline parsimony cache...", taskId); @@ -490,32 +657,99 @@ private TransientDatabaseResultsManager CreateResultsManager(string outputFolder #region Search - private void ProcessTransientDatabase(DbForTask transientDb, string outputFolder, string taskId) + /// + /// A transient database after the LOAD stage (producer): proteins/peptides ready, fragments fetched. + /// Carried through the bounded channel to a searcher (consumer) so loading overlaps with searching. + /// + private sealed class LoadedTransientDatabase + { + public DbForTask TransientDb = null!; + public string DbName = null!; + public string DbOutputFolder = null!; + public List NestedIds = null!; + public List TransientProteins = null!; + public List<(IBioPolymerWithSetMods Peptide, List Fragments)>? PrecomputedPeptides; + public HashSet TransientProteinAccessions = null!; + } + + /// Builds the .msl candidate pre-filter priors (scan masses + RTs, learned precursor/RT + /// calibration). The scan arrays are cached once (identical for every database). + private MslPeptideReader.CandidatePriors BuildCandidatePriors() + { + var sortedScanMasses = _sortedScanMasses ?? Array.ConvertAll(AllSortedMs2Scans, s => s.PrecursorMass); + var scanRetentionTimes = _scanRetentionTimes ?? Array.ConvertAll(AllSortedMs2Scans, s => s.RetentionTime); + var cal = _mslCalibration ?? new MslCandidateCalibration( + CommonParameters.PrecursorMassTolerance.Value, 1, 0, double.PositiveInfinity); + return new MslPeptideReader.CandidatePriors( + sortedScanMasses, scanRetentionTimes, + precursorTolPpm: cal.PrecursorTolPpm, rtSlope: cal.RtSlope, + rtIntercept: cal.RtIntercept, rtWindowMin: cal.RtWindowMin); + } + + /// + /// Builds a prepared database from a merged-index db-group's candidate peptides (no file load). The + /// returned database is searched INDEPENDENTLY by the consumer, identical to the per-file path. + /// + private LoadedTransientDatabase? BuildLoadedFromCandidates(string dbName, + List<(IBioPolymerWithSetMods Peptide, List Fragments)> candidates, string outputFolder, string taskId) { if (GlobalVariables.StopLoops) - return; + return null; + + List nestedIds = [taskId, dbName]; + bool shouldProcess = !_resultsManager!.HasCachedResults(dbName) || ParallelSearchParameters.OverwriteTransientSearchOutputs; + if (!shouldProcess) + { + ReportProgress(new(100, $"Skipping {dbName} - results already exist in cache", nestedIds)); + UpdateProgress(TotalDatabases, taskId); + return null; + } + + var transientProteins = candidates.Select(p => p.Peptide.Parent).Distinct().ToList(); + return new LoadedTransientDatabase + { + // Synthetic per-group identity (name only); no per-database file in merged mode. + TransientDb = new DbForTask(dbName, false), + DbName = dbName, + DbOutputFolder = Path.Combine(outputFolder, dbName), + NestedIds = nestedIds, + TransientProteins = transientProteins, + PrecomputedPeptides = candidates, + TransientProteinAccessions = new HashSet(transientProteins.Select(p => p.Accession)), + }; + } + + /// + /// PRODUCER stage: skip/overwrite handling + load the transient database (FASTA digest set, or the + /// .msl mass+RT-filtered candidate peptides with their fragments). Returns null when the database is + /// skipped (cached) or the run is stopping. I/O-bound; runs ahead of the searchers. + /// + private LoadedTransientDatabase? LoadTransientDatabaseForPipeline(DbForTask transientDb, string outputFolder, string taskId) + { + if (GlobalVariables.StopLoops) + return null; - string dbName = Path.GetFileNameWithoutExtension(transientDb.FilePath); - string dbOutputFolder = Path.Combine(outputFolder, dbName); - List nestedIds = [taskId, dbName]; + string dbName = Path.GetFileNameWithoutExtension(transientDb.FilePath); + string dbOutputFolder = Path.Combine(outputFolder, dbName); + List nestedIds = [taskId, dbName]; Status($"Processing {dbName}...", nestedIds); - // Check if we should skip or overwrite this database - bool shouldProcess = !_resultsManager!.HasCachedResults(dbName) || - ParallelSearchParameters.OverwriteTransientSearchOutputs; + // Check if we should skip or overwrite this database + bool shouldProcess = !_resultsManager!.HasCachedResults(dbName) || + ParallelSearchParameters.OverwriteTransientSearchOutputs; if (!shouldProcess) { ReportProgress(new(100, $"Skipping {dbName} - results already exist in cache", nestedIds)); UpdateProgress(TotalDatabases, taskId); - return; + return null; } ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseStarted, dbName, $"Processing {dbName}...", DashboardDatabaseProcessingProgress); - // Handle overwrite scenario + // Handle overwrite scenario if (ParallelSearchParameters.OverwriteTransientSearchOutputs && _resultsManager.HasCachedResults(dbName)) { Status($"Overwriting existing results for {dbName}...", nestedIds); @@ -527,55 +761,99 @@ private void ProcessTransientDatabase(DbForTask transientDb, string outputFolder } } - if (!Directory.Exists(dbOutputFolder)) - Directory.CreateDirectory(dbOutputFolder); + // NOTE: the output folder is created lazily by the writer, and ONLY for databases that produced + // transient PSMs — at 1000s of databases the vast majority match nothing, so creating a folder + + // header-only files for each was a large fraction of the runtime. See WriteCompletedDatabaseOutputsAsync. Status($"Loading transient database {dbName}...", nestedIds); ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName, $"Loading transient database {dbName}...", DashboardDatabaseLoadingProgress); - // Load transient database - var transientProteins = LoadTransientDatabase(transientDb, nestedIds, taskId); + // Load transient database. FASTA/XML -> proteins (digested during search); a .msl spectral + // library -> precomputed peptides paired with their stored (float) fragments (the search + // iterates these and matches the stored fragments directly, skipping digestion+fragmentation). + List transientProteins; + List<(IBioPolymerWithSetMods Peptide, List Fragments)>? precomputedPeptides = null; + bool isMslLibrary = transientDb.FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase); + if (isMslLibrary) + { + // S3: index-only candidate pre-filter on precursor mass AND Chronologer-iRT (learned from the + // base search, S1). See BuildCandidatePriors. + precomputedPeptides = MslPeptideReader.ReadPeptides(transientDb.FilePath, dbName, BuildCandidatePriors()); + // shared parent proteins (one per accession) back the accession filter and counts + transientProteins = precomputedPeptides.Select(p => p.Peptide.Parent).Distinct().ToList(); + } + else + { + transientProteins = LoadTransientDatabase(transientDb, nestedIds, taskId); + } - if (GlobalVariables.StopLoops) - { - return; - } + if (GlobalVariables.StopLoops) + return null; - // Create HashSet of transient bioPolymer accessions for later filtering - var transientProteinAccessions = new HashSet( - transientProteins.Select(p => p.Accession)); + return new LoadedTransientDatabase + { + TransientDb = transientDb, + DbName = dbName, + DbOutputFolder = dbOutputFolder, + NestedIds = nestedIds, + TransientProteins = transientProteins, + PrecomputedPeptides = precomputedPeptides, + // Create HashSet of transient bioPolymer accessions for later filtering + TransientProteinAccessions = new HashSet(transientProteins.Select(p => p.Accession)), + }; + } + + /// + /// CONSUMER stage: search the loaded database against the shared spectra and run post-analysis, then + /// hand the results to the (already-pipelined) write channel. CPU-bound; overlaps with loaders. + /// + private void SearchLoadedTransientDatabase(LoadedTransientDatabase loaded, string taskId) + { + if (GlobalVariables.StopLoops) + return; + + DbForTask transientDb = loaded.TransientDb; + string dbName = loaded.DbName; + string dbOutputFolder = loaded.DbOutputFolder; + List nestedIds = loaded.NestedIds; + List transientProteins = loaded.TransientProteins; + HashSet transientProteinAccessions = loaded.TransientProteinAccessions; Status($"Searching {dbName} ({transientProteins.Count} transient proteins)...", nestedIds); ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName, $"Searching {dbName} ({transientProteins.Count} transient proteins)...", DashboardDatabaseSearchStartProgress, DashboardDatabaseSearchStartProgress, DashboardDatabaseSearchEndProgress); - // Reuse baseline PSMs with copy-on-write in peptide/proteoform mode. - SpectralMatch[] psmArray = BaseSearchPsms.ToArray(); - PerformSearch(transientProteins, psmArray, nestedIds, out HashSet updatedPsmIndexes, out int transientPeptideCount, useCopyOnWrite: true); + // Reuse baseline PSMs with copy-on-write in peptide/proteoform mode. + SpectralMatch[] psmArray = BaseSearchPsms.ToArray(); + long searchStart = Stopwatch.GetTimestamp(); + PerformSearch(transientProteins, psmArray, nestedIds, out HashSet updatedPsmIndexes, out int transientPeptideCount, useCopyOnWrite: true, precomputedPeptides: loaded.PrecomputedPeptides); + Interlocked.Add(ref _searchEngineTicks, Stopwatch.GetTimestamp() - searchStart); Status($"Performing post-search analysis for {dbName}...", nestedIds); ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName, $"Performing post-search analysis for {dbName}...", DashboardDatabasePostSearchProgress); - int totalProteins = BaseBioPolymers.Count + transientProteins.Count; - - // Process database through unified manager (handles analysis + statistical caching) - var (analysisContext, dbResults) = PerformPostSearchAnalysis( - psmArray, - dbOutputFolder, - nestedIds, - dbName, - totalProteins, - transientProteinAccessions, - transientPeptideCount, - transientProteins, - transientDb, - updatedPsmIndexes - ).GetAwaiter().GetResult(); - - _completedDatabaseWriteChannel!.Writer.WriteAsync((analysisContext, dbResults)).AsTask().GetAwaiter().GetResult(); + int totalProteins = BaseBioPolymers.Count + transientProteins.Count; + + // Process database through unified manager (handles analysis + statistical caching) + long postAnalysisStart = Stopwatch.GetTimestamp(); + var (analysisContext, dbResults) = PerformPostSearchAnalysis( + psmArray, + dbOutputFolder, + nestedIds, + dbName, + totalProteins, + transientProteinAccessions, + transientPeptideCount, + transientProteins, + transientDb, + updatedPsmIndexes + ).GetAwaiter().GetResult(); + Interlocked.Add(ref _postAnalysisTicks, Stopwatch.GetTimestamp() - postAnalysisStart); + + _completedDatabaseWriteChannel!.Writer.WriteAsync((analysisContext, dbResults)).AsTask().GetAwaiter().GetResult(); // Cleanup transient proteins to free memory transientProteins.Clear(); @@ -583,21 +861,136 @@ private void ProcessTransientDatabase(DbForTask transientDb, string outputFolder ReportProgress(new(100, $"Finished analysis for {dbName}; queued output writing", nestedIds)); } + /// Writes a per-database <db>_PROCESS_ERROR.txt diagnostic; never throws. + private static void WriteTransientProcessError(DbForTask transientDb, string outputFolder, Exception ex) + { + try + { + string dbName = Path.GetFileNameWithoutExtension(transientDb.FilePath); + File.WriteAllText(Path.Combine(outputFolder, dbName + "_PROCESS_ERROR.txt"), ex.ToString()); + } + catch { } + } + + /// Calibration learned from the base search and applied to the .msl candidate pre-filter. + private readonly struct MslCandidateCalibration + { + public readonly double PrecursorTolPpm; // symmetric precursor tolerance (covers offset+spread) + public readonly double RtSlope, RtIntercept; // observedRT = RtSlope*iRT + RtIntercept + public readonly double RtWindowMin; // +/- observed-RT window (k * residual SD); +inf = no RT filter + public MslCandidateCalibration(double tolPpm, double slope, double intercept, double rtWindowMin) + { PrecursorTolPpm = tolPpm; RtSlope = slope; RtIntercept = intercept; RtWindowMin = rtWindowMin; } + } + + /// + /// S1 — learn the .msl candidate pre-filter priors from confident base-search PSMs: + /// • precursor tolerance = |mean| + 3·SD of the precursor mass error (clamped to the configured tol), + /// • RT calibration = OLS of observed RT on Chronologer iRT, with a ±2·residualSD window. + /// On any failure (too few PSMs, predictor/model unavailable) returns a SAFE fallback: the configured + /// precursor tolerance and NO RT filter (RtWindowMin = +inf) so nothing is lost. + /// + private MslCandidateCalibration LearnMslCalibration(List baselinePsms, string taskId) + { + double configuredTolPpm = CommonParameters.PrecursorMassTolerance.Value; + var fallback = new MslCandidateCalibration(configuredTolPpm, 1, 0, double.PositiveInfinity); + try + { + // Confident, unambiguous target PSMs. + var conf = baselinePsms.Where(p => p != null && !p.IsDecoy + && p.PsmFdrInfo != null && p.PsmFdrInfo.QValue < 0.01 + && p.BestMatchingBioPolymersWithSetMods.Count() == 1) + .ToList(); + if (conf.Count < 50) + { + Status($"S1: only {conf.Count} confident base PSMs — using safe fallback (no RT filter).", taskId); + return fallback; + } + + // Precursor mass error (ppm) and the (peptide, observedRT) pairs for the RT regression. + var ppmErrors = new List(conf.Count); + var obsByPeptide = new Dictionary(ReferenceEqualityComparer.Instance); + var peptides = new List(conf.Count); + foreach (var psm in conf) + { + var pep = psm.BestMatchingBioPolymersWithSetMods.First().SpecificBioPolymer; + double mass = pep.MonoisotopicMass; + if (mass <= 0) continue; + ppmErrors.Add((psm.ScanPrecursorMass - mass) / mass * 1e6); + if (pep is IRetentionPredictable rp && !obsByPeptide.ContainsKey(rp)) + { + obsByPeptide[rp] = psm.ScanRetentionTime; + peptides.Add(rp); + } + } + + double pMean = ppmErrors.Average(); + double pSd = Math.Sqrt(ppmErrors.Sum(e => (e - pMean) * (e - pMean)) / ppmErrors.Count); + // Trust the calibration: the real precursor accuracy (|mean|+3σ of the confident base PSMs) IS + // the tolerance, regardless of the looser configured/standard value. This is applied to BOTH the + // candidate pre-filter AND the transient search's acceptor (see PerformSearch), so they are + // consistent — peptides outside it are rejected as loose-tolerance false positives, not lost. + double tolPpm = Math.Max(2.0, Math.Abs(pMean) + 3.0 * pSd); + + // Chronologer iRT for the confident peptides, then OLS observedRT vs iRT. + using var predictor = RetentionTimePredictorFactory.Create(PredictorType.Chronologer); + var predRt = new List(peptides.Count); + var obsRt = new List(peptides.Count); + foreach (var r in predictor.PredictRetentionTimeEquivalents(peptides, + maxThreads: Math.Max(1, Environment.ProcessorCount - 2))) + { + if (r.PredictedValue is null) continue; + if (obsByPeptide.TryGetValue(r.Peptide, out double rt)) { predRt.Add(r.PredictedValue.Value); obsRt.Add(rt); } + } + if (predRt.Count < 50) + { + Status($"S1: only {predRt.Count} RT predictions — precursor tol {tolPpm:F1} ppm, no RT filter.", taskId); + return new MslCandidateCalibration(tolPpm, 1, 0, double.PositiveInfinity); + } + + int n = predRt.Count; + double mx = predRt.Average(), my = obsRt.Average(); + double sxy = 0, sxx = 0; + for (int i = 0; i < n; i++) { double dx = predRt[i] - mx; sxy += dx * (obsRt[i] - my); sxx += dx * dx; } + double slope = sxx > 0 ? sxy / sxx : 1.0; + double intercept = my - slope * mx; + double residSs = 0; + for (int i = 0; i < n; i++) { double e = obsRt[i] - (slope * predRt[i] + intercept); residSs += e * e; } + double residSd = Math.Sqrt(residSs / n); + double rtWindow = Math.Max(1.0, 2.0 * residSd); // ±2σ; floor 1 min + + Status($"S1 learned: precursor ±{tolPpm:F1} ppm (μ={pMean:F2},σ={pSd:F2}); " + + $"RT obs={slope:F3}·iRT+{intercept:F2} residSD={residSd:F2}min → window ±{rtWindow:F1}min (n={n}).", taskId); + return new MslCandidateCalibration(tolPpm, slope, intercept, rtWindow); + } + catch (Exception ex) + { + Status($"S1 calibration failed ({ex.GetType().Name}: {ex.Message}); safe fallback (no RT filter).", taskId); + return fallback; + } + } + /// /// Populates and returns the spectral match array using classic search engine /// - private void PerformSearch(List proteinsToSearch, SpectralMatch[] spectralMatchArray, List nestedIds, out HashSet updatedPsmIndexes, out int peptidesSearched, bool useCopyOnWrite = false) + private void PerformSearch(List proteinsToSearch, SpectralMatch[] spectralMatchArray, List nestedIds, out HashSet updatedPsmIndexes, out int peptidesSearched, bool useCopyOnWrite = false, List<(IBioPolymerWithSetMods Peptide, List Fragments)> precomputedPeptides = null) { + // For a .msl search, use the LEARNED precursor tolerance (S1) so the search engine and the + // candidate pre-filter agree — the calibration is trusted over the looser configured tolerance. + MzLibUtil.Tolerance precursorTolerance = CommonParameters.PrecursorMassTolerance; + if (precomputedPeptides != null && _mslCalibration.HasValue) + precursorTolerance = new MzLibUtil.PpmTolerance(_mslCalibration.Value.PrecursorTolPpm); + var massDiffAcceptor = GetMassDiffAcceptor( - CommonParameters.PrecursorMassTolerance, + precursorTolerance, SearchParameters.MassDiffAcceptorType, SearchParameters.CustomMdac); - // Run the classic search engine + // Run the classic search engine. When precomputedPeptides is supplied (from a .msl library), + // the engine iterates those peptides directly instead of digesting proteinsToSearch. var searchEngine = new TransientClassicSearchEngine( spectralMatchArray, AllSortedMs2Scans, VariableModifications, FixedModifications, proteinsToSearch, massDiffAcceptor, CommonParameters, - FileSpecificParameters, nestedIds, copyOnWriteEnabled: useCopyOnWrite); + FileSpecificParameters, nestedIds, copyOnWriteEnabled: useCopyOnWrite, precomputedPeptides: precomputedPeptides); var results = searchEngine.Run(); updatedPsmIndexes = (results as TransientSearchEngineResults)!.UpdatedSpectralMatchIndexes; @@ -750,14 +1143,23 @@ private void PerformSearch(List proteinsToSearch, SpectralMatch[] s private void InitializeCompletedDatabaseWriter() { + // The completed-database output (per-db folder + PSM/peptide/results files) was written by ONE + // consumer draining a capacity-2 channel, so at 1000s of databases the parallel searchers blocked + // on a single serial writer (low CPU utilization). Run several writer consumers — each writes a + // different database's folder, so it's safe (the shared checkpoint/progress bookkeeping is locked) + // — and widen the channel so searchers don't stall waiting for a write slot. + int writerCount = Math.Max(2, Environment.ProcessorCount / 4); _completedDatabaseWriteChannel = Channel.CreateBounded<(TransientDatabaseContext Context, TransientDatabaseMetrics Metrics)>( - new BoundedChannelOptions(2) + new BoundedChannelOptions(Math.Max(writerCount * 4, 32)) { - SingleReader = true, + SingleReader = false, // multiple parallel writer consumers SingleWriter = false, FullMode = BoundedChannelFullMode.Wait }); - _completedDatabaseWriterTask = Task.Run(RunCompletedDatabaseWriterLoopAsync); + var writers = new Task[writerCount]; + for (int i = 0; i < writerCount; i++) + writers[i] = Task.Run(RunCompletedDatabaseWriterLoopAsync); + _completedDatabaseWriterTask = Task.WhenAll(writers); } private void CompleteCompletedDatabaseWriter() @@ -796,6 +1198,26 @@ private async Task WriteCompletedDatabaseOutputsAsync(TransientDatabaseContext c List nestedIds = context.NestedIds; bool writeAllResults = !ParallelSearchParameters.WriteTransientResultsOnly; + // Skip per-database output entirely when the database produced no transient PSMs. At 1000s of + // databases the vast majority match nothing; writing a folder + header-only files for each was a + // dominant cost. The database is still recorded in the cross-database summary/checkpoint below, so + // it counts as processed (and resume still works — HasCachedResults reads the checkpoint, not disk). + bool hasTransientOutput = (context.TransientPsms?.Count ?? 0) > 0; + if (!hasTransientOutput) + { + _resultsManager!.AppendCheckpoint(metrics); + MarkDatabaseCompleted(); + UpdateProgress(TotalDatabases, nestedIds[0]); + ReportDatabaseDashboard(nestedIds[0], ParallelSearchDashboardUpdateKind.DatabaseFinished, dbName, + $"Finished {dbName} (no transient hits)", DashboardDatabaseFinishedProgress); + ReportProgress(new(100, $"Finished {dbName} (no transient hits)", nestedIds)); + return; + } + + // Create the output folder lazily — only databases with transient hits get one. + if (!Directory.Exists(outputFolder)) + Directory.CreateDirectory(outputFolder); + Status($"Writing results for {dbName}...", nestedIds); ReportDatabaseDashboard(nestedIds[0], ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName, $"Writing results for {dbName}...", DashboardDatabaseWritingProgress); diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs b/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs index 5fcc31afd..85225e0f3 100644 --- a/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs +++ b/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs @@ -121,6 +121,9 @@ private static void SetFamilyFeature(double[] vec, int index, double pValue) vec[index] = NegLog10WithSentinel(pValue, AnomalySentinels.MissingPValue); } + /// Replaces NaN/Infinity with 0 so feature vectors are always valid IsolationForest input. + private static double SanitizeFinite(double v) => double.IsNaN(v) || double.IsInfinity(v) ? 0.0 : v; + private static double NegLog10WithSentinel(double value, double sentinel) { if (double.IsNaN(value) || double.IsInfinity(value) || value <= 0) @@ -144,11 +147,15 @@ private static List ScaleAndBuild( for (int f = 0; f < featureCount; f++) { - var values = rawFeatures.Select(r => r.RawFeatures[f]).ToList(); + // Sanitize raw values first (a non-null Inf effect size, or a 0/0 ratio, can slip through the + // per-feature sentinels) so the robust-scaling statistics aren't corrupted by NaN/Infinity. + var values = rawFeatures.Select(r => SanitizeFinite(r.RawFeatures[f])).ToList(); var scaledValues = ScaleFeature(values).ToList(); + // And sanitize the scaled output — IsolationForestInput rejects any NaN/Infinity, and with + // thousands of (mostly empty) databases a single stray value would abort the whole analysis. for (int i = 0; i < rawFeatures.Count; i++) - scaled[i][f] = scaledValues[i]; + scaled[i][f] = SanitizeFinite(scaledValues[i]); } var result = new List(rawFeatures.Count); diff --git a/MetaMorpheus/TaskLayer/TaskLayer.csproj b/MetaMorpheus/TaskLayer/TaskLayer.csproj index e2b793186..04ded7be9 100644 --- a/MetaMorpheus/TaskLayer/TaskLayer.csproj +++ b/MetaMorpheus/TaskLayer/TaskLayer.csproj @@ -1,4 +1,4 @@ - + net8.0 @@ -39,7 +39,7 @@ - + diff --git a/MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs b/MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs new file mode 100644 index 000000000..0401b8e27 --- /dev/null +++ b/MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs @@ -0,0 +1,111 @@ +using Chemistry; +using NUnit.Framework; +using Priors = EngineLayer.ParallelSearch.MslPeptideReader.CandidatePriors; +using Reader = EngineLayer.ParallelSearch.MslPeptideReader; + +namespace Test.ParallelSearchTask; + +/// +/// Tests the .msl candidate pre-filter — the learned mass + RT gate that decides which library entries are +/// worth fragmenting/scoring — plus the merged-index "db|accession" parse and the lower-bound helper. +/// +[TestFixture] +public class MslCandidateFilterTests +{ + private const int Charge = 2; + private const double Mz = 500.0; + + // Build priors whose scan list has exactly one mass-matching scan (the middle one) at a chosen RT. + private static Priors PriorsWithMatch(double matchingScanRt, double tolPpm = 10, double rtWindowMin = 5, + double rtSlope = 1, double rtIntercept = 0) + { + double neutral = Mz.ToMass(Charge); + var scanMasses = new[] { neutral - 100.0, neutral, neutral + 100.0 }; // only the middle one is within tol + var scanRts = new[] { 0.0, matchingScanRt, 0.0 }; + return new Priors(scanMasses, scanRts, tolPpm, rtSlope, rtIntercept, rtWindowMin); + } + + [Test] + public void IsCandidate_MassMatch_UnpredictedRt_KeepsOnMassAlone() + { + var priors = PriorsWithMatch(matchingScanRt: 20.0); + // irt == 0 -> RT filter disabled; mass matches -> candidate. + bool cand = Reader.IsCandidate(Mz, Charge, 0f, priors, out bool massMatched); + Assert.That(cand, Is.True); + Assert.That(massMatched, Is.True); + } + + [Test] + public void IsCandidate_NoMassMatch_NotCandidate() + { + var priors = PriorsWithMatch(matchingScanRt: 20.0); + // m/z + 25 -> neutral + 50, which lands in the 100-Da gap between scan masses -> no mass match. + bool cand = Reader.IsCandidate(Mz + 25.0, Charge, 0f, priors, out bool massMatched); + Assert.That(cand, Is.False); + Assert.That(massMatched, Is.False); + } + + [Test] + public void IsCandidate_MassMatch_RtInWindow_IsCandidate() + { + var priors = PriorsWithMatch(matchingScanRt: 20.0); // slope 1, intercept 0 -> predRt == iRT + bool cand = Reader.IsCandidate(Mz, Charge, 20.0f, priors, out bool massMatched); // predRt 20 vs scan RT 20 + Assert.That(cand, Is.True); + Assert.That(massMatched, Is.True); + } + + [Test] + public void IsCandidate_MassMatch_RtOutOfWindow_NotCandidate_ButMassMatched() + { + var priors = PriorsWithMatch(matchingScanRt: 20.0, rtWindowMin: 5); + // predRt 50 vs the only mass-matching scan's RT 20 -> |50-20|=30 > 5 -> RT fails. + bool cand = Reader.IsCandidate(Mz, Charge, 50.0f, priors, out bool massMatched); + Assert.That(cand, Is.False); + Assert.That(massMatched, Is.True, "mass matched even though RT did not"); + } + + [Test] + public void IsCandidate_MassTolerance_BoundaryBehaviour() + { + double neutral = Mz.ToMass(Charge); + // The effective window is (ppm * mass) + a 0.01 Da margin. A scan 0.03 Da away is outside a 10 ppm + // window (~0.02 Da total) but inside a 40 ppm one (~0.05 Da total). + double offset = 0.03; + var priorsNarrow = new Priors(new[] { neutral + offset }, new[] { 0.0 }, precursorTolPpm: 10, rtSlope: 1, rtIntercept: 0, rtWindowMin: 5); + Assert.That(Reader.IsCandidate(Mz, Charge, 0f, priorsNarrow, out _), Is.False); + + var priorsWide = new Priors(new[] { neutral + offset }, new[] { 0.0 }, precursorTolPpm: 40, rtSlope: 1, rtIntercept: 0, rtWindowMin: 5); + Assert.That(Reader.IsCandidate(Mz, Charge, 0f, priorsWide, out _), Is.True); + } + + [Test] + public void LowerBound_ReturnsFirstIndexAtOrAboveValue() + { + var sorted = new[] { 1.0, 3.0, 5.0, 9.0 }; + Assert.Multiple(() => + { + Assert.That(Reader.LowerBound(sorted, 0.0), Is.EqualTo(0)); // below all + Assert.That(Reader.LowerBound(sorted, 3.0), Is.EqualTo(1)); // exact match + Assert.That(Reader.LowerBound(sorted, 4.0), Is.EqualTo(2)); // between -> next-higher index + Assert.That(Reader.LowerBound(sorted, 9.0), Is.EqualTo(3)); // exact last + Assert.That(Reader.LowerBound(sorted, 10.0), Is.EqualTo(4)); // above all -> length + }); + } + + [Test] + public void ParseDbTagAndAccession_SplitsOnFirstBar() + { + Assert.Multiple(() => + { + Assert.That(Reader.ParseDbTagAndAccession("UP000464024|P0DTC9"), Is.EqualTo(("UP000464024", "P0DTC9"))); + // Only the FIRST bar splits; later bars stay in the accession. + Assert.That(Reader.ParseDbTagAndAccession("db|acc|extra"), Is.EqualTo(("db", "acc|extra"))); + // No bar -> UNKNOWN db, whole string is the accession. + Assert.That(Reader.ParseDbTagAndAccession("loneAccession"), Is.EqualTo(("UNKNOWN", "loneAccession"))); + // Empty accession after the bar -> "_UNKNOWN". + Assert.That(Reader.ParseDbTagAndAccession("db|"), Is.EqualTo(("db", "db_UNKNOWN"))); + // Null/empty input. + Assert.That(Reader.ParseDbTagAndAccession(null), Is.EqualTo(("UNKNOWN", "UNKNOWN_UNKNOWN"))); + }); + } +} diff --git a/MetaMorpheus/Test/Test.csproj b/MetaMorpheus/Test/Test.csproj index 37d8f8a6d..4ae5cd0f0 100644 --- a/MetaMorpheus/Test/Test.csproj +++ b/MetaMorpheus/Test/Test.csproj @@ -1,4 +1,4 @@ - + net8.0-windows @@ -24,7 +24,7 @@ - + From fda60e74ff9e0e9ca6c24135635734b622ab8b84 Mon Sep 17 00:00:00 2001 From: trishorts Date: Mon, 15 Jun 2026 12:31:30 -0500 Subject: [PATCH 2/4] chore(csproj): restore UTF-8 BOM dropped alongside the mzLib bump The .msl dependency bump stripped the byte-order mark from six project files as unrelated churn; restore it so each csproj diff is just the mzLib version line. Addresses manual review (C5) on nbollis/MetaMorpheus#9. Co-Authored-By: Claude Opus 4.8 (1M context) --- MetaMorpheus/CMD/CMD.csproj | 2 +- MetaMorpheus/EngineLayer/EngineLayer.csproj | 2 +- MetaMorpheus/GUI/GUI.csproj | 2 +- MetaMorpheus/GuiFunctions/GuiFunctions.csproj | 2 +- MetaMorpheus/TaskLayer/TaskLayer.csproj | 2 +- MetaMorpheus/Test/Test.csproj | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/MetaMorpheus/CMD/CMD.csproj b/MetaMorpheus/CMD/CMD.csproj index 94cfdf487..f459d84c5 100644 --- a/MetaMorpheus/CMD/CMD.csproj +++ b/MetaMorpheus/CMD/CMD.csproj @@ -1,4 +1,4 @@ - + Exe diff --git a/MetaMorpheus/EngineLayer/EngineLayer.csproj b/MetaMorpheus/EngineLayer/EngineLayer.csproj index 1baac2c68..987ce2d8f 100644 --- a/MetaMorpheus/EngineLayer/EngineLayer.csproj +++ b/MetaMorpheus/EngineLayer/EngineLayer.csproj @@ -1,4 +1,4 @@ - + net8.0 diff --git a/MetaMorpheus/GUI/GUI.csproj b/MetaMorpheus/GUI/GUI.csproj index 658d8ec0e..c9da4341c 100644 --- a/MetaMorpheus/GUI/GUI.csproj +++ b/MetaMorpheus/GUI/GUI.csproj @@ -1,4 +1,4 @@ - + WinExe diff --git a/MetaMorpheus/GuiFunctions/GuiFunctions.csproj b/MetaMorpheus/GuiFunctions/GuiFunctions.csproj index 7dd533525..ad079f318 100644 --- a/MetaMorpheus/GuiFunctions/GuiFunctions.csproj +++ b/MetaMorpheus/GuiFunctions/GuiFunctions.csproj @@ -1,4 +1,4 @@ - + net8.0-windows diff --git a/MetaMorpheus/TaskLayer/TaskLayer.csproj b/MetaMorpheus/TaskLayer/TaskLayer.csproj index 04ded7be9..a1fa9df37 100644 --- a/MetaMorpheus/TaskLayer/TaskLayer.csproj +++ b/MetaMorpheus/TaskLayer/TaskLayer.csproj @@ -1,4 +1,4 @@ - + net8.0 diff --git a/MetaMorpheus/Test/Test.csproj b/MetaMorpheus/Test/Test.csproj index 4ae5cd0f0..db8a42080 100644 --- a/MetaMorpheus/Test/Test.csproj +++ b/MetaMorpheus/Test/Test.csproj @@ -1,4 +1,4 @@ - + net8.0-windows From 2cc6b1d780ff480806262907240ab083c2cea6c1 Mon Sep 17 00:00:00 2001 From: trishorts Date: Mon, 15 Jun 2026 12:31:40 -0500 Subject: [PATCH 3/4] fix(parallelsearch): stop abandoning the loader when a transient search faults If a consumer (search) threw, Parallel.ForEach propagated out of the loop and the producer's GetResult() was skipped: the load task could stay blocked on a full bounded queue until disposal, and its exception went unobserved. Add a CancellationTokenSource the consumer cancels on fault so the producer's Add unwinds promptly, and always observe the producer task in finally. Also replace the MM_PARALLELSEARCH_MERGED environment variable with a real ParallelSearchParameters.UseMergedTransientLibrary flag (parameter/TOML/CLI surface, reproducible) and document where the multi-writer bookkeeping is synchronized (AppendToFile / MarkDatabaseCompleted / UpdateProgress locks). Addresses manual review (C1, A4, C2) on nbollis/MetaMorpheus#9. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../ParallelSearchParameters.cs | 10 +++++++ .../ParallelSearch/ParallelSearchTask.cs | 30 ++++++++++++++----- .../ParallelSearchParametersTests.cs | 1 + 3 files changed, 34 insertions(+), 7 deletions(-) diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs index 7bd49d408..1412683fd 100644 --- a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs +++ b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs @@ -16,6 +16,16 @@ public class ParallelSearchParameters : SearchParameters public bool CompressTransientSearchOutputs { get; set; } = false; public string? DeNovoMappingDataFilePath { get; set; } = null; + /// + /// When true, a single .msl entry in is treated as a merged index: one + /// file holds many databases (each entry tagged "db|accession"). Merged mode loads the file once, runs the + /// candidate filter once, and emits one searchable database per source db-group, so parallelism is driven + /// by the database count INSIDE the file rather than the file-list count. Replaces the former + /// MM_PARALLELSEARCH_MERGED environment variable so the behavior is on the parameter/TOML/CLI + /// surface and reproducible. Only takes effect when the single transient database ends in .msl. + /// + public bool UseMergedTransientLibrary { get; set; } = false; + #region Follow-Up Search Parameters /// diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs index 18fa371b3..33f81d0d6 100644 --- a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs +++ b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs @@ -188,7 +188,7 @@ protected override MyTaskResults RunSpecific(string outputFolder, // MERGED-INDEX mode: one .msl file holds many databases (entries tagged "db|accession"). The number // of databases to search comes from inside the file, not the (count==1) file list, so parallelism // must not be capped by the file count. - bool mergedMode = Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_MERGED") == "1" + bool mergedMode = ParallelSearchParameters.UseMergedTransientLibrary && ParallelSearchParameters.TransientDatabases.Count == 1 && ParallelSearchParameters.TransientDatabases[0].FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase); @@ -207,6 +207,9 @@ protected override MyTaskResults RunSpecific(string outputFolder, // The result WRITE channel was already pipelined; this pipelines the read side too. The bounded // capacity caps how many loaded databases are held in memory at once. var swTransientLoop = Stopwatch.StartNew(); + // Abort signal: if a CONSUMER (search) faults, cancel so the producer stops blocking on Add into a + // bounded queue that is no longer being drained, instead of being abandoned until queue disposal. + using (var loadAbort = new System.Threading.CancellationTokenSource()) using (var loadedQueue = new System.Collections.Concurrent.BlockingCollection( boundedCapacity: Math.Max(2, databaseParallelism))) { @@ -233,7 +236,7 @@ protected override MyTaskResults RunSpecific(string outputFolder, { if (GlobalVariables.StopLoops) return; var loaded = BuildLoadedFromCandidates(kvp.Key, kvp.Value, outputFolder, taskId); - if (loaded != null) loadedQueue.Add(loaded); + if (loaded != null) loadedQueue.Add(loaded, loadAbort.Token); }); } else @@ -254,7 +257,7 @@ protected override MyTaskResults RunSpecific(string outputFolder, throw; } if (loaded != null) - loadedQueue.Add(loaded); + loadedQueue.Add(loaded, loadAbort.Token); }); } } @@ -283,12 +286,22 @@ protected override MyTaskResults RunSpecific(string outputFolder, } }); } + catch + { + // A searcher faulted and Parallel.ForEach stopped draining the queue. Signal the producer so a + // load thread blocked on the bounded Add unwinds promptly instead of being abandoned; then rethrow. + loadAbort.Cancel(); + throw; + } finally { CompleteCompletedDatabaseWriter(); + // ALWAYS observe the producer task: on the happy path this surfaces a load exception; on the + // failure path it prevents an unobserved task / abandoned loader. When WE triggered the abort + // above, the producer's resulting cancellation is a consequence, not the root cause — swallow it. + try { producer.GetAwaiter().GetResult(); } + catch (Exception) when (loadAbort.IsCancellationRequested) { } } - - producer.GetAwaiter().GetResult(); // surface any producer (load) exception } swTransientLoop.Stop(); LogPhaseTimingBreakdown(taskId, outputFolder, swInit.Elapsed, swTransientLoop.Elapsed, databaseParallelism); @@ -1146,8 +1159,11 @@ private void InitializeCompletedDatabaseWriter() // The completed-database output (per-db folder + PSM/peptide/results files) was written by ONE // consumer draining a capacity-2 channel, so at 1000s of databases the parallel searchers blocked // on a single serial writer (low CPU utilization). Run several writer consumers — each writes a - // different database's folder, so it's safe (the shared checkpoint/progress bookkeeping is locked) - // — and widen the channel so searchers don't stall waiting for a write slot. + // different database's folder, so file writes never collide — and widen the channel so searchers + // don't stall waiting for a write slot. The three pieces of shared bookkeeping the consumers touch + // are each synchronized: AppendCheckpoint -> ParallelSearchResultCache.AppendToFile (lock _writeLock + // + _cacheLock over a ConcurrentDictionary), MarkDatabaseCompleted (lock _dashboardLock), and + // UpdateProgress (lock _progressLock). int writerCount = Math.Max(2, Environment.ProcessorCount / 4); _completedDatabaseWriteChannel = Channel.CreateBounded<(TransientDatabaseContext Context, TransientDatabaseMetrics Metrics)>( new BoundedChannelOptions(Math.Max(writerCount * 4, 32)) diff --git a/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs b/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs index 7205f30f9..274f55c2e 100644 --- a/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs +++ b/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs @@ -20,6 +20,7 @@ public void Constructor_SetsParallelSearchDefaults() Assert.That(parameters.MaxSearchesInParallel, Is.EqualTo(4)); Assert.That(parameters.WriteTransientResultsOnly, Is.True); Assert.That(parameters.WriteTransientSpectralLibrary, Is.False); + Assert.That(parameters.UseMergedTransientLibrary, Is.False); Assert.That(parameters.DoParsimony, Is.True); Assert.That(parameters.NoOneHitWonders, Is.True); Assert.That(parameters.MassDiffAcceptorType, Is.EqualTo(MassDiffAcceptorType.Exact)); From 9f6f5d606fa61e71c14bfc8739155b68ecb2ef67 Mon Sep 17 00:00:00 2001 From: trishorts Date: Mon, 15 Jun 2026 13:03:10 -0500 Subject: [PATCH 4/4] refactor(parallelsearch): remove the backend-swappable scorer abstraction The ISpectralScorer / SpectralScorerProvider / ScoringBatch / FragmentMatch / IScoringSink / SpectralScoringData layer existed to host a GPU fragment-matching backend. GPU benchmarks slower and there is no concrete second backend, so the abstraction is pure surface area. Collapse the transient search inner loop back onto the shared MatchFragmentIons / CalculatePeptideScore so MatchedFragmentIon is again the type at the engine boundary; the .msl precomputed-fragments path and the shared-base-PSM / copy-on-write optimizations are unchanged. Behavior matches the original pre-abstraction engine (same closest-envelope match, 1 + intensity/TIC scoring, count-based ScoreCutoff). Addresses manual review (A1-A3) on nbollis/MetaMorpheus#9. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../Scoring/CpuSpectralScorer.cs | 80 ------ .../ParallelSearch/Scoring/ISpectralScorer.cs | 95 ------- .../Scoring/SpectralScorerProvider.cs | 26 -- .../Scoring/SpectralScoringData.cs | 147 ----------- .../TransientClassicSearchEngine.cs | 236 ++---------------- 5 files changed, 23 insertions(+), 561 deletions(-) delete mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs delete mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs delete mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs delete mode 100644 MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs deleted file mode 100644 index 3c2e6c078..000000000 --- a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/CpuSpectralScorer.cs +++ /dev/null @@ -1,80 +0,0 @@ -using System; -using Chemistry; -using MzLibUtil; - -namespace EngineLayer.ParallelSearch.Scoring -{ - /// - /// CPU implementation of and the correctness baseline for the - /// GPU path. Reproduces the original TransientClassicSearchEngine matching exactly: for each - /// theoretical fragment, find the closest experimental envelope (binary search), accept it if - /// it is within the product tolerance and its charge does not exceed the precursor charge. - /// - public sealed class CpuSpectralScorer : ISpectralScorer - { - private readonly SpectralScoringData _data; - private readonly Tolerance _productTolerance; - private FragmentMatch[] _matchBuffer = new FragmentMatch[64]; - - public CpuSpectralScorer(SpectralScoringData data, Tolerance productTolerance) - { - _data = data ?? throw new ArgumentNullException(nameof(data)); - _productTolerance = productTolerance ?? throw new ArgumentNullException(nameof(productTolerance)); - } - - public string BackendDescription => "CPU (binary search per fragment)"; - - public void ScoreBatch(ScoringBatch batch, IScoringSink sink) - { - if (batch == null) throw new ArgumentNullException(nameof(batch)); - if (sink == null) throw new ArgumentNullException(nameof(sink)); - - for (int w = 0; w < batch.WorkItemCount; w++) - { - int slot = batch.WorkPeptideSlot[w]; - int scanIndex = batch.WorkScanIndex[w]; - - int fragStart = batch.PeptideFragmentOffsets[slot]; - int fragEnd = batch.PeptideFragmentOffsets[slot + 1]; - int fragCount = fragEnd - fragStart; - - if (_matchBuffer.Length < fragCount) - _matchBuffer = new FragmentMatch[Math.Max(fragCount, _matchBuffer.Length * 2)]; - - int precursorCharge = _data.ScanPrecursorCharges[scanIndex]; - int matchCount = 0; - - for (int local = 0; local < fragCount; local++) - { - double theoreticalMass = batch.FragmentNeutralMasses[fragStart + local]; - - // unknown fragment mass; rare, for sequences with unknown amino acids - if (double.IsNaN(theoreticalMass)) - continue; - - int idx = _data.GetClosestFragmentIndex(scanIndex, theoreticalMass); - if (idx < 0) - continue; - - double experimentalMass = _data.FragmentMonoMasses[idx]; - int experimentalCharge = _data.FragmentCharges[idx]; - - if (_productTolerance.Within(experimentalMass, theoreticalMass) - && Math.Abs(experimentalCharge) <= Math.Abs(precursorCharge)) - { - _matchBuffer[matchCount++] = new FragmentMatch( - local, - experimentalMass.ToMz(experimentalCharge), - _data.FragmentIntensities[idx], - experimentalCharge); - } - } - - if (matchCount > 0) - sink.AcceptWorkItem(w, _matchBuffer, matchCount); - } - } - - public void Dispose() { /* nothing to release on the CPU path */ } - } -} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs deleted file mode 100644 index 521a1e5c4..000000000 --- a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/ISpectralScorer.cs +++ /dev/null @@ -1,95 +0,0 @@ -using System; - -namespace EngineLayer.ParallelSearch.Scoring -{ - /// - /// Abstraction over the hot inner step of the transient classic search: for each - /// (peptide, candidate-scan) work item, decide which theoretical fragments match an - /// experimental fragment within tolerance and the charge filter. - /// - /// Deliberately narrow: the scorer does ONLY the matching (the expensive binary-search + - /// tolerance work that is identical across all 30k databases). Score computation, the - /// ScoreCutoff gate, MatchedFragmentIon construction and PSM updates all stay in the - /// engine. - /// - /// Implementation: - /// CpuSpectralScorer — binary search per fragment on the CPU. - /// - /// One instance per thread (not safe for concurrent ScoreBatch calls on the same instance). - /// - public interface ISpectralScorer : IDisposable - { - /// - /// Match every work item in , reporting each item's matched - /// fragments to . Items with no matches need not be reported. - /// - void ScoreBatch(ScoringBatch batch, IScoringSink sink); - - /// Human-readable backend description (for startup logging). - string BackendDescription { get; } - } - - /// One matched theoretical-vs-experimental fragment pairing. - public readonly struct FragmentMatch - { - /// Index of the theoretical product within its peptide's product list. - public readonly int LocalProductIndex; - /// Experimental m/z = monoisotopic mass converted at the envelope charge. - public readonly double ExperimentalMz; - /// Representative (first-peak) intensity of the matched experimental envelope. - public readonly double ExperimentalIntensity; - /// Charge of the matched experimental envelope. - public readonly int ExperimentalCharge; - - public FragmentMatch(int localProductIndex, double experimentalMz, double experimentalIntensity, int experimentalCharge) - { - LocalProductIndex = localProductIndex; - ExperimentalMz = experimentalMz; - ExperimentalIntensity = experimentalIntensity; - ExperimentalCharge = experimentalCharge; - } - } - - /// - /// Receives match results one work item at a time. The array is a - /// reusable buffer owned by the scorer — the sink must consume entries [0, matchCount) during - /// the call and not retain the array. - /// - public interface IScoringSink - { - void AcceptWorkItem(int workIndex, FragmentMatch[] matches, int matchCount); - } - - /// - /// A batch of (peptide, candidate-scan) work items to match, in struct-of-arrays form. - /// Theoretical fragment masses are concatenated per peptide slot (CSR via - /// ); NaN masses are kept in place so - /// aligns with the peptide's product list (the - /// scorer skips NaN, which never matches). - /// - public sealed class ScoringBatch - { - /// Number of work items. - public int WorkItemCount; - - /// Peptide slot for each work item (index into PeptideFragmentOffsets). - public int[] WorkPeptideSlot = Array.Empty(); - /// Candidate scan index for each work item (index into the SpectralScoringData). - public int[] WorkScanIndex = Array.Empty(); - - /// Number of distinct peptide slots in this batch. - public int PeptideCount; - /// CSR offsets into ; length = PeptideCount + 1. - public int[] PeptideFragmentOffsets = new int[1]; - /// Concatenated theoretical product neutral masses for all peptide slots. - public double[] FragmentNeutralMasses = Array.Empty(); - - /// Reset counters for reuse without reallocating the backing arrays. - public void Clear() - { - WorkItemCount = 0; - PeptideCount = 0; - PeptideFragmentOffsets[0] = 0; - } - } -} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs deleted file mode 100644 index f80f56112..000000000 --- a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScorerProvider.cs +++ /dev/null @@ -1,26 +0,0 @@ -using MzLibUtil; - -namespace EngineLayer.ParallelSearch.Scoring -{ - /// - /// Hands out an to each search partition (thread). - /// Each thread gets its own lightweight scorer (per-thread match buffer, no shared state). - /// - public interface ISpectralScorerProvider : System.IDisposable - { - /// Scorer for the calling partition. Do NOT dispose the returned scorer; the provider owns it. - ISpectralScorer GetScorer(); - string BackendDescription { get; } - } - - /// CPU provider: a fresh (cheap) CpuSpectralScorer per partition. - public sealed class CpuScorerProvider : ISpectralScorerProvider - { - private readonly SpectralScoringData _data; - private readonly Tolerance _tolerance; - public CpuScorerProvider(SpectralScoringData data, Tolerance tolerance) { _data = data; _tolerance = tolerance; } - public ISpectralScorer GetScorer() => new CpuSpectralScorer(_data, _tolerance); - public string BackendDescription => "CPU (binary search per fragment)"; - public void Dispose() { } - } -} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs b/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs deleted file mode 100644 index 77061dc5c..000000000 --- a/MetaMorpheus/EngineLayer/ParallelSearch/Scoring/SpectralScoringData.cs +++ /dev/null @@ -1,147 +0,0 @@ -using System; -using System.Linq; -using Chemistry; -using MassSpectrometry; - -namespace EngineLayer.ParallelSearch.Scoring -{ - /// - /// Read-only, flattened (struct-of-arrays) view of the experimental MS2 spectra used by - /// the transient parallel search. Built ONCE from the shared, precursor-mass-sorted - /// array and reused, read-only, across every - /// transient database search. - /// - /// This is the data the GPU keeps resident in VRAM: experimental fragment masses are - /// concatenated into one sorted-per-scan array (CSR layout via ), - /// with parallel charge/intensity arrays and per-scan TIC / precursor charge. The CPU scorer - /// reconstructs s from the same data; the GPU scorer reads only - /// the flat arrays. - /// - /// Per-scan fragment slices are sorted ascending by monoisotopic mass (the - /// constructor guarantees ExperimentalFragments is sorted), - /// which is what makes the per-fragment closest-mass lookup a binary search. - /// - public sealed class SpectralScoringData - { - /// The original scans, kept for CPU-side MatchedFragmentIon reconstruction. - public Ms2ScanWithSpecificMass[] Scans { get; } - - /// Precursor masses, ascending (mirrors ClassicSearchEngine.MyScanPrecursorMasses). - public double[] ScanPrecursorMasses { get; } - - /// CSR offsets into the flat fragment arrays; length = Scans.Length + 1. - public int[] ScanFragmentOffsets { get; } - - /// Experimental fragment monoisotopic masses, concatenated; each scan's slice is sorted ascending. - public double[] FragmentMonoMasses { get; } - - /// Charge of each experimental fragment envelope (parallel to FragmentMonoMasses). - public int[] FragmentCharges { get; } - - /// Representative intensity (first peak) of each experimental fragment envelope. - public double[] FragmentIntensities { get; } - - /// Total ion current per scan (scoring normalization). - public double[] ScanTotalIonCurrents { get; } - - /// Precursor charge per scan (fragment-charge acceptance filter). - public int[] ScanPrecursorCharges { get; } - - public int ScanCount => Scans.Length; - - private SpectralScoringData( - Ms2ScanWithSpecificMass[] scans, double[] scanPrecursorMasses, int[] scanFragmentOffsets, - double[] fragmentMonoMasses, int[] fragmentCharges, double[] fragmentIntensities, - double[] scanTotalIonCurrents, int[] scanPrecursorCharges) - { - Scans = scans; - ScanPrecursorMasses = scanPrecursorMasses; - ScanFragmentOffsets = scanFragmentOffsets; - FragmentMonoMasses = fragmentMonoMasses; - FragmentCharges = fragmentCharges; - FragmentIntensities = fragmentIntensities; - ScanTotalIonCurrents = scanTotalIonCurrents; - ScanPrecursorCharges = scanPrecursorCharges; - } - - /// - /// Flatten the (already precursor-mass-sorted) scans into struct-of-arrays form. - /// One-time cost amortized over all transient database searches. - /// - public static SpectralScoringData Build(Ms2ScanWithSpecificMass[] sortedScans) - { - if (sortedScans == null) throw new ArgumentNullException(nameof(sortedScans)); - - int scanCount = sortedScans.Length; - var precursorMasses = new double[scanCount]; - var offsets = new int[scanCount + 1]; - var tics = new double[scanCount]; - var precursorCharges = new int[scanCount]; - - int totalFragments = 0; - for (int i = 0; i < scanCount; i++) - { - var frags = sortedScans[i].ExperimentalFragments; - totalFragments += frags?.Length ?? 0; - } - - var monoMasses = new double[totalFragments]; - var charges = new int[totalFragments]; - var intensities = new double[totalFragments]; - - int write = 0; - for (int i = 0; i < scanCount; i++) - { - var scan = sortedScans[i]; - offsets[i] = write; - precursorMasses[i] = scan.PrecursorMass; - tics[i] = scan.TotalIonCurrent; - precursorCharges[i] = scan.PrecursorCharge; - - var frags = scan.ExperimentalFragments; - if (frags != null) - { - for (int f = 0; f < frags.Length; f++) - { - var env = frags[f]; - monoMasses[write] = env.MonoisotopicMass; - charges[write] = env.Charge; - // Mirrors the existing engine: closestExperimentalMass.Peaks.First().intensity - intensities[write] = env.Peaks.First().intensity; - write++; - } - } - } - offsets[scanCount] = write; - - return new SpectralScoringData(sortedScans, precursorMasses, offsets, - monoMasses, charges, intensities, tics, precursorCharges); - } - - /// - /// Index of the experimental fragment whose mass is closest to - /// within scan 's slice, or -1 if the scan has no fragments. - /// Reproduces Ms2ScanWithSpecificMass.GetClosestFragmentMass over the flattened slice. - /// - public int GetClosestFragmentIndex(int scanIndex, double mass) - { - int lo = ScanFragmentOffsets[scanIndex]; - int hi = ScanFragmentOffsets[scanIndex + 1]; - int len = hi - lo; - if (len == 0) - return -1; - - int index = Array.BinarySearch(FragmentMonoMasses, lo, len, mass); - if (index >= 0) - return index; - - index = ~index; - - if (index == hi) - return index - 1; - if (index == lo || mass - FragmentMonoMasses[index - 1] > FragmentMonoMasses[index] - mass) - return index; - return index - 1; - } - } -} diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs b/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs index c0a91b0be..e57a78116 100644 --- a/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs +++ b/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs @@ -2,15 +2,11 @@ using System.Collections.Concurrent; using System.Collections.Generic; using System.Linq; -using System.Runtime.CompilerServices; using System.Threading; using System.Threading.Tasks; -using Chemistry; using EngineLayer.ClassicSearch; using EngineLayer.FdrAnalysis; -using EngineLayer.ParallelSearch.Scoring; using EngineLayer.Util; -using MzLibUtil; using Omics; using Omics.Fragmentation; using Omics.Modifications; @@ -23,7 +19,7 @@ namespace EngineLayer.ParallelSearch /// The base PSM list is shared across threads, and only the PSMs that are updated by the transient search are cloned and modified, while the rest of the PSMs remain shared and read-only. /// Features that are removed from classic search engine for runtime efficiency and memory efficiency: /// - Some detailed logging - /// - Mass Differance acceptor is always exact. + /// - Mass Differance acceptor is always exact. /// public class TransientClassicSearchEngine : ClassicSearchEngine { @@ -57,13 +53,6 @@ public TransientClassicSearchEngine(SpectralMatch[] globalPsms, Ms2ScanWithSpeci } } - // The flattened experimental spectra depend only on the (shared) scan array, so cache by - // its reference: every transient engine over the same file set reuses one build. - private static readonly ConditionalWeakTable _scoringDataCache = new(); - - private static SpectralScoringData GetOrBuildScoringData(Ms2ScanWithSpecificMass[] sortedScans) - => _scoringDataCache.GetValue(sortedScans, SpectralScoringData.Build); - protected override MetaMorpheusEngineResults RunSpecific() { double proteinsSearched = 0; @@ -75,23 +64,13 @@ protected override MetaMorpheusEngineResults RunSpecific() bool usePrecomputedPeptides = _precomputedPeptides != null && _precomputedPeptides.Count > 0; if (Proteins.Any() || usePrecomputedPeptides) { - - // Experimental spectra are identical and read-only across every transient database - // search, so flatten them once (struct-of-arrays) and reuse across all peptides/threads. - var scoringData = GetOrBuildScoringData(ArrayOfSortedMS2Scans); - using var scorerProvider = new CpuScorerProvider( - scoringData, CommonParameters.ProductMassTolerance); - Status("ParallelSearch scoring backend: " + scorerProvider.BackendDescription); - - // Shared per-peptide work: fragment (in double), queue candidate scans, flush. The - // scorer does ONLY the fragment matching (the GPU-able hot step); the accumulator - // builds MatchedFragmentIons, applies the score cutoff, scores, and updates PSMs. - // A peptide's work items are queued contiguously so per-scan candidate ordering is - // preserved (bit-identical results). - // precomputedFragments != null (from a .msl library) skips the Fragment() call and matches - // the library's stored fragments directly — the fragmentation-time savings. - void ProcessOnePeptide(IBioPolymerWithSetMods peptide, TransientScoringBatch batch, - ISpectralScorer scorer, List peptideTheorProducts, List precomputedFragments = null) + // Match one peptide against every candidate scan and update PSMs directly, using the shared + // MatchFragmentIons / CalculatePeptideScore so MatchedFragmentIon stays the type at the engine + // boundary. The two memory/runtime wins over ClassicSearchEngine are kept: the base PSM list is + // shared (only updated PSMs are cloned — see AddPeptideCandidateToPsm), and a .msl library + // supplies precomputed fragments so the Fragment() call is skipped entirely. + void ProcessOnePeptide(IBioPolymerWithSetMods peptide, List peptideTheorProducts, + List precomputedFragments = null) { Interlocked.Increment(ref peptideCounter); @@ -107,33 +86,27 @@ void ProcessOnePeptide(IBioPolymerWithSetMods peptide, TransientScoringBatch bat products = peptideTheorProducts; } - int slot = -1; foreach (ScanWithIndexAndNotchInfo scan in GetAcceptableScans(peptide.MonoisotopicMass, SearchMode)) { - if (slot < 0) - slot = batch.BeginPeptide(peptide, products); - batch.AddWorkItem(slot, scan.ScanIndex, scan.Notch); + Ms2ScanWithSpecificMass theScan = ArrayOfSortedMS2Scans[scan.ScanIndex]; + List matchedIons = MatchFragmentIons(theScan, products, CommonParameters); + double thisScore = CalculatePeptideScore(theScan.TheScan, matchedIons); + AddPeptideCandidateToPsm(scan, thisScore, peptide, matchedIons); } - - if (batch.ShouldFlush) - batch.Flush(scorer); } - // GPU = one shared thread-safe scorer (resident spectra); CPU = a cheap per-partition - // scorer. The provider owns the scorer's lifetime. + // Protein (FASTA) peptide source: digest each protein and search each peptide. Action processProteinRange = (start, end) => { - var scorer = scorerProvider.GetScorer(); - var batch = new TransientScoringBatch(this, scoringData); List peptideTheorProducts = new(); for (int i = start; i < end; i++) { - if (GlobalVariables.StopLoops) { batch.Flush(scorer); return; } + if (GlobalVariables.StopLoops) return; // digest each protein into peptides and search each peptide in all spectra within precursor mass tolerance foreach (var specificBioPolymer in Proteins[i].Digest(CommonParameters.DigestionParams, FixedModifications, VariableModifications)) - ProcessOnePeptide(specificBioPolymer, batch, scorer, peptideTheorProducts); + ProcessOnePeptide(specificBioPolymer, peptideTheorProducts); // report search progress (proteins searched so far out of total proteins in database) proteinsSearched++; @@ -144,25 +117,19 @@ void ProcessOnePeptide(IBioPolymerWithSetMods peptide, TransientScoringBatch bat ReportProgress(new ProgressEventArgs(percentProgress, "Performing classic search... ", NestedIds)); } } - - batch.Flush(scorer); }; // Library (.msl) peptide source: iterate precomputed peptides directly, no digestion. Action processPeptideRange = (start, end) => { - var scorer = scorerProvider.GetScorer(); - var batch = new TransientScoringBatch(this, scoringData); List peptideTheorProducts = new(); for (int i = start; i < end; i++) { - if (GlobalVariables.StopLoops) { batch.Flush(scorer); return; } + if (GlobalVariables.StopLoops) return; var (peptide, fragments) = _precomputedPeptides[i]; - ProcessOnePeptide(peptide, batch, scorer, peptideTheorProducts, fragments); + ProcessOnePeptide(peptide, peptideTheorProducts, fragments); } - - batch.Flush(scorer); }; int itemCount = usePrecomputedPeptides ? _precomputedPeptides.Count : Proteins.Count; @@ -196,6 +163,10 @@ void ProcessOnePeptide(IBioPolymerWithSetMods peptide, TransientScoringBatch bat private void AddPeptideCandidateToPsm(ScanWithIndexAndNotchInfo scan, double thisScore, IBioPolymerWithSetMods peptide, List matchedIons) { + // matched ion count below the cutoff is not a candidate (mirrors ClassicSearchEngine's ScoreCutoff gate) + if (matchedIons.Count < CommonParameters.ScoreCutoff) + return; + int scanIndex = scan.ScanIndex; if (_singleThreadMode) @@ -294,8 +265,8 @@ private SpectralMatch EnsureWritablePsm(int scanIndex, SpectralMatch existingPsm if (!_copyOnWriteEnabled || existingPsm == null || UpdatedIndexes.ContainsKey(scanIndex)) { return existingPsm; - } - + } + SpectralMatches[scanIndex] = ClonePsmForWrite(existingPsm); return SpectralMatches[scanIndex]; } @@ -348,167 +319,6 @@ private int GetFirstScanWithMassOverOrEqual(double minimum) // index of the first element that is larger than value return index; } - - /// - /// Per-thread accumulator that batches (peptide, candidate-scan) work for the - /// and, on flush, turns each item's matched fragments into a - /// PSM update — reproducing the original per-scan logic (HashSet dedup, ScoreCutoff gate, - /// (1 + intensity/TIC) scoring) exactly, just deferred to flush time. Work items are queued - /// in digestion order with each peptide's items contiguous, so PSM update order (and thus - /// tie resolution) matches the unbatched engine. - /// - private sealed class TransientScoringBatch : IScoringSink - { - private const int WorkItemFlushThreshold = 16384; - private const int PeptideFlushThreshold = 4096; - - private readonly TransientClassicSearchEngine _engine; - private readonly SpectralScoringData _data; - private readonly ScoringBatch _batch = new(); - // Reused across work items / peptides to keep the batched path allocation-free in the - // hot loop (the unbatched engine reused one cleared HashSet and never copied products). - private readonly HashSet _matchedIonScratch = new(); - private IBioPolymerWithSetMods[] _slotPeptides = new IBioPolymerWithSetMods[1024]; - private int[] _slotProductOffset = new int[1025]; - private Product[] _productPool = new Product[4096]; - private int[] _workNotch = new int[1024]; - private int _fragmentWrite; - private int _productWrite; - - public TransientScoringBatch(TransientClassicSearchEngine engine, SpectralScoringData data) - { - _engine = engine; - _data = data; - } - - public bool ShouldFlush => - _batch.WorkItemCount >= WorkItemFlushThreshold || _batch.PeptideCount >= PeptideFlushThreshold; - - /// Register a peptide and its theoretical products; returns its batch slot. - public int BeginPeptide(IBioPolymerWithSetMods peptide, List products) - { - int slot = _batch.PeptideCount; - int count = products.Count; - - EnsureSlotCapacity(slot + 1); - _slotPeptides[slot] = peptide; - _slotProductOffset[slot] = _productWrite; - - // Snapshot the products (the caller reuses/clears its list) into a pooled array, - // and mirror their neutral masses into the flat batch buffer for the scorer. - EnsureProductPoolCapacity(_productWrite + count); - EnsureFragmentCapacity(_fragmentWrite + count); - for (int k = 0; k < count; k++) - { - Product p = products[k]; - _productPool[_productWrite + k] = p; - _batch.FragmentNeutralMasses[_fragmentWrite + k] = p.NeutralMass; - } - _productWrite += count; - _fragmentWrite += count; - - EnsureOffsetCapacity(slot + 2); - _batch.PeptideFragmentOffsets[slot + 1] = _fragmentWrite; - _slotProductOffset[slot + 1] = _productWrite; - _batch.PeptideCount = slot + 1; - return slot; - } - - public void AddWorkItem(int slot, int scanIndex, int notch) - { - int w = _batch.WorkItemCount; - EnsureWorkCapacity(w + 1); - _batch.WorkPeptideSlot[w] = slot; - _batch.WorkScanIndex[w] = scanIndex; - _workNotch[w] = notch; - _batch.WorkItemCount = w + 1; - } - - public void Flush(ISpectralScorer scorer) - { - if (_batch.WorkItemCount > 0) - scorer.ScoreBatch(_batch, this); - Reset(); - } - - private void Reset() - { - _batch.Clear(); - _fragmentWrite = 0; - _productWrite = 0; - } - - void IScoringSink.AcceptWorkItem(int workIndex, FragmentMatch[] matches, int matchCount) - { - int slot = _batch.WorkPeptideSlot[workIndex]; - int scanIndex = _batch.WorkScanIndex[workIndex]; - int notch = _workNotch[workIndex]; - int productBase = _slotProductOffset[slot]; - - HashSet matchedFragmentIons = _matchedIonScratch; - matchedFragmentIons.Clear(); - for (int k = 0; k < matchCount; k++) - { - FragmentMatch m = matches[k]; - matchedFragmentIons.Add(new MatchedFragmentIon( - _productPool[productBase + m.LocalProductIndex], m.ExperimentalMz, m.ExperimentalIntensity, m.ExperimentalCharge)); - } - - if (matchedFragmentIons.Count < _engine.CommonParameters.ScoreCutoff) - return; - - double tic = _data.ScanTotalIonCurrents[scanIndex]; - double score = 0; - foreach (var ion in matchedFragmentIons) - { - if (ion.NeutralTheoreticalProduct.ProductType != ProductType.D) - score += 1 + ion.Intensity / tic; - } - - var matchedIons = matchedFragmentIons.ToList(); - _engine.AddPeptideCandidateToPsm( - new ScanWithIndexAndNotchInfo(notch, scanIndex), score, _slotPeptides[slot], matchedIons); - } - - private void EnsureSlotCapacity(int neededSlots) - { - if (_slotPeptides.Length < neededSlots) - { - int n = Math.Max(neededSlots, _slotPeptides.Length * 2); - Array.Resize(ref _slotPeptides, n); - Array.Resize(ref _slotProductOffset, n + 1); - } - } - - private void EnsureProductPoolCapacity(int needed) - { - if (_productPool.Length < needed) - Array.Resize(ref _productPool, Math.Max(needed, Math.Max(4096, _productPool.Length * 2))); - } - - private void EnsureFragmentCapacity(int needed) - { - if (_batch.FragmentNeutralMasses.Length < needed) - Array.Resize(ref _batch.FragmentNeutralMasses, Math.Max(needed, Math.Max(1024, _batch.FragmentNeutralMasses.Length * 2))); - } - - private void EnsureOffsetCapacity(int needed) - { - if (_batch.PeptideFragmentOffsets.Length < needed) - Array.Resize(ref _batch.PeptideFragmentOffsets, Math.Max(needed, _batch.PeptideFragmentOffsets.Length * 2)); - } - - private void EnsureWorkCapacity(int needed) - { - if (_batch.WorkPeptideSlot.Length < needed) - { - int n = Math.Max(needed, Math.Max(1024, _batch.WorkPeptideSlot.Length * 2)); - Array.Resize(ref _batch.WorkPeptideSlot, n); - Array.Resize(ref _batch.WorkScanIndex, n); - Array.Resize(ref _workNotch, n); - } - } - } } public class TransientSearchEngineResults(