From 1bf317f42a4fdd2deb9d5d5a029ffd6f7d71c248 Mon Sep 17 00:00:00 2001 From: kferjaoui Date: Mon, 3 Aug 2026 11:39:02 +0200 Subject: [PATCH] CUDA: double-precision device path, local-max fix, pedestal introspection - DEVICE_PED_TYPE for device pedestal/variance (shipped double/double to match the double CPU ClusterFinder). - Local-max suppression in Test 1/Test 3: non-peak pixels no longer store or update, mirroring ClusterFinder's `value < max -> continue`. - device_pedestal()/device_noise() accessors (+ bindings) for the in-kernel decision-time pedestal. --- include/aare/ClusterFinderCUDA.hpp | 96 ++++++++++++++++++++++----- include/aare/clusterfinder_kernel.cuh | 58 ++++++++++------ python/src/bind_ClusterFinderCUDA.hpp | 23 +++++++ 3 files changed, 142 insertions(+), 35 deletions(-) diff --git a/include/aare/ClusterFinderCUDA.hpp b/include/aare/ClusterFinderCUDA.hpp index bf8fa3d1..98c2a7f8 100644 --- a/include/aare/ClusterFinderCUDA.hpp +++ b/include/aare/ClusterFinderCUDA.hpp @@ -20,9 +20,10 @@ template struct StreamContext { cudaStream_t stream = nullptr; FRAME_TYPE *d_frame = nullptr; - float *d_pd_mean = nullptr; // always float on device; host stays double - float *d_pd_sum = nullptr; - float *d_pd_sum2 = nullptr; + // Device pedestal precision is set by DEVICE_PED_TYPE in the kernel header. + device::DEVICE_PED_TYPE *d_pd_mean = nullptr; + device::DEVICE_PED_TYPE *d_pd_sum = nullptr; + device::DEVICE_PED_TYPE *d_pd_sum2 = nullptr; uint8_t *d_output = nullptr; // [uint32_t count | ClusterType clusters[max]] }; @@ -170,9 +171,12 @@ class ClusterFinderCUDA { CUDA_CHECK( cudaStreamCreateWithFlags(&sc.stream, cudaStreamNonBlocking)); CUDA_CHECK(cudaMalloc(&sc.d_frame, m_image_bytes)); - CUDA_CHECK(cudaMalloc(&sc.d_pd_mean, m_image_size * sizeof(float))); - CUDA_CHECK(cudaMalloc(&sc.d_pd_sum, m_image_size * sizeof(float))); - CUDA_CHECK(cudaMalloc(&sc.d_pd_sum2, m_image_size * sizeof(float))); + CUDA_CHECK(cudaMalloc( + &sc.d_pd_mean, m_image_size * sizeof(device::DEVICE_PED_TYPE))); + CUDA_CHECK(cudaMalloc( + &sc.d_pd_sum, m_image_size * sizeof(device::DEVICE_PED_TYPE))); + CUDA_CHECK(cudaMalloc( + &sc.d_pd_sum2, m_image_size * sizeof(device::DEVICE_PED_TYPE))); CUDA_CHECK(cudaMalloc(&sc.d_output, m_output_bytes_per_frame)); } @@ -266,6 +270,67 @@ class ClusterFinderCUDA { NDArray pedestal() { return m_pedestal.mean(); } NDArray noise() { return m_pedestal.std(); } + /** + * @brief Device pedestal MEAN for one stream — the pedestal the kernel + * actually reads and updates in place every frame. This differs from + * pedestal() (the host pedestal, advanced only by + * push_pedestal_frame): it carries the in-kernel running update, so + * it is the baseline an accept/reject decision was really made + * against. Reading it right BEFORE a find_clusters() call gives the + * state that call will decide with (the kernel updates at frame + * end). With the single-frame path every frame lands on stream 0. + */ + NDArray device_pedestal(int stream = 0) { + if (m_pedestal_dirty) { + sync_pedestal_to_device(); + m_pedestal_dirty = false; + } + auto &sc = v_sc.at(static_cast(stream)); + CUDA_CHECK(cudaStreamSynchronize(sc.stream)); + using DPT = device::DEVICE_PED_TYPE; + std::vector h_mean(m_image_size); + CUDA_CHECK(cudaMemcpy(h_mean.data(), sc.d_pd_mean, + m_image_size * sizeof(DPT), + cudaMemcpyDeviceToHost)); + NDArray out( + {static_cast(nrows), static_cast(ncols)}); + for (size_t i = 0; i < m_image_size; ++i) + out.data()[i] = static_cast(h_mean[i]); + return out; + } + + /** + * @brief Device pedestal RMS for one stream, computed exactly as the kernel + * does: sqrt(max(sum2/n - mean^2, 0)). See device_pedestal(). + */ + NDArray device_noise(int stream = 0) { + if (m_pedestal_dirty) { + sync_pedestal_to_device(); + m_pedestal_dirty = false; + } + auto &sc = v_sc.at(static_cast(stream)); + CUDA_CHECK(cudaStreamSynchronize(sc.stream)); + using DPT = device::DEVICE_PED_TYPE; + std::vector h_mean(m_image_size), h_sum2(m_image_size); + CUDA_CHECK(cudaMemcpy(h_mean.data(), sc.d_pd_mean, + m_image_size * sizeof(DPT), + cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(h_sum2.data(), sc.d_pd_sum2, + m_image_size * sizeof(DPT), + cudaMemcpyDeviceToHost)); + const double n = static_cast(m_pedestal.n_samples()); + NDArray out( + {static_cast(nrows), static_cast(ncols)}); + for (size_t i = 0; i < m_image_size; ++i) { + double var = + static_cast(h_sum2[i]) / n - + static_cast(h_mean[i]) * static_cast(h_mean[i]); + out.data()[i] = + static_cast(std::sqrt(std::max(var, 0.0))); + } + return out; + } + /** * @brief Move clusters out of the internal ClusterVector, optionally * reallocating the internal one with the same capacity. @@ -580,18 +645,19 @@ class ClusterFinderCUDA { NDArray h_sum = m_pedestal.get_sum(); NDArray h_sum2 = m_pedestal.get_sum2(); - // Host accumulates in double for precision; cast to float for device - // to halve global-memory bandwidth and eliminate FP64 arithmetic. - std::vector f_mean(m_image_size); - std::vector f_sum(m_image_size); - std::vector f_sum2(m_image_size); + // Host accumulates in double; cast to the device pedestal precision + // (DEVICE_PED_TYPE — float or double per the kernel-header toggle). + using DPT = device::DEVICE_PED_TYPE; + std::vector f_mean(m_image_size); + std::vector f_sum(m_image_size); + std::vector f_sum2(m_image_size); for (size_t i = 0; i < m_image_size; ++i) { - f_mean[i] = static_cast(h_mean.data()[i]); - f_sum[i] = static_cast(h_sum.data()[i]); - f_sum2[i] = static_cast(h_sum2.data()[i]); + f_mean[i] = static_cast(h_mean.data()[i]); + f_sum[i] = static_cast(h_sum.data()[i]); + f_sum2[i] = static_cast(h_sum2.data()[i]); } - const size_t bytes = m_image_size * sizeof(float); + const size_t bytes = m_image_size * sizeof(DPT); for (auto &sc : v_sc) { CUDA_CHECK(cudaMemcpyAsync(sc.d_pd_mean, f_mean.data(), bytes, cudaMemcpyHostToDevice, sc.stream)); diff --git a/include/aare/clusterfinder_kernel.cuh b/include/aare/clusterfinder_kernel.cuh index 2686bf72..d57b7db1 100644 --- a/include/aare/clusterfinder_kernel.cuh +++ b/include/aare/clusterfinder_kernel.cuh @@ -6,17 +6,25 @@ namespace aare::device { -// Implementing mixed precision for shared memory and stencil arithmetic -using COMPUTE_TYPE = float; +// Device arithmetic precision. +// COMPUTE_TYPE : per-frame stencil arithmetic (shared-memory tile, sums). +// DEVICE_PED_TYPE : device pedestal storage + running variance update. +// Shipped as double/double so the GPU result matches the double-precision CPU +// ClusterFinder to the floating-point floor. float/float is ~2x faster but +// reintroduces a small near-threshold mismatch; a build-time toggle for that +// trade-off is planned. +using COMPUTE_TYPE = double; +using DEVICE_PED_TYPE = double; template , typename FRAME_TYPE = uint16_t, typename = std::enable_if_t::value>> __global__ void find_clusters_in_single_frame( - const FRAME_TYPE *__restrict__ d_frame, float *__restrict__ d_pd_mean, - float *__restrict__ d_pd_sum, float *__restrict__ d_pd_sum2, - const uint32_t n_pd_samples, const COMPUTE_TYPE m_nSigma, - const size_t nrows, const size_t ncols, + const FRAME_TYPE *__restrict__ d_frame, + DEVICE_PED_TYPE *__restrict__ d_pd_mean, + DEVICE_PED_TYPE *__restrict__ d_pd_sum, + DEVICE_PED_TYPE *__restrict__ d_pd_sum2, const uint32_t n_pd_samples, + const COMPUTE_TYPE m_nSigma, const size_t nrows, const size_t ncols, // const uint64_t frame_number, ClusterType *d_clusters, uint32_t *d_cluster_count, const uint32_t max_clusters) { @@ -174,11 +182,13 @@ __global__ void find_clusters_in_single_frame( // Per-pixel variance from global pedestal arrays // Variance = rms^2 = E[X^2] - E[X]^2 // NOTE: Keep thresholds squared to avoid one sqrtf() per pixel. - float mean_px = d_pd_mean[global_tid]; - float var_px = d_pd_sum2[global_tid] / static_cast(n_pd_samples) - - mean_px * mean_px; - float rms_sq = fmaxf(var_px, 0.0f); - float nSig_sq_rms_sq = m_nSigma * m_nSigma * rms_sq; + DEVICE_PED_TYPE mean_px = d_pd_mean[global_tid]; + DEVICE_PED_TYPE var_px = + d_pd_sum2[global_tid] / static_cast(n_pd_samples) - + mean_px * mean_px; + COMPUTE_TYPE rms_sq = static_cast( + var_px > DEVICE_PED_TYPE{0} ? var_px : DEVICE_PED_TYPE{0}); + COMPUTE_TYPE nSig_sq_rms_sq = m_nSigma * m_nSigma * rms_sq; // Pedestal-subtracted value of the center pixel (already in shmem) COMPUTE_TYPE val_pixel = shmem[shmem_tid]; @@ -193,7 +203,7 @@ __global__ void find_clusters_in_single_frame( // Stencil reduction: total, max, quadrant sums COMPUTE_TYPE total = 0.0f; - COMPUTE_TYPE max_val = -HUGE_VALF; + COMPUTE_TYPE max_val = -HUGE_VAL; // double inf; narrows to float if needed // // Quandrants // PEDESTAL_TYPE tl = PEDESTAL_TYPE{0}; // top-left quadrant (ir<=0, @@ -209,7 +219,7 @@ __global__ void find_clusters_in_single_frame( COMPUTE_TYPE val = shmem[shmem_tid + ir * shmem_stride + ic]; total += val; - max_val = fmaxf(max_val, val); + max_val = val > max_val ? val : max_val; // // Quadrant accumulation (pixels on the axes contribute to two // quadrants) if (ir <= 0 && ic <= 0) tl += val; if (ir <= 0 && ic @@ -255,8 +265,15 @@ __global__ void find_clusters_in_single_frame( // Test 3: total significance (only if tests 1 & 2 didn't fire) if (!is_photon) { - if (total > 0.0f && - total * total > static_cast(pow2_c3) * nSig_sq_rms_sq) { + if (total > 0.0f && total * total > static_cast(pow2_c3) * + nSig_sq_rms_sq) { + // Local-max suppression, mirroring ClusterFinder's `value == max` + // store gate: only the peak pixel of the window records the + // cluster, so an extended charge-shared event yields one cluster, + // not one per pixel. Return (not fall-through) because the serial + // total branch pushes no pedestal for these non-max pixels. + if (val_pixel < max_val) + return; is_photon = true; } } @@ -267,10 +284,11 @@ __global__ void find_clusters_in_single_frame( // frame simultaneously. So the updated pedestal will only be used starting // from the next frame. -> This avoids a/serialization and b/global mem I/O. if (!is_photon && valid_pixel) { - float raw_val = static_cast(d_frame[global_tid]); - float sum = d_pd_sum[global_tid]; - float sum2 = d_pd_sum2[global_tid]; - float n = static_cast(n_pd_samples); + DEVICE_PED_TYPE raw_val = + static_cast(d_frame[global_tid]); + DEVICE_PED_TYPE sum = d_pd_sum[global_tid]; + DEVICE_PED_TYPE sum2 = d_pd_sum2[global_tid]; + DEVICE_PED_TYPE n = static_cast(n_pd_samples); sum += raw_val - sum / n; sum2 += raw_val * raw_val - sum2 / n; @@ -297,7 +315,7 @@ __global__ void find_clusters_in_single_frame( for (int ic = -col_radius; ic <= col_radius; ++ic) { COMPUTE_TYPE val = shmem[shmem_tid + ir * shmem_stride + ic]; if constexpr (std::is_integral_v) - clusterData[idx] = static_cast(lroundf(val)); + clusterData[idx] = static_cast(lround(val)); else clusterData[idx] = static_cast(val); idx++; diff --git a/python/src/bind_ClusterFinderCUDA.hpp b/python/src/bind_ClusterFinderCUDA.hpp index fc6f160d..2d555394 100644 --- a/python/src/bind_ClusterFinderCUDA.hpp +++ b/python/src/bind_ClusterFinderCUDA.hpp @@ -66,6 +66,29 @@ void define_ClusterFinderCUDA(py::module &m, const std::string &typestr) { return return_image_data(arr); }) + .def( + "device_pedestal", + [](CF &self, int stream) { + auto pd = new NDArray{}; + *pd = self.device_pedestal(stream); + return return_image_data(pd); + }, + py::arg("stream") = 0, + R"(Device pedestal MEAN for a stream — the pedestal the kernel +actually decides with and updates each frame, unlike `pedestal` (the frozen +host pedestal). Read it before find_clusters() for the decision-time state.)") + + .def( + "device_noise", + [](CF &self, int stream) { + auto arr = new NDArray{}; + *arr = self.device_noise(stream); + return return_image_data(arr); + }, + py::arg("stream") = 0, + R"(Device pedestal RMS for a stream, computed as the kernel does: +sqrt(max(sum2/n - mean^2, 0)). Counterpart to `noise` for the device pedestal.)") + .def( "steal_clusters", [](CF &self, bool realloc_same_capacity) {