Merge branch '242_low_spectrum_edge_warning' into 'dev'

Fix resolving power at spectrum edges without warning spam

See merge request mass-spectrometry/corems!258

Author: kheal

corems/encapsulation/factory/processingSetting.py

@@ -585,7 +585,7 @@ class MassSpectrumSetting:
         Maximum m/z to use for peak picking. Default is 1200.0.
     picking_point_extrapolate : int, optional
         How many data points (in each direction) to extrapolate the mz axis and 0 pad the abundance axis. Default is 3.
-        Recommend 3 for reduced profile data or if peak picking faults
+        Recommend not changing this value unless necessary.
     calib_minimize_method : str, optional
         Minimization method to use for calibration. Default is 'Powell'.
     calib_pol_order : int, optional

corems/mass_spectrum/calc/PeakPicking.py

@@ -396,7 +396,9 @@ def calculate_resolving_power(self, intes, massa, current_index):
         Returns
         -------
         float
-            The resolving power of the peak.
+            The resolving power of the peak, or ``nan`` if a half-maximum
+            crossing does not exist on both sides (spectrum edge truncates
+            the peak).
 
         Notes
         --------
@@ -408,37 +410,19 @@ def calculate_resolving_power(self, intes, massa, current_index):
 
         peak_height = intes[current_index]
         target_peak_height = peak_height / 2
+        n_points = len(intes)
 
-        peak_height_minus = peak_height
-        peak_height_plus = peak_height
-
-        # There are issues when a peak is at the high or low limit of a spectrum in finding its local minima and maxima
-        # This solution will return nan for resolving power when a peak is possibly too close to an edge to avoid the issue
-
-        if current_index < 5:
-            warnings.warn("peak at low spectrum edge, returning no resolving power")
-            return nan
-        elif abs(current_index - len(intes)) < 5:
-            warnings.warn("peak at high spectrum edge, returning no resolving power")
+        index_minus = current_index
+        while index_minus > 0 and intes[index_minus] >= target_peak_height:
+            index_minus -= 1
+        if intes[index_minus] >= target_peak_height:
             return nan
-        else:
-            pass
 
-        index_minus = current_index
-        while peak_height_minus >= target_peak_height:
-            index_minus = index_minus - 1
-            if index_minus < 0:
-                warnings.warn(
-                    "Res. calc. warning - peak index minus adjacent to spectrum edge \n \
-                        Zeroing the first 5 data points of abundance. Peaks at spectrum edge may be incorrectly reported \n \
-                        Perhaps try to increase picking_point_extrapolate (e.g. to 3)"
-                )
-                # Pad the first 5 data points with zeros and restart the loop
-                intes[:5] = 0
-                peak_height_minus = target_peak_height
-                index_minus = current_index
-            else:
-                peak_height_minus = intes[index_minus]
+        index_plus = current_index
+        while index_plus < n_points - 1 and intes[index_plus] >= target_peak_height:
+            index_plus += 1
+        if intes[index_plus] >= target_peak_height:
+            return nan
 
         if self.mspeaks_settings.legacy_centroid_polyfit:
             x = [massa[index_minus], massa[index_minus + 1]]
@@ -459,23 +443,6 @@ def calculate_resolving_power(self, intes, massa, current_index):
             y_intercept = target_peak_height
         massa1 = (y_intercept - b) / a
 
-        index_plus = current_index
-        while peak_height_plus >= target_peak_height:
-            index_plus = index_plus + 1
-
-            try:
-                peak_height_plus = intes[index_plus]
-            except IndexError:
-                warnings.warn(
-                    "Res. calc. warning - peak index plus adjacent to spectrum edge \n \
-                        Zeroing the last 5 data points of abundance. Peaks at spectrum edge may be incorrectly reported\
-                        Perhaps try to increase picking_point_extrapolate (e.g. to 3)"
-                )
-                # Pad the first 5 data points with zeros and restart the loop
-                intes[-5:] = 0
-                peak_height_plus = target_peak_height
-                index_plus = current_index
-
         if self.mspeaks_settings.legacy_centroid_polyfit:
             x = [massa[index_plus], massa[index_plus - 1]]
             y = [intes[index_plus], intes[index_plus - 1]]

tests/test_mass_spectrum.py

@@ -1,3 +1,10 @@
+import warnings
+from types import SimpleNamespace
+
+from numpy import allclose, arange, array, copy, isfinite, isnan
+
+from corems.mass_spectrum.calc.PeakPicking import PeakPicking
+
 # Tests for adpodization methods
 def test_hamming(bruker_transient):
     """Test the creation of a mass spectrum object with the Hamming apodization method"""
@@ -103,4 +110,98 @@ def test_mass_spectrum_properties(mass_spectrum_ftms):
     res = mass_spectrum_ftms.get_mz_and_abundance_peaks_tuples()
     assert len(res) == len(mass_spectrum_ftms)
     mass_spectrum_ftms.get_masses_count_by_nominal_mass()
-    mass_spectrum_ftms.resolving_power_calc(12, 1)
+    mass_spectrum_ftms.resolving_power_calc(12, 1)
+
+
+def _resolving_power_picker():
+    picker = PeakPicking()
+    picker.mspeaks_settings = SimpleNamespace(
+        legacy_resolving_power=True,
+        legacy_centroid_polyfit=False,
+    )
+    return picker
+
+
+def _spectrum_edge_warnings(caught):
+    needles = (
+        "peak at low spectrum edge",
+        "peak at high spectrum edge",
+        "Zeroing the first 5 data points",
+        "Zeroing the last 5 data points",
+        "peak index minus adjacent to spectrum edge",
+        "peak index plus adjacent to spectrum edge",
+    )
+    return [
+        str(w.message)
+        for w in caught
+        if any(n in str(w.message) for n in needles)
+    ]
+
+
+def test_resolving_power_when_fwhm_is_measurable():
+    """Interior and padded-edge peaks get finite RP; abundance is not mutated."""
+    picker = _resolving_power_picker()
+    cases = (
+        (
+            array([100.0, 100.1, 100.2, 100.3, 100.4]),
+            array([0.0, 50.0, 100.0, 50.0, 0.0]),
+            2,
+            334.0,
+        ),
+        (
+            100.0 + 0.1 * arange(10),
+            array([0.0, 0.0, 0.0, 100.0, 50.0, 0.0, 10.0, 10.0, 10.0, 10.0]),
+            3,
+            None,
+        ),
+        (
+            100.0 + 0.1 * arange(10),
+            array([10.0, 10.0, 10.0, 10.0, 0.0, 50.0, 100.0, 50.0, 0.0, 0.0]),
+            6,
+            None,
+        ),
+    )
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        for massa, intes, apex, expected in cases:
+            original = copy(intes)
+            rp = picker.calculate_resolving_power(intes, massa, apex)
+            assert isfinite(rp)
+            assert rp > 0
+            if expected is not None:
+                assert allclose(rp, expected)
+            assert allclose(intes, original)
+    assert _spectrum_edge_warnings(caught) == []
+
+
+def test_resolving_power_nan_when_fwhm_truncated():
+    """Missing half-max on either side yields silent nan and does not zero abundance."""
+    picker = _resolving_power_picker()
+    cases = (
+        (
+            100.0 + 0.1 * arange(8),
+            array([100.0, 50.0, 0.0, 10.0, 10.0, 10.0, 10.0, 10.0]),
+            0,
+        ),
+        (
+            100.0 + 0.1 * arange(8),
+            array([10.0, 10.0, 10.0, 10.0, 10.0, 0.0, 50.0, 100.0]),
+            7,
+        ),
+        (
+            100.0 + 0.1 * arange(12),
+            array(
+                [80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 100.0, 50.0, 0.0, 10.0, 10.0, 10.0]
+            ),
+            6,
+        ),
+    )
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        for massa, intes, apex in cases:
+            original = copy(intes)
+            rp = picker.calculate_resolving_power(intes, massa, apex)
+            assert isnan(rp)
+            assert allclose(intes, original)
+    assert _spectrum_edge_warnings(caught) == []
+