Rudimentary tests and a fix they already found

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;

tests/test_pf9.py

+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