diff --git a/src/calng/kernels/peakfinder9_gpu.cu b/src/calng/kernels/peakfinder9_gpu.cu index fda8259ea46d47e157d0481621540662dd412b77..fe2cb05c144de4bc91894cfedd25a0112d0d6367 100644 --- a/src/calng/kernels/peakfinder9_gpu.cu +++ b/src/calng/kernels/peakfinder9_gpu.cu @@ -154,8 +154,10 @@ extern "C" __global__ void pf9(const unsigned short num_frames, // whole peak should have sufficent SNR float peak_weighted_ss; float peak_weighted_fs; + // start at 0 because center is in loop as layer 0 + // note: this means min_snr_peak_pixels should be <= min_snr_biggest_pixel float peak_total_mass = 0; - unsigned num_peak_pixels = 1; // will expand from center pixel, count it already + unsigned num_peak_pixels = 0; { float peak_weighted_ss_nom = 0; float peak_weighted_fs_nom = 0; diff --git a/tests/test_pf9.py b/tests/test_pf9.py new file mode 100644 index 0000000000000000000000000000000000000000..235790fe55c01b6856676861490458aace8f3744 --- /dev/null +++ b/tests/test_pf9.py @@ -0,0 +1,131 @@ +import functools + +from calng.correction_addons import peakfinder9 +import cupy as cp +import numpy as np +from numpy.testing import assert_allclose + + +# we don't need to be e-6 or e-7 precise +assert_allclose = functools.partial(assert_allclose, rtol=1e-2, atol=1e-2) + + +class NotAHash(dict): + set = dict.__setitem__ + merge = dict.update + has = dict.__contains__ + + +class Pf9Wrapper: + def __init__(self, **params): + self.output = NotAHash() + default_params = { + "windowRadius": 2, # assumed default + "maxPeaks": 500, + "minPeakValueOverNeighbors": 10, # assumed default + "minSnrMaxPixel": 5, + "minSnrPeakPixels": 4, + "minSnrWholePeak": 6, + "minPeakPixels": 1, + "maxPeakPixels": 500, + "minSigma": 5, + "blockX": 1, + "blockY": 1, + "blockZ": 64, + } + print(default_params | params) + if (unknown_params := params.keys() - default_params.keys()): + raise ValueError(f"Unknown PF9 parameters {unknown_params}") + self.pf9 = peakfinder9.Peakfinder9( + None, # device + "", # prefix + NotAHash(default_params | params), + ) + + def run(self, data): + self.pf9.post_correction( + None, # train id + cp.asarray(data)[np.newaxis], + None, # cell table + None, # pulse table + self.output, # hash + ) + num_peaks = self.output["peakfinding.numPeaks"][0] + return ( + ( + self.output["peakfinding.peakX"][0, :num_peaks], + self.output["peakfinding.peakY"][0, :num_peaks], + ) + ) + + +def test_small_peak(): + pf = Pf9Wrapper() + data = cp.random.random_sample((100, 100), dtype=cp.float32) + + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 0 + + data[40, 60] = 100 + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 1 + assert_allclose(peak_x, 40) + assert_allclose(peak_y, 60) + + data[40, 61] = 95 + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 1 + assert_allclose(peak_x, 40) + assert 60 < peak_y[0] < 61 + + +def test_window_radius(): + pf = Pf9Wrapper(windowRadius=5) + data = cp.random.random_sample((100, 100), dtype=cp.float32) + data[50:55, 60] = 100 + + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 1 + assert 50 <= peak_x[0] <= 55 + assert_allclose(peak_y, 60) + + data[55:60, 60] = 100 + # we cannot pick a valid candidate as none will exceed the window border + peak_x, peak_y = pf.run(data) + print(peak_x, peak_y) + assert peak_x.size == peak_y.size == 0 + + # with bigger window, we can + pf = Pf9Wrapper(windowRadius=20) + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 1 + assert 50 <= peak_x[0] <= 60 + assert_allclose(peak_y, 60) + + # but if middle pixel stands out, it's also fine + pf = Pf9Wrapper(windowRadius=5) + data[55, 60] = 200 + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 1 + assert 50 <= peak_x[0] <= 60 + assert_allclose(peak_y, 60) + + +def test_min_max_peak_pixels(): + data = cp.random.random_sample((100, 100), dtype=cp.float32) + + # radius 5 means we look at 2*5+1 square i.e. up to 121 pixels + pf = Pf9Wrapper(windowRadius=5, minPeakPixels=5, maxPeakPixels=10) + data[20:22, 70:72] = 100 # 4 pixels not enough + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 0 + + data[20:23, 70:73] = 100 # 9 pixels just right + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 1 + assert 20 <= peak_x[0] < 23 + assert 70 <= peak_y[0] < 73 + + data[20:25, 70:75] = 100 # 25 pixels too much + peak_x, peak_y = pf.run(data) + assert peak_x.size == peak_y.size == 0