Add isotopically-unresolved (average-mass) envelope support

mzLib/MassSpectrometry/Deconvolution/AverageResidue/AverageResidue.cs

@@ -16,6 +16,42 @@ public abstract class AverageResidue : IEquatable<AverageResidue>
     public abstract double[] GetAllTheoreticalIntensities(int index);
     public abstract double GetDiffToMonoisotopic(int index);
 
+    // Terminology for the offsets composed from this model:
+    //  • most-abundant offset = GetDiffToMonoisotopic(GetMostIntenseMassIndex(mono)) — the gap from the
+    //    monoisotopic mass to the single tallest (most abundant) isotopologue. Used by the resolved
+    //    most-abundant search; composed from the abstract methods above, so no dedicated method is needed.
+    //  • average / centroid offset = GetAverageOffset(mono) below — the gap from the monoisotopic mass to
+    //    the intensity-weighted mean of the whole envelope. Used for isotopically unresolved (high-mass)
+    //    species, where the observed precursor is a centroid rather than a resolved most-abundant peak.
+
+    /// <summary>
+    /// The number of daltons to add to a monoisotopic mass to reach the intensity-weighted average
+    /// (centroid) mass of the isotopic envelope, according to this average-residue model. Used for
+    /// isotopically unresolved (high-mass) species, where the observed precursor is a centroid rather
+    /// than a resolved most-abundant peak.
+    /// </summary>
+    /// <param name="monoisotopicMass">A monoisotopic mass, in daltons.</param>
+    public double GetAverageOffset(double monoisotopicMass)
+    {
+        int index = GetMostIntenseMassIndex(monoisotopicMass);
+        double[] masses = GetAllTheoreticalMasses(index);
+        double[] intensities = GetAllTheoreticalIntensities(index);
+
+        double totalIntensity = 0;
+        double weightedMassSum = 0;
+        for (int i = 0; i < masses.Length; i++)
+        {
+            totalIntensity += intensities[i];
+            weightedMassSum += masses[i] * intensities[i];
+        }
+        double averageMass = totalIntensity > 0 ? weightedMassSum / totalIntensity : masses[0];
+
+        // masses/intensities are sorted by intensity (most intense first); the model's monoisotopic
+        // mass is the most-intense mass minus its diff-to-monoisotopic.
+        double monoisotopicModelMass = masses[0] - GetDiffToMonoisotopic(index);
+        return averageMass - monoisotopicModelMass;
+    }
+
     #region IEquatable<AverageResidue>
 
     public bool Equals(AverageResidue? other)

mzLib/MassSpectrometry/MzSpectra/IsotopicEnvelope.cs

@@ -6,15 +6,95 @@
 
 namespace MassSpectrometry
 {
+    /// <summary>
+    /// Whether the isotopic peaks of an envelope were individually resolved by deconvolution.
+    /// A <b>resolved</b> envelope has a well-defined most-abundant isotopic peak, so candidate
+    /// selection uses <see cref="IsotopicEnvelope.MostAbundantObservedNeutralMass"/>. An <b>unresolved</b>
+    /// envelope (typically high-mass, &gt;~35 kDa) has no individually resolved peaks, so its
+    /// intensity-weighted <see cref="IsotopicEnvelope.AverageObservedMass"/> (centroid) is used instead.
+    /// </summary>
+    public enum EnvelopeResolution
+    {
+        Resolved,
+        Unresolved
+    }
+
     public class IsotopicEnvelope : IHasMass, IEquatable<IsotopicEnvelope>
     {
         public readonly List<(double mz, double intensity)> Peaks;
         public double MonoisotopicMass { get; private set; }
 
         /// <summary>
-        /// Mass of most abundant observed isotopic peak, not accounting for addition or subtraction or protons due to ESI charge state induction
+        /// Most abundant observed isotopic peak as m/z × |charge|, i.e. <b>not</b> proton-corrected to a
+        /// neutral mass (it does not subtract the charge-carrier protons added during ESI). This is a
+        /// charge-scaled m/z, not a neutral mass; for the proton-corrected neutral mass use
+        /// <see cref="MostAbundantObservedNeutralMass"/>. Retained for existing deconvolution-quality tests.
+        ///
+        /// <b>Sentinel:</b> this value is only meaningful for an envelope that carries observed isotopic
+        /// peaks. Constructors that have no observed envelope (e.g. a neutral mass read from a
+        /// pre-deconvoluted file) set it to <c>-1</c> to mark "no most-abundant peak available", rather
+        /// than a misleading zero or a synthetic value.
+        /// </summary>
+        public double MostAbundantObservedIsotopicMass { get; private set; }
+
+        /// <summary>
+        /// The <b>most-abundant observed neutral mass</b>: the neutral mass of the single most intense
+        /// (tallest) observed isotopic peak of this charge-state envelope — the highest-signal, directly
+        /// measured point of the envelope. The most intense peak's m/z is converted to a neutral mass via
+        /// the envelope charge (i.e. proton-corrected, <c>mz.ToMass(Charge)</c>), so it is directly
+        /// comparable to a theoretical proteoform neutral mass. It is the precursor mass used for candidate
+        /// selection in most-abundant mode. Computed from <see cref="Peaks"/> and <see cref="Charge"/>
+        /// (both readonly).
+        ///
+        /// Terminology (the three masses are distinct):
+        ///  • <b>most-abundant (observed neutral) mass</b> — this property; the tallest single isotopic
+        ///    peak, proton-corrected to a neutral mass.
+        ///  • <b>monoisotopic mass</b> (<see cref="MonoisotopicMass"/>) — the all-light-isotope mass; for
+        ///    large proteoforms this isotopologue is rare and often undetectable.
+        ///  • <b>average / centroid mass</b> — the intensity-weighted mean over the whole envelope, used
+        ///    for isotopically unresolved (high-mass) species. That property is added separately, with the
+        ///    unresolved-envelope work, and is not part of this resolved most-abundant feature.
+        /// This is the neutral-mass form of <see cref="MostAbundantObservedIsotopicMass"/> (which is the
+        /// SAME tallest peak expressed as m/z × |charge|, NOT proton-corrected): the two differ by exactly
+        /// <c>|charge| × ProtonMass</c>. Deriving it from that single source means the two values can never
+        /// disagree, and it carries the same <c>-1</c> "no observed peak" sentinel.
         /// </summary>
-        internal double MostAbundantObservedIsotopicMass { get; private set; }
+        public double MostAbundantObservedNeutralMass =>
+            MostAbundantObservedIsotopicMass < 0 ? -1 : MostAbundantObservedIsotopicMass - Charge * Constants.ProtonMass;
+
+        /// <summary>
+        /// The <b>average (centroid) mass</b>: the intensity-weighted mean neutral mass over all observed
+        /// peaks of the envelope (proton-corrected). This is the mass-spectrometry "average mass" — the
+        /// abundance-weighted counterpart to the monoisotopic mass — and is distinct from the
+        /// <see cref="MostAbundantObservedNeutralMass"/> (the single tallest peak). It is the precursor mass used
+        /// for candidate selection when the envelope is isotopically <see cref="EnvelopeResolution.Unresolved"/>
+        /// (high-mass species with no individually resolved most-abundant peak). Shares the same -1
+        /// "no observed envelope" sentinel as <see cref="MostAbundantObservedIsotopicMass"/>.
+        /// </summary>
+        public double AverageObservedMass
+        {
+            get
+            {
+                if (MostAbundantObservedIsotopicMass < 0 || Peaks is not { Count: > 0 })
+                    return -1;
+
+                double totalIntensity = Peaks.Sum(p => p.intensity);
+                double centroidMz = totalIntensity > 0
+                    ? Peaks.Sum(p => p.mz * p.intensity) / totalIntensity
+                    : Peaks.MaxBy(p => p.intensity).mz;
+                return centroidMz.ToMass(Charge);
+            }
+        }
+
+        /// <summary>
+        /// Whether deconvolution individually resolved the isotopic peaks of this envelope. Defaults to
+        /// <see cref="EnvelopeResolution.Resolved"/>; deconvolution algorithms that detect unresolved
+        /// (high-mass) envelopes set this to <see cref="EnvelopeResolution.Unresolved"/> via
+        /// <see cref="SetResolution"/>, which switches candidate selection from the most-abundant peak to
+        /// the average (centroid) mass.
+        /// </summary>
+        public EnvelopeResolution Resolution { get; private set; } = EnvelopeResolution.Resolved;
+
         public readonly int Charge;
         public readonly double TotalIntensity;
         public readonly int PrecursorId;
@@ -60,6 +140,9 @@ public IsotopicEnvelope(double monoisotopicMass, double intensity, int charge)
             TotalIntensity = intensity;
             Score = double.MaxValue;
             Peaks = [(monoisotopicMass.ToMz(charge), intensity)];
+            // A neutral mass read from a pre-deconvoluted file has no observed isotopic envelope, so
+            // there is no most-abundant observed peak to report — use the -1 sentinel.
+            MostAbundantObservedIsotopicMass = -1;
         }
 
         /// <summary>
@@ -82,6 +165,17 @@ public IsotopicEnvelope(int id, List<(double mz, double intensity)> peaks, doubl
             MostAbundantObservedIsotopicMass = peaks.MaxBy(p => p.intensity).mz * Math.Abs(chargestate);
         }
 
+        /// <summary>
+        /// Sets the resolution state of this envelope. Deconvolution algorithms call this when they
+        /// determine an envelope's isotopic peaks were not individually resolved (e.g. high-mass
+        /// species), switching candidate selection to the average (centroid) mass. Defaults to
+        /// <see cref="EnvelopeResolution.Resolved"/>.
+        /// </summary>
+        public void SetResolution(EnvelopeResolution resolution)
+        {
+            Resolution = resolution;
+        }
+
         /// <summary>
         /// Sets the optional generic deconvolution score for this envelope. The generic score is
         /// produced by <see cref="DeconvolutionScorer"/> and is comparable across deconvolution
@@ -140,6 +234,7 @@ public bool Equals(IsotopicEnvelope other)
             if (other is null) return false;
             if (ReferenceEquals(this, other)) return true;
             if (Charge != other.Charge || Peaks.Count != other.Peaks.Count) return false;
+            if (Resolution != other.Resolution) return false;
             if (Math.Abs(TotalIntensity - other.TotalIntensity) >= 0.001) return false;
             if (Math.Abs(MonoisotopicMass - other.MonoisotopicMass) >= 0.001) return false;
             if (Math.Abs(MostAbundantObservedIsotopicMass - other.MostAbundantObservedIsotopicMass) >= 0.001) return false;
@@ -176,6 +271,7 @@ public override int GetHashCode()
                 hash = hash * 23 + TotalIntensity.GetHashCode();
                 hash = hash * 23 + MonoisotopicMass.GetHashCode();
                 hash = hash * 23 + MostAbundantObservedIsotopicMass.GetHashCode();
+                hash = hash * 23 + Resolution.GetHashCode();
                 return hash;
             }
         }

mzLib/Test/Deconvolution/TestMostAbundantMass.cs

@@ -0,0 +1,149 @@
+using System.Collections.Generic;
+using System.Linq;
+using Chemistry;
+using MassSpectrometry;
+using NUnit.Framework;
+
+namespace Test
+{
+    /// <summary>
+    /// Unit tests for the resolved "most abundant mass" precursor-selection support (Strategy B).
+    ///
+    /// Terminology pinned by these tests:
+    ///  • most-abundant mass = the neutral mass of the single most intense (tallest) isotopic peak,
+    ///    proton-corrected (<see cref="IsotopicEnvelope.MostAbundantObservedNeutralMass"/>);
+    ///  • most-abundant offset = GetDiffToMonoisotopic(GetMostIntenseMassIndex(mono)) on the averagine
+    ///    model — the gap from the monoisotopic mass to that tallest isotopologue.
+    /// The intensity-weighted average (centroid) mass and the isotopically-unresolved path are a
+    /// separate change and are tested with that work, not here. All tests build synthetic envelopes
+    /// from the Averagine model — no deconvolution is run.
+    /// </summary>
+    [TestFixture]
+    public sealed class TestMostAbundantMass
+    {
+        private static readonly AverageResidue Model = new Averagine();
+
+        // Most-abundant offset for a monoisotopic mass = the averagine diff-to-monoisotopic of the
+        // nearest mass bin (the mass-keyed composition consumers use in place of a dedicated method).
+        private static double MostAbundantOffset(double monoMass) => Model.GetDiffToMonoisotopic(Model.GetMostIntenseMassIndex(monoMass));
+
+        /// <summary>
+        /// Builds a perfect synthetic envelope: peaks at exact theoretical m/z with
+        /// Averagine-proportional intensities. (Same construction as TestDeconvolutionScorerUnit.)
+        /// </summary>
+        private static List<(double mz, double intensity)> BuildPerfectPeaks(double monoMass, int charge, double baseIntens = 1e6)
+        {
+            int avgIdx = Model.GetMostIntenseMassIndex(monoMass);
+            double[] rawMasses = Model.GetAllTheoreticalMasses(avgIdx);
+            double[] rawIntens = Model.GetAllTheoreticalIntensities(avgIdx);
+
+            var sorted = rawMasses.Zip(rawIntens).OrderBy(p => p.First).ToArray();
+
+            double isotopeStep = Constants.C13MinusC12 / charge;
+            double monoMz = monoMass.ToMz(charge);
+            var peaks = new List<(double mz, double intensity)>();
+            for (int n = 0; n < sorted.Length; n++)
+            {
+                double intensity = baseIntens * sorted[n].Second;
+                if (intensity < baseIntens * 0.001) continue;
+                peaks.Add((monoMz + n * isotopeStep, intensity));
+            }
+            return peaks;
+        }
+
+        private static IsotopicEnvelope BuildPerfectEnvelope(double monoMass, int charge, double baseIntens = 1e6)
+        {
+            var peaks = BuildPerfectPeaks(monoMass, charge, baseIntens);
+            return new IsotopicEnvelope(0, peaks, monoMass, charge, peaks.Sum(p => p.intensity), 0.999);
+        }
+
+        // ── AverageResidue most-abundant offset ─────────────────────────────────
+
+        [Test]
+        public void MostAbundantOffset_IsNonNegativeAndNonDecreasingWithMass()
+        {
+            double[] masses = { 500, 2000, 5000, 10000, 20000, 40000 };
+            double prev = double.NegativeInfinity;
+            foreach (double m in masses)
+            {
+                double offset = MostAbundantOffset(m);
+                Assert.That(offset, Is.GreaterThanOrEqualTo(-1e-6)); // ~0 at tiny mass (mono IS most abundant)
+                Assert.That(offset, Is.GreaterThanOrEqualTo(prev - 1e-6), $"offset decreased at mass {m}");
+                prev = offset;
+            }
+        }
+
+        [Test]
+        public void MostAbundantOffset_IsNearZeroForSmallMass()
+        {
+            // A small peptide's monoisotopic peak is (nearly) the most abundant.
+            Assert.That(MostAbundantOffset(500), Is.LessThan(0.5));
+        }
+
+        [Test]
+        public void MostAbundantOffset_GrowsRoughlyOneNeutronPer1600Da()
+        {
+            // ~1 13C neutron (~1.00235 Da) per ~1.6 kDa. Assert the offset at 16 kDa is in a
+            // physically plausible band (≈ 9–11 Da) rather than an exact value.
+            double offset = MostAbundantOffset(16000);
+            Assert.That(offset, Is.GreaterThan(8.0).And.LessThan(12.0));
+        }
+
+        // ── IsotopicEnvelope most-abundant observed mass ────────────────────────
+
+        [Test]
+        public void MostAbundantObservedNeutralMass_IsProtonCorrectedNeutralMass()
+        {
+            const int charge = 10;
+            var env = BuildPerfectEnvelope(15000, charge);
+
+            // The proton-corrected neutral mass equals the most intense peak's m/z .ToMass(charge)...
+            double mostIntenseMz = env.Peaks.MaxBy(p => p.intensity).mz;
+            Assert.That(env.MostAbundantObservedNeutralMass, Is.EqualTo(mostIntenseMz.ToMass(charge)).Within(1e-6));
+
+            // ...and it differs from the un-proton-corrected mz*|z| field by exactly z proton masses.
+            Assert.That(env.MostAbundantObservedIsotopicMass - env.MostAbundantObservedNeutralMass,
+                Is.EqualTo(charge * Constants.ProtonMass).Within(1e-6));
+        }
+
+        [Test]
+        public void MostAbundantObservedNeutralMass_MatchesMonoPlusAveragineOffset()
+        {
+            // Ties the two features together: the observed most-abundant neutral mass of a perfect
+            // envelope ≈ candidate monoisotopic + averagine most-abundant offset.
+            const double mono = 12000;
+            const int charge = 12;
+            var env = BuildPerfectEnvelope(mono, charge);
+
+            double predicted = mono + MostAbundantOffset(mono);
+            Assert.That(env.MostAbundantObservedNeutralMass, Is.EqualTo(predicted).Within(0.15));
+        }
+
+        [Test]
+        public void DeconvolutionConstructor_ComputesMostAbundantObservedNeutralMass()
+        {
+            // The 5-arg mzLib-deconvolution constructor must compute the most-abundant observed mass.
+            const double mono = 12000;
+            const int charge = 12;
+            var peaks = BuildPerfectPeaks(mono, charge);
+            var env = new IsotopicEnvelope(peaks, mono, charge, peaks.Sum(p => p.intensity), 0.5);
+
+            double mostIntenseMz = peaks.MaxBy(p => p.intensity).mz;
+            Assert.That(env.MostAbundantObservedNeutralMass, Is.EqualTo(mostIntenseMz.ToMass(charge)).Within(1e-6));
+        }
+
+        [Test]
+        public void FileReadEnvelope_HasNoMostAbundantPeak_ReturnsSentinel()
+        {
+            // The file-read constructor carries a neutral mass but no observed isotopic envelope, so the
+            // most-abundant observed mass is undefined: both the m/z×|charge| form and the proton-corrected
+            // form report the -1 sentinel rather than a synthetic value.
+            const double mono = 8000;
+            const int charge = 8;
+            var env = new IsotopicEnvelope(mono, 1e6, charge);
+
+            Assert.That(env.MostAbundantObservedIsotopicMass, Is.EqualTo(-1));
+            Assert.That(env.MostAbundantObservedNeutralMass, Is.EqualTo(-1));
+        }
+    }
+}