Merge pull request #266 from lbl-anp/edit-peakfinder

update peakfinder to suppress double peaks and update plots

becquerel/core/peakfinder.py

@@ -302,19 +302,20 @@ def add_peak(self, xpeak):
     def plot(self, facecolor="red", linecolor="red", alpha=0.5, peaks=True):
         """Plot the peak signal-to-noise ratios calculated using the kernel."""
         bin_edges = self.spectrum.bin_edges_raw
+        bin_centers = self.spectrum.bin_centers_raw
 
         if facecolor is not None:
-            plt.fill_between(bin_edges[:-1], self.snr, 0, color=facecolor, alpha=alpha)
+            plt.fill_between(bin_centers, self.snr, 0, color=facecolor, alpha=alpha)
         if linecolor is not None:
-            plt.plot(bin_edges[:-1], self.snr, "-", color=linecolor)
+            plt.plot(bin_centers, self.snr, "-", color=linecolor)
         if peaks:
             for cent, snr, fwhm in zip(self.centroids, self.snrs, self.fwhms):
                 plt.plot([cent] * 2, [0, snr], "b-", lw=1.5)
                 plt.plot(cent, snr, "bo")
                 plt.plot(
                     [cent - fwhm / 2, cent + fwhm / 2], [snr / 2] * 2, "b-", lw=1.5
                 )
-        plt.xlim(0, bin_edges.max())
+        plt.xlim(0, bin_centers.max())
         plt.ylim(0)
         plt.xlabel("x")
         plt.ylabel("SNR")
@@ -389,15 +390,9 @@ def find_peaks(self, xmin=None, xmax=None, min_snr=2, max_num=40):
         max_num = int(max_num)
         if max_num < 1:
             raise PeakFinderError("Must keep at least 1 peak, not {}".format(max_num))
-        # calculate the first derivative and second derivatives of the SNR
-        d1 = (self.snr[2:] - self.snr[:-2]) / 2
-        d1 = np.append(0, d1)
-        d1 = np.append(d1, 0)
-        d2 = self.snr[2:] - 2 * self.snr[1:-1] + self.snr[:-2]
-        d2 = np.append(0, d2)
-        d2 = np.append(d2, 0)
+
         # find maxima
-        peak = (d1[2:] < 0) & (d1[:-2] > 0) & (d2[1:-1] < 0)
+        peak = (self.snr[:-2] < self.snr[1:-1]) & (self.snr[1:-1] >= self.snr[2:])
         peak = np.append(False, peak)
         peak = np.append(peak, False)
         # select peaks using SNR and centroid criteria
@@ -415,3 +410,4 @@ def find_peaks(self, xmin=None, xmax=None, min_snr=2, max_num=40):
         self.backgrounds = self.backgrounds[-max_num:]
         # sort by centroid
         self.sort_by(self.centroids)
+        self.peak = peak

tests/autocal_test.py

@@ -87,13 +87,45 @@ def test_peakfinder():
     kernel = bq.GaussianPeakFilter(500, 50, fwhm_at_0=10)
     finder = bq.PeakFinder(spec1, kernel)
     finder.find_peak(500, min_snr=3.0)
-    assert np.isclose(finder.channels[0], 485.5)
+    assert np.isclose(finder.centroids[0], 485.5)
     finder.reset()
     finder.find_peaks()
-    assert len(finder.channels) == 9
+    assert len(finder.centroids) == 9
     finder.reset()
     finder.find_peaks(min_snr=0.5, xmin=50, xmax=1000, max_num=10)
-    assert len(finder.channels) == 10
+    assert len(finder.centroids) == 10
+
+
+def test_peakfinder_double_peaks():
+    """Test that the new peakfinder logic correctly gives only one peak, even
+    when min_sep is not specified, while old logic gives two peaks.
+    """
+
+    # Toy data: gaussian lineshape
+    x = np.arange(31)
+    y = 1.5 * np.exp(-0.5 * (x[:-1] - 15) ** 2 / (2 * 3 ** 2))
+    spec = bq.Spectrum(bin_edges_raw=x, counts=y)
+    kernel = bq.GaussianPeakFilter(15, 7.5, fwhm_at_0=7.5)
+    finder = bq.PeakFinder(spec, kernel)
+
+    # OLD code from PeakFinder.find_peaks
+    d1 = (finder.snr[2:] - finder.snr[:-2]) / 2
+    d1 = np.append(0, d1)
+    d1 = np.append(d1, 0)
+    d2 = finder.snr[2:] - 2 * finder.snr[1:-1] + finder.snr[:-2]
+    d2 = np.append(0, d2)
+    d2 = np.append(d2, 0)
+    peak_old = (d1[2:] < 0) & (d1[:-2] > 0) & (d2[1:-1] < 0)
+    peak_old = np.append(False, peak_old)
+    peak_old = np.append(peak_old, False)
+    assert peak_old.sum() == 2
+
+    # NEW peakfind
+    finder.find_peaks(min_snr=1.0)
+    assert len(finder.centroids) == 1
+    assert finder.centroids[0] in finder.spectrum.bin_centers_raw[peak_old]
+    peak_new = np.isclose(finder.spectrum.bin_centers_raw, finder.centroids[0])
+    assert peak_new.sum() == 1
 
 
 def test_peakfinder_exceptions():
@@ -192,7 +224,7 @@ def test_autocal_spec1():
     finder = bq.PeakFinder(spec1, kernel)
     cal = bq.AutoCalibrator(finder)
     finder.find_peaks(min_snr=1, xmin=50)
-    assert len(cal.peakfinder.channels) == 10
+    assert len(cal.peakfinder.centroids) == 10
     cal.fit(
         REQUIRED,
         optional=OPTIONAL,
@@ -210,7 +242,7 @@ def test_autocal_one_line():
     finder = bq.PeakFinder(spec1, kernel)
     cal = bq.AutoCalibrator(finder)
     finder.find_peaks(min_snr=8, xmin=50)
-    assert len(cal.peakfinder.channels) == 1
+    assert len(cal.peakfinder.centroids) == 1
     cal.fit(
         [1460.82],
         gain_range=[2.5, 4.0],
@@ -230,7 +262,7 @@ def test_autocal_exceptions():
     cal = bq.AutoCalibrator(finder)
     # only one peak found, but multiple energies required
     finder.find_peaks(min_snr=10, xmin=50)
-    assert len(cal.peakfinder.channels) == 1
+    assert len(cal.peakfinder.centroids) == 1
     with pytest.raises(bq.AutoCalibratorError):
         cal.fit(
             REQUIRED,
@@ -240,7 +272,7 @@ def test_autocal_exceptions():
     # multiple peaks found, but only one energy given
     finder.reset()
     finder.find_peaks(min_snr=1, xmin=50)
-    assert len(cal.peakfinder.channels) == 10
+    assert len(cal.peakfinder.centroids) == 10
     with pytest.raises(bq.AutoCalibratorError):
         cal.fit(
             [1460.82],
@@ -250,7 +282,7 @@ def test_autocal_exceptions():
     # more energies required than peaks found
     finder.reset()
     finder.find_peaks(min_snr=4, xmin=50)
-    assert len(cal.peakfinder.channels) == 4
+    assert len(cal.peakfinder.centroids) == 4
     with pytest.raises(bq.AutoCalibratorError):
         cal.fit(
             [238.63, 351.93, 609.32, 1460.82, 2614.3],
@@ -265,7 +297,7 @@ def test_autocal_no_fit():
     finder = bq.PeakFinder(spec1, kernel)
     cal = bq.AutoCalibrator(finder)
     finder.find_peaks(min_snr=2, xmin=50)
-    assert len(cal.peakfinder.channels) == 8
+    assert len(cal.peakfinder.centroids) == 8
     with pytest.raises(bq.AutoCalibratorError):
         cal.fit(
             REQUIRED,
@@ -282,7 +314,7 @@ def test_autocal_plot():
     finder = bq.PeakFinder(spec1, kernel)
     cal = bq.AutoCalibrator(finder)
     finder.find_peaks(min_snr=1, xmin=50)
-    assert len(cal.peakfinder.channels) == 10
+    assert len(cal.peakfinder.centroids) == 10
     cal.fit(
         REQUIRED,
         optional=OPTIONAL,
@@ -299,15 +331,15 @@ def test_autocal_spec2():
     kernel = bq.GaussianPeakFilter(3700, 10, 5)
     finder = bq.PeakFinder(spec2, kernel)
     cal = bq.AutoCalibrator(finder)
-    cal.peakfinder.find_peaks(min_snr=15, xmin=1000)
-    assert len(cal.peakfinder.channels) == 12
+    cal.peakfinder.find_peaks(min_snr=20, xmin=1000)
+    assert len(cal.peakfinder.centroids) == 9
     cal.fit(
         REQUIRED,
         optional=OPTIONAL,
         gain_range=[0.1, 0.6],
         de_max=5.0,
     )
-    assert len(cal.fit_channels) == 6
+    assert len(cal.fit_channels) == 3
     assert np.isclose(cal.gain, 0.3785, rtol=1e-2)
 
 
@@ -317,7 +349,7 @@ def test_autocal_spec3():
     finder = bq.PeakFinder(spec3, kernel)
     cal = bq.AutoCalibrator(finder)
     cal.peakfinder.find_peaks(min_snr=3, xmin=100)
-    assert len(cal.peakfinder.channels) == 8
+    assert len(cal.peakfinder.centroids) == 8
     # this fit succeeds but misidentifies the lines
     cal.fit(
         [609.32, 1460.82],
@@ -344,7 +376,7 @@ def test_autocal_spec4():
     finder = bq.PeakFinder(spec4, kernel)
     cal = bq.AutoCalibrator(finder)
     cal.peakfinder.find_peaks(min_snr=3, xmin=100)
-    assert len(cal.peakfinder.channels) == 7
+    assert len(cal.peakfinder.centroids) == 7
     cal.fit(
         [356.0129, 661.657, 1460.82],
         optional=[911.20, 1120.294, 1620.50, 1764.49, 2118.514, 2614.3],