From 38205ea5349b13520066661e3cd127ccc2791da7 Mon Sep 17 00:00:00 2001 From: Jayson Vavrek Date: Fri, 23 Apr 2021 13:57:28 -0700 Subject: [PATCH 1/6] update peakfinder to suppress double peaks and update plots --- becquerel/core/peakfinder.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/becquerel/core/peakfinder.py b/becquerel/core/peakfinder.py index a4f56286..b0a794cd 100644 --- a/becquerel/core/peakfinder.py +++ b/becquerel/core/peakfinder.py @@ -302,11 +302,12 @@ 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) @@ -314,7 +315,7 @@ def plot(self, facecolor="red", linecolor="red", alpha=0.5, peaks=True): 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[:-2] < self.snr[1:-1]) 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 From e31e8a4783dbc1684f1c99d71a2c5034382f8edf Mon Sep 17 00:00:00 2001 From: Jayson Vavrek Date: Thu, 6 May 2021 11:16:31 -0700 Subject: [PATCH 2/6] add a test of old/new peakfinder logic for double peak issue --- tests/autocal_test.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/tests/autocal_test.py b/tests/autocal_test.py index a6abaac9..a04fe240 100644 --- a/tests/autocal_test.py +++ b/tests/autocal_test.py @@ -99,6 +99,38 @@ def test_peakfinder(): assert len(finder.channels) == 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.) + 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(): """Test PeakFinder exceptions.""" kernel = bq.GaussianPeakFilter(500, 50, fwhm_at_0=10) From 7448007023711494840dc8afbd7954b9c78412c6 Mon Sep 17 00:00:00 2001 From: Jayson Vavrek Date: Thu, 6 May 2021 11:16:44 -0700 Subject: [PATCH 3/6] format --- tests/autocal_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/autocal_test.py b/tests/autocal_test.py index a04fe240..a3c67830 100644 --- a/tests/autocal_test.py +++ b/tests/autocal_test.py @@ -124,7 +124,7 @@ def test_peakfinder_double_peaks(): assert peak_old.sum() == 2 # NEW peakfind - finder.find_peaks(min_snr=1.) + 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]) From bbdadd93df6f85bc43e3871795af29ceb5ce5033 Mon Sep 17 00:00:00 2001 From: Jayson Vavrek Date: Thu, 6 May 2021 12:50:50 -0700 Subject: [PATCH 4/6] adjust a slightly brittle autocal test --- tests/autocal_test.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/autocal_test.py b/tests/autocal_test.py index a3c67830..7211c67e 100644 --- a/tests/autocal_test.py +++ b/tests/autocal_test.py @@ -334,15 +334,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.channels) == 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) From 972bec6ffec10e65f67bc01450890cfc0cccd62e Mon Sep 17 00:00:00 2001 From: Jayson Vavrek Date: Thu, 6 May 2021 12:52:23 -0700 Subject: [PATCH 5/6] eliminate warnings on channels vs centroids --- tests/autocal_test.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/tests/autocal_test.py b/tests/autocal_test.py index 7211c67e..4c5e593b 100644 --- a/tests/autocal_test.py +++ b/tests/autocal_test.py @@ -90,13 +90,13 @@ 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(): @@ -227,7 +227,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, @@ -245,7 +245,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], @@ -265,7 +265,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, @@ -275,7 +275,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], @@ -285,7 +285,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], @@ -300,7 +300,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, @@ -317,7 +317,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, @@ -335,7 +335,7 @@ def test_autocal_spec2(): finder = bq.PeakFinder(spec2, kernel) cal = bq.AutoCalibrator(finder) cal.peakfinder.find_peaks(min_snr=20, xmin=1000) - assert len(cal.peakfinder.channels) == 9 + assert len(cal.peakfinder.centroids) == 9 cal.fit( REQUIRED, optional=OPTIONAL, @@ -352,7 +352,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], @@ -379,7 +379,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], From 9796068259c4a34d102c4e06b48a7c5e1e74485f Mon Sep 17 00:00:00 2001 From: Jayson Vavrek Date: Thu, 6 May 2021 13:32:34 -0700 Subject: [PATCH 6/6] small order change in peakfinder --- becquerel/core/peakfinder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/becquerel/core/peakfinder.py b/becquerel/core/peakfinder.py index b0a794cd..a7f95dd4 100644 --- a/becquerel/core/peakfinder.py +++ b/becquerel/core/peakfinder.py @@ -392,7 +392,7 @@ def find_peaks(self, xmin=None, xmax=None, min_snr=2, max_num=40): raise PeakFinderError("Must keep at least 1 peak, not {}".format(max_num)) # find maxima - peak = (self.snr[2:] <= self.snr[1:-1]) & (self.snr[:-2] < self.snr[1:-1]) + 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