detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-09-03 04:20:43 +02:00
Code Issues Actions Packages Releases Activity

Refactor ClusterFinderCUDA
Build on RHEL9 / build (push) Successful in 3m22s
Details
Build on RHEL8 / build (push) Successful in 3m28s
Details
Run tests using data on local RHEL8 / build (push) Successful in 3m37s
Details

Browse Source 
Rework the multi-stream pipeline to eliminate per-frame sync barriers and
fix the D2H staging architecture.

Sync reduction:
- Replace one cudaStreamSynchronize per frame with one per stream per batch,
  cutting synchronisation calls from O(n_frames x n_streams) to O(n_streams)
- Introduce a unified per-frame D2H output layout [uint32_t count | clusters[max]]
  stored in a single class-level lazy-allocated pinned pool (h_output_pinned),
  replacing the per-stream separate cluster/count device buffers
- Move CUDA event pool from per-stream fixed-size to per-frame-slot lazy-allocated,
  enabling correct kernel timing across any batch size

Pinned H2D without CPU-side copy:
- Add register_input_buffer(ptr, bytes) / unregister_input_buffer() wrapping
  cudaHostRegister so callers can pin their existing batch buffer once; all
  find_clusters_batched() slices then transfer at DMA speed (~22 GB/s) instead
  of ~15 GB/s for pageable, with no extra memcpy or WC-memory penalty

Result (RTX 4090, 400x400 uint16, 3x3 clusters, batch=2000, 5 streams):
  Before: ~34 µs/frame  ->  After: ~28 µs/frame  (−18 %)
This commit is contained in:
kferjaoui
2026-05-18 16:30:13 +02:00
parent 5cfc888567
commit 6a12e3de24
8 changed files with 1497 additions and 290 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/aare/ClusterFinderCUDA.hpp
+183 -193
View File
@@ -7,31 +7,21 @@
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <limits>
#include <stdexcept>
namespace aare {
// Per-stream device resources
// Per-stream device resources (device-side only; all pinned host staging is
// class-level)
template <typename ClusterType, typename FRAME_TYPE, typename PEDESTAL_TYPE>
struct StreamContext {
cudaStream_t stream = nullptr; // handle to the stream
cudaStream_t stream = nullptr;
FRAME_TYPE *d_frame = nullptr;
PEDESTAL_TYPE *d_pd_mean = nullptr;
PEDESTAL_TYPE *d_pd_sum = nullptr;
PEDESTAL_TYPE *d_pd_sum2 = nullptr;
ClusterType *d_clusters = nullptr;
uint32_t *d_cluster_count = nullptr;
// Pinned host staging buffers. These make cudaMemcpyAsync real async DMA
// transfers even when the caller's NDView points to pageable memory.
FRAME_TYPE *h_frame = nullptr;
uint32_t *h_cluster_count = nullptr;
ClusterType *h_clusters = nullptr;
cudaEvent_t kernel_start = nullptr;
cudaEvent_t kernel_stop = nullptr;
uint8_t *d_output = nullptr; // [uint32_t count | ClusterType clusters[max]]
};
template <typename ClusterType = Cluster<int32_t, 3, 3>,
@@ -46,15 +36,28 @@ class ClusterFinderCUDA {
static constexpr int col_radius = ClusterType::cluster_size_x / 2;
static constexpr int row_radius = ClusterType::cluster_size_y / 2;
size_t m_output_pinned_capacity =
0; // # frames currently allocated in h_output_pinned
void *h_output_pinned = nullptr;
// Pointer registered via cudaHostRegister by the caller (for pinned H2D
// speed). The class tracks it only to unregister in the destructor if the
// caller forgets.
void *m_registered_input = nullptr;
Shape<2> m_shape;
size_t nrows;
size_t ncols;
size_t m_image_size; // nrows * ncols
int n_streams;
size_t m_capacity;
size_t m_image_bytes;
size_t m_cluster_bytes;
int n_streams;
size_t m_max_clusters_per_frame;
// Per-frame output layout helpers
size_t m_output_bytes_per_frame; // sizeof(uint32_t) + max *
// sizeof(ClusterType), aligned
size_t m_clusters_offset; // offset of cluster array within output block
COMPUTE_TYPE m_nSigma;
Pedestal<PEDESTAL_TYPE> m_pedestal;
@@ -67,6 +70,11 @@ class ClusterFinderCUDA {
float m_total_kernel_ms = 0.0f;
size_t m_frames_processed = 0;
// Per-frame kernel timing event pool. Sized lazily to the largest batch
// seen; one event pair per frame slot.
std::vector<cudaEvent_t> m_kernel_start_pool;
std::vector<cudaEvent_t> m_kernel_stop_pool;
// Kernel parameters
dim3 grid;
dim3 block;
@@ -79,30 +87,32 @@ class ClusterFinderCUDA {
* @param m_image_size shape of the detector frame (rows, cols)
* @param nSigma threshold in units of per-pixel pedestal
* std
* @param capacity device-side cluster buffer size per stream
* @param max_clusters_per_frame tight upper bound on clusters/frame for
* fixed-size D2H
* @param n_streams number of CUDA streams for multi-frame
* overlap
*/
ClusterFinderCUDA(Shape<2> shape_, COMPUTE_TYPE nSigma = 5.0,
size_t capacity = 1000000, int n_streams_ = 1)
size_t max_clusters_per_frame = 2048, int n_streams_ = 5)
: m_shape(shape_), nrows(shape_[0]), ncols(shape_[1]),
m_image_size(nrows * ncols), n_streams(n_streams_),
m_capacity(capacity), m_nSigma(nSigma),
m_pedestal(shape_[0], shape_[1]), m_clusters(capacity) {
m_max_clusters_per_frame(max_clusters_per_frame), m_nSigma(nSigma),
m_pedestal(shape_[0], shape_[1]), m_clusters(max_clusters_per_frame) {
if (n_streams_ <= 0) {
throw std::invalid_argument(
"ClusterFinderCUDA: n_streams must be > 0");
}
if (capacity >
if (max_clusters_per_frame >
static_cast<size_t>(std::numeric_limits<uint32_t>::max())) {
throw std::invalid_argument(
"ClusterFinderCUDA: capacity must fit in uint32_t");
"ClusterFinderCUDA: max_clusters_per_frame must fit in "
"uint32_t");
}
if (capacity == 0) {
if (max_clusters_per_frame == 0) {
throw std::invalid_argument(
"ClusterFinderCUDA: capacity must be > 0");
"ClusterFinderCUDA: max_clusters_per_frame must be > 0");
}
// Grid/Block dimensions
@@ -118,15 +128,22 @@ class ClusterFinderCUDA {
sizeof(COMPUTE_TYPE);
m_image_bytes = m_image_size * sizeof(FRAME_TYPE);
m_cluster_bytes = m_capacity * sizeof(ClusterType);
// Output block layout: [count][padding to ClusterType
// alignment][clusters]
constexpr size_t cluster_align = alignof(ClusterType);
const size_t count_bytes = sizeof(uint32_t);
// next multiple of cluster_align
m_clusters_offset =
(count_bytes + cluster_align - 1) & ~(cluster_align - 1);
m_output_bytes_per_frame =
m_clusters_offset + m_max_clusters_per_frame * sizeof(ClusterType);
v_sc.resize(n_streams);
for (int k = 0; k < n_streams; ++k) {
auto &sc = v_sc[k];
CUDA_CHECK(
cudaStreamCreateWithFlags(&sc.stream, cudaStreamNonBlocking));
CUDA_CHECK(cudaEventCreate(&sc.kernel_start));
CUDA_CHECK(cudaEventCreate(&sc.kernel_stop));
CUDA_CHECK(cudaMalloc(&sc.d_frame, m_image_bytes));
CUDA_CHECK(cudaMalloc(&sc.d_pd_mean,
m_image_size * sizeof(PEDESTAL_TYPE)));
@@ -134,19 +151,7 @@ class ClusterFinderCUDA {
cudaMalloc(&sc.d_pd_sum, m_image_size * sizeof(PEDESTAL_TYPE)));
CUDA_CHECK(cudaMalloc(&sc.d_pd_sum2,
m_image_size * sizeof(PEDESTAL_TYPE)));
CUDA_CHECK(cudaMalloc(&sc.d_clusters, m_cluster_bytes));
CUDA_CHECK(cudaMalloc(&sc.d_cluster_count, sizeof(uint32_t)));
CUDA_CHECK(cudaMallocHost(reinterpret_cast<void **>(&sc.h_frame),
m_image_bytes));
CUDA_CHECK(
cudaMallocHost(reinterpret_cast<void **>(&sc.h_cluster_count),
sizeof(uint32_t)));
if (m_cluster_bytes > 0) {
CUDA_CHECK(
cudaMallocHost(reinterpret_cast<void **>(&sc.h_clusters),
m_cluster_bytes));
}
CUDA_CHECK(cudaMalloc(&sc.d_output, m_output_bytes_per_frame));
}
}
@@ -154,7 +159,6 @@ class ClusterFinderCUDA {
for (auto &sc : v_sc) {
if (sc.stream)
cudaStreamSynchronize(sc.stream);
if (sc.d_frame)
cudaFree(sc.d_frame);
if (sc.d_pd_mean)
@@ -163,25 +167,26 @@ class ClusterFinderCUDA {
cudaFree(sc.d_pd_sum);
if (sc.d_pd_sum2)
cudaFree(sc.d_pd_sum2);
if (sc.d_clusters)
cudaFree(sc.d_clusters);
if (sc.d_cluster_count)
cudaFree(sc.d_cluster_count);
if (sc.h_frame)
cudaFreeHost(sc.h_frame);
if (sc.h_clusters)
cudaFreeHost(sc.h_clusters);
if (sc.h_cluster_count)
cudaFreeHost(sc.h_cluster_count);
if (sc.kernel_start)
cudaEventDestroy(sc.kernel_start);
if (sc.kernel_stop)
cudaEventDestroy(sc.kernel_stop);
if (sc.d_output)
cudaFree(sc.d_output);
if (sc.stream)
cudaStreamDestroy(sc.stream);
}
for (auto e : m_kernel_start_pool)
if (e)
cudaEventDestroy(e);
for (auto e : m_kernel_stop_pool)
if (e)
cudaEventDestroy(e);
m_kernel_start_pool.clear();
m_kernel_stop_pool.clear();
// free pinned memory
if (h_output_pinned)
cudaFreeHost(h_output_pinned);
if (m_registered_input)
cudaHostUnregister(m_registered_input);
}
// Non-copyable, non-movable
@@ -193,6 +198,30 @@ class ClusterFinderCUDA {
void set_nSigma(COMPUTE_TYPE nSigma) { m_nSigma = nSigma; }
COMPUTE_TYPE get_nSigma() const { return m_nSigma; }
/**
* @brief Pin an existing host buffer so that find_clusters_batched
* transfers it at full PCIe bandwidth (~22 GB/s) instead of going through
* the CUDA driver's internal staging (~15 GB/s for pageable memory).
*
* Call once before the processing loop (not per-frame). The buffer must
* cover the largest NDView you will pass to find_clusters_batched.
* Call unregister_input_buffer() when done, or the destructor will clean
* up.
*/
void register_input_buffer(void *ptr, size_t bytes) {
if (m_registered_input)
CUDA_CHECK(cudaHostUnregister(m_registered_input));
CUDA_CHECK(cudaHostRegister(ptr, bytes, cudaHostRegisterDefault));
m_registered_input = ptr;
}
void unregister_input_buffer() {
if (m_registered_input) {
CUDA_CHECK(cudaHostUnregister(m_registered_input));
m_registered_input = nullptr;
}
}
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);
m_pedestal_dirty = true;
@@ -221,61 +250,20 @@ class ClusterFinderCUDA {
}
/**
* @brief Find clusters in a single frame, appending them to the internal
* ClusterVector.
* @brief Find clusters in a single frame, appending results to the internal
* ClusterVector (accessible via steal_clusters).
* Delegates to find_clusters_batched to avoid duplicating the GPU
* pipeline.
*/
void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
if (m_pedestal_dirty) { // need to update the pedestal on the gpu
sync_pedestal_to_device();
m_pedestal_dirty = false;
}
auto &sc = v_sc[0];
const uint32_t n_pd_samples =
static_cast<uint32_t>(m_pedestal.n_samples());
// First, CPU copies frame into a reusable pinned buffer
std::memcpy(sc.h_frame, frame.data(), m_image_bytes);
// Reset cluster counter
CUDA_CHECK(cudaMemsetAsync(sc.d_cluster_count, 0, sizeof(uint32_t),
sc.stream));
// Upload frame
CUDA_CHECK(cudaMemcpyAsync(sc.d_frame, sc.h_frame, m_image_bytes,
cudaMemcpyHostToDevice, sc.stream));
// Timed Kernel launch
CUDA_CHECK(cudaEventRecord(sc.kernel_start, sc.stream));
device::find_clusters_in_single_frame<ClusterType, FRAME_TYPE,
PEDESTAL_TYPE>
<<<grid, block, shmem_bytes, sc.stream>>>(
sc.d_frame, sc.d_pd_mean, sc.d_pd_sum, sc.d_pd_sum2,
n_pd_samples, m_nSigma, nrows, ncols, sc.d_clusters,
sc.d_cluster_count, static_cast<uint32_t>(m_capacity));
CUDA_CHECK(cudaEventRecord(sc.kernel_stop, sc.stream));
CUDA_CHECK(cudaGetLastError());
// Read back cluster count into pinned buffer
CUDA_CHECK(cudaMemcpyAsync(sc.h_cluster_count, sc.d_cluster_count,
sizeof(uint32_t), cudaMemcpyDeviceToHost,
sc.stream));
// Synchronize to ensure count is available before the CPU reads
// clusters
CUDA_CHECK(cudaStreamSynchronize(sc.stream));
record_kernel_time(sc);
// Clamp to max in case of overflow
uint32_t n_found = *sc.h_cluster_count;
n_found = std::min(n_found, static_cast<uint32_t>(m_capacity));
// Read back clusters
NDView<FRAME_TYPE, 3> batch(
frame.data(),
{1, static_cast<ssize_t>(nrows), static_cast<ssize_t>(ncols)});
auto results = find_clusters_batched(batch, frame_number);
m_clusters.set_frame_number(frame_number);
if (n_found > 0) {
append_device_clusters_to(m_clusters, sc, n_found);
}
auto &cv = results[0];
for (size_t i = 0; i < cv.size(); ++i)
m_clusters.push_back(cv[i]);
}
/**
@@ -293,80 +281,94 @@ class ClusterFinderCUDA {
m_pedestal_dirty = false;
}
const size_t n_frames = frames.shape(0);
const size_t n_frames_batch = frames.shape(0);
const uint32_t n_pd_samples =
static_cast<uint32_t>(m_pedestal.n_samples());
// Lazy grow D2H output staging buffer (one slot per frame)
if (n_frames_batch > m_output_pinned_capacity) {
if (h_output_pinned)
CUDA_CHECK(cudaFreeHost(h_output_pinned));
CUDA_CHECK(cudaMallocHost(
&h_output_pinned, n_frames_batch * m_output_bytes_per_frame));
m_output_pinned_capacity = n_frames_batch;
}
ensure_event_pool(n_frames_batch);
std::vector<ClusterVector<ClusterType>> results;
results.reserve(n_frames);
for (size_t i = 0; i < n_frames; ++i) {
results.reserve(n_frames_batch);
for (size_t i = 0; i < n_frames_batch; ++i) {
results.emplace_back();
results.back().set_frame_number(first_frame + i);
}
const size_t n_rounds = (n_frames + n_streams - 1) / n_streams;
for (size_t round = 0; round < n_rounds; ++round) {
// Launch all frames round-robin across streams.
// If the caller has called register_input_buffer() on frames.data(),
// H2D runs at pinned DMA bandwidth (~22 GB/s); otherwise the CUDA
// driver stages it internally (~15 GB/s for pageable memory).
for (size_t frame_idx = 0; frame_idx < n_frames_batch; ++frame_idx) {
auto &sc = v_sc[frame_idx % n_streams];
// Launch phase: fan out kernels on all streams for this round
for (int k = 0; k < n_streams; ++k) {
// OOB guard
const size_t frame_idx = round * n_streams + k;
if (frame_idx >= n_frames)
continue;
const FRAME_TYPE *h_src = frames.data() + frame_idx * m_image_size;
auto *d_cluster_count = reinterpret_cast<uint32_t *>(sc.d_output);
auto &sc_k = v_sc[k];
const FRAME_TYPE *h_src =
frames.data() + frame_idx * m_image_size;
CUDA_CHECK(cudaMemsetAsync(d_cluster_count, 0, sizeof(uint32_t),
sc.stream));
CUDA_CHECK(cudaMemcpyAsync(sc.d_frame, h_src, m_image_bytes,
cudaMemcpyHostToDevice, sc.stream));
std::memcpy(sc_k.h_frame, h_src, m_image_bytes);
auto *d_clusters = reinterpret_cast<ClusterType *>(
sc.d_output + m_clusters_offset);
CUDA_CHECK(
cudaEventRecord(m_kernel_start_pool[frame_idx], sc.stream));
device::find_clusters_in_single_frame<ClusterType, FRAME_TYPE,
PEDESTAL_TYPE>
<<<grid, block, shmem_bytes, sc.stream>>>(
sc.d_frame, sc.d_pd_mean, sc.d_pd_sum, sc.d_pd_sum2,
n_pd_samples, m_nSigma, nrows, ncols, d_clusters,
d_cluster_count,
static_cast<uint32_t>(m_max_clusters_per_frame));
CUDA_CHECK(
cudaEventRecord(m_kernel_stop_pool[frame_idx], sc.stream));
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaMemsetAsync(sc_k.d_cluster_count, 0,
sizeof(uint32_t), sc_k.stream));
CUDA_CHECK(
cudaMemcpyAsync(sc_k.d_frame, sc_k.h_frame, m_image_bytes,
cudaMemcpyHostToDevice, sc_k.stream));
CUDA_CHECK(cudaEventRecord(sc_k.kernel_start, sc_k.stream));
device::find_clusters_in_single_frame<ClusterType, FRAME_TYPE,
PEDESTAL_TYPE>
<<<grid, block, shmem_bytes, sc_k.stream>>>(
sc_k.d_frame, sc_k.d_pd_mean, sc_k.d_pd_sum,
sc_k.d_pd_sum2, n_pd_samples, m_nSigma, nrows, ncols,
sc_k.d_clusters, sc_k.d_cluster_count,
static_cast<uint32_t>(m_capacity));
CUDA_CHECK(cudaEventRecord(sc_k.kernel_stop, sc_k.stream));
CUDA_CHECK(cudaGetLastError());
// Queue count D2H immediately after the kernel
CUDA_CHECK(cudaMemcpyAsync(
sc_k.h_cluster_count, sc_k.d_cluster_count,
sizeof(uint32_t), cudaMemcpyDeviceToHost, sc_k.stream));
}
// Drain phase: fan in results from all streams
for (int k = 0; k < n_streams; ++k) {
const size_t frame_idx = round * n_streams + k;
if (frame_idx >= n_frames)
continue;
auto &sc_k = v_sc[k];
// Wait for memset -> H2D -> kernel -> count D2H
CUDA_CHECK(cudaStreamSynchronize(sc_k.stream));
record_kernel_time(sc_k);
uint32_t n_found = *sc_k.h_cluster_count;
n_found = std::min<uint32_t>(n_found,
static_cast<uint32_t>(m_capacity));
if (n_found > 0) {
append_device_clusters_to(results[frame_idx], sc_k,
n_found);
}
}
void *h_slot = static_cast<char *>(h_output_pinned) +
frame_idx * m_output_bytes_per_frame;
CUDA_CHECK(cudaMemcpyAsync(h_slot, sc.d_output,
m_output_bytes_per_frame,
cudaMemcpyDeviceToHost, sc.stream));
}
// Sync once per stream
const int streams_used =
std::min<int>(n_streams, static_cast<int>(n_frames_batch));
for (int k = 0; k < streams_used; ++k)
CUDA_CHECK(cudaStreamSynchronize(v_sc[k].stream));
// Drain: fan in results from pinned D2H output pool
for (size_t frame_idx = 0; frame_idx < n_frames_batch; ++frame_idx) {
const void *h_slot = static_cast<const char *>(h_output_pinned) +
frame_idx * m_output_bytes_per_frame;
uint32_t n_found = *reinterpret_cast<const uint32_t *>(h_slot);
n_found = std::min<uint32_t>(
n_found, static_cast<uint32_t>(m_max_clusters_per_frame));
if (n_found > 0) {
const auto *src = reinterpret_cast<const ClusterType *>(
static_cast<const char *>(h_slot) + m_clusters_offset);
for (uint32_t i = 0; i < n_found; ++i)
results[frame_idx].push_back(src[i]);
}
float kernel_ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&kernel_ms,
m_kernel_start_pool[frame_idx],
m_kernel_stop_pool[frame_idx]));
m_total_kernel_ms += kernel_ms;
}
m_frames_processed += n_frames_batch;
return results;
}
@@ -407,28 +409,16 @@ class ClusterFinderCUDA {
CUDA_CHECK(cudaStreamSynchronize(sc.stream));
}
/**
* Copy n_found clusters from sc.d_clusters into the given ClusterVector
* and block on the transfer.
*/
void append_device_clusters_to(ClusterVector<ClusterType> &cv, SC &sc,
uint32_t n_found) {
CUDA_CHECK(cudaMemcpyAsync(sc.h_clusters, sc.d_clusters,
n_found * sizeof(ClusterType),
cudaMemcpyDeviceToHost, sc.stream));
CUDA_CHECK(cudaStreamSynchronize(sc.stream));
for (uint32_t i = 0; i < n_found; ++i)
cv.push_back(sc.h_clusters[i]);
}
void record_kernel_time(SC &sc) {
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, sc.kernel_start, sc.kernel_stop));
m_total_kernel_ms += ms;
m_frames_processed++;
void ensure_event_pool(size_t n_frames) {
const size_t old_size = m_kernel_start_pool.size();
if (n_frames <= old_size)
return;
m_kernel_start_pool.resize(n_frames);
m_kernel_stop_pool.resize(n_frames);
for (size_t i = old_size; i < n_frames; ++i) {
CUDA_CHECK(cudaEventCreate(&m_kernel_start_pool[i]));
CUDA_CHECK(cudaEventCreate(&m_kernel_stop_pool[i]));
}
}
};
include/aare/ClusterFinderCUDA_old.hpp
+435
View File
@@ -0,0 +1,435 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/ClusterFinder.hpp"
#include "aare/clusterfinder_kernel.cuh"
#include "aare/utils/cuda_check.cuh"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <limits>
#include <stdexcept>
namespace aare {
// Per-stream device resources
template <typename ClusterType, typename FRAME_TYPE, typename PEDESTAL_TYPE>
struct StreamContext {
cudaStream_t stream = nullptr; // handle to the stream
FRAME_TYPE *d_frame = nullptr;
PEDESTAL_TYPE *d_pd_mean = nullptr;
PEDESTAL_TYPE *d_pd_sum = nullptr;
PEDESTAL_TYPE *d_pd_sum2 = nullptr;
ClusterType *d_clusters = nullptr;
uint32_t *d_cluster_count = nullptr;
// Pinned host staging buffers. These make cudaMemcpyAsync real async DMA
// transfers even when the caller's NDView points to pageable memory.
FRAME_TYPE *h_frame = nullptr;
uint32_t *h_cluster_count = nullptr;
ClusterType *h_clusters = nullptr;
cudaEvent_t kernel_start = nullptr;
cudaEvent_t kernel_stop = nullptr;
};
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinderCUDA {
using COMPUTE_TYPE =
device::COMPUTE_TYPE; // match the kernel's internal precision
static constexpr int BLOCK_X = 16;
static constexpr int BLOCK_Y = 16;
static constexpr int col_radius = ClusterType::cluster_size_x / 2;
static constexpr int row_radius = ClusterType::cluster_size_y / 2;
Shape<2> m_shape;
size_t nrows;
size_t ncols;
size_t m_image_size; // nrows * ncols
int n_streams;
size_t m_capacity;
size_t m_image_bytes;
size_t m_cluster_bytes;
COMPUTE_TYPE m_nSigma;
Pedestal<PEDESTAL_TYPE> m_pedestal;
ClusterVector<ClusterType> m_clusters;
bool m_pedestal_dirty = true;
using SC = StreamContext<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>;
std::vector<SC> v_sc;
float m_total_kernel_ms = 0.0f;
size_t m_frames_processed = 0;
// Kernel parameters
dim3 grid;
dim3 block;
size_t shmem_bytes;
public:
/**
* @brief Construct a ClusterFinderCUDA
*
* @param m_image_size shape of the detector frame (rows, cols)
* @param nSigma threshold in units of per-pixel pedestal
* std
* @param capacity device-side cluster buffer size per stream
* @param n_streams number of CUDA streams for multi-frame
* overlap
*/
ClusterFinderCUDA(Shape<2> shape_, COMPUTE_TYPE nSigma = 5.0,
size_t capacity = 1000000, int n_streams_ = 1)
: m_shape(shape_), nrows(shape_[0]), ncols(shape_[1]),
m_image_size(nrows * ncols), n_streams(n_streams_),
m_capacity(capacity), m_nSigma(nSigma),
m_pedestal(shape_[0], shape_[1]), m_clusters(capacity) {
if (n_streams_ <= 0) {
throw std::invalid_argument(
"ClusterFinderCUDA: n_streams must be > 0");
}
if (capacity >
static_cast<size_t>(std::numeric_limits<uint32_t>::max())) {
throw std::invalid_argument(
"ClusterFinderCUDA: capacity must fit in uint32_t");
}
if (capacity == 0) {
throw std::invalid_argument(
"ClusterFinderCUDA: capacity must be > 0");
}
// Grid/Block dimensions
block = dim3(BLOCK_X, BLOCK_Y);
grid = dim3((static_cast<unsigned int>(ncols) + BLOCK_X - 1) / BLOCK_X,
(static_cast<unsigned int>(nrows) + BLOCK_Y - 1) / BLOCK_Y);
// Shared memory: one tile of (BLOCK_X + 2*col_radius) x (BLOCK_Y +
// 2*row_radius) elements
// Mixed precision used -> shmem takes COMPUTE_TYPE = floats (not
// PEDESTAL_TYPE)
shmem_bytes = (BLOCK_X + 2 * col_radius) * (BLOCK_Y + 2 * row_radius) *
sizeof(COMPUTE_TYPE);
m_image_bytes = m_image_size * sizeof(FRAME_TYPE);
m_cluster_bytes = m_capacity * sizeof(ClusterType);
v_sc.resize(n_streams);
for (int k = 0; k < n_streams; ++k) {
auto &sc = v_sc[k];
CUDA_CHECK(
cudaStreamCreateWithFlags(&sc.stream, cudaStreamNonBlocking));
CUDA_CHECK(cudaEventCreate(&sc.kernel_start));
CUDA_CHECK(cudaEventCreate(&sc.kernel_stop));
CUDA_CHECK(cudaMalloc(&sc.d_frame, m_image_bytes));
CUDA_CHECK(cudaMalloc(&sc.d_pd_mean,
m_image_size * sizeof(PEDESTAL_TYPE)));
CUDA_CHECK(
cudaMalloc(&sc.d_pd_sum, m_image_size * sizeof(PEDESTAL_TYPE)));
CUDA_CHECK(cudaMalloc(&sc.d_pd_sum2,
m_image_size * sizeof(PEDESTAL_TYPE)));
CUDA_CHECK(cudaMalloc(&sc.d_clusters, m_cluster_bytes));
CUDA_CHECK(cudaMalloc(&sc.d_cluster_count, sizeof(uint32_t)));
CUDA_CHECK(cudaMallocHost(reinterpret_cast<void **>(&sc.h_frame),
m_image_bytes));
CUDA_CHECK(
cudaMallocHost(reinterpret_cast<void **>(&sc.h_cluster_count),
sizeof(uint32_t)));
if (m_cluster_bytes > 0) {
CUDA_CHECK(
cudaMallocHost(reinterpret_cast<void **>(&sc.h_clusters),
m_cluster_bytes));
}
}
}
~ClusterFinderCUDA() {
for (auto &sc : v_sc) {
if (sc.stream)
cudaStreamSynchronize(sc.stream);
if (sc.d_frame)
cudaFree(sc.d_frame);
if (sc.d_pd_mean)
cudaFree(sc.d_pd_mean);
if (sc.d_pd_sum)
cudaFree(sc.d_pd_sum);
if (sc.d_pd_sum2)
cudaFree(sc.d_pd_sum2);
if (sc.d_clusters)
cudaFree(sc.d_clusters);
if (sc.d_cluster_count)
cudaFree(sc.d_cluster_count);
if (sc.h_frame)
cudaFreeHost(sc.h_frame);
if (sc.h_clusters)
cudaFreeHost(sc.h_clusters);
if (sc.h_cluster_count)
cudaFreeHost(sc.h_cluster_count);
if (sc.kernel_start)
cudaEventDestroy(sc.kernel_start);
if (sc.kernel_stop)
cudaEventDestroy(sc.kernel_stop);
if (sc.stream)
cudaStreamDestroy(sc.stream);
}
}
// Non-copyable, non-movable
ClusterFinderCUDA(const ClusterFinderCUDA &) = delete;
ClusterFinderCUDA &operator=(const ClusterFinderCUDA &) = delete;
ClusterFinderCUDA(ClusterFinderCUDA &&) = delete;
ClusterFinderCUDA &operator=(ClusterFinderCUDA &&) = delete;
void set_nSigma(COMPUTE_TYPE nSigma) { m_nSigma = nSigma; }
COMPUTE_TYPE get_nSigma() const { return m_nSigma; }
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);
m_pedestal_dirty = true;
}
void clear_pedestal() {
m_pedestal.clear();
m_pedestal_dirty = true;
}
NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
/**
* @brief Move clusters out of the internal ClusterVector, optionally
* reallocating the internal one with the same capacity.
*/
ClusterVector<ClusterType>
steal_clusters(bool realloc_same_capacity = false) {
ClusterVector<ClusterType> tmp = std::move(m_clusters);
if (realloc_same_capacity)
m_clusters = ClusterVector<ClusterType>(tmp.capacity());
else
m_clusters = ClusterVector<ClusterType>{};
return tmp;
}
/**
* @brief Find clusters in a single frame, appending them to the internal
* ClusterVector.
*/
void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
if (m_pedestal_dirty) { // need to update the pedestal on the gpu
sync_pedestal_to_device();
m_pedestal_dirty = false;
}
auto &sc = v_sc[0];
const uint32_t n_pd_samples =
static_cast<uint32_t>(m_pedestal.n_samples());
// First, CPU copies frame into a reusable pinned buffer
std::memcpy(sc.h_frame, frame.data(), m_image_bytes);
// Reset cluster counter
CUDA_CHECK(cudaMemsetAsync(sc.d_cluster_count, 0, sizeof(uint32_t),
sc.stream));
// Upload frame
CUDA_CHECK(cudaMemcpyAsync(sc.d_frame, sc.h_frame, m_image_bytes,
cudaMemcpyHostToDevice, sc.stream));
// Timed Kernel launch
CUDA_CHECK(cudaEventRecord(sc.kernel_start, sc.stream));
device::find_clusters_in_single_frame<ClusterType, FRAME_TYPE,
PEDESTAL_TYPE>
<<<grid, block, shmem_bytes, sc.stream>>>(
sc.d_frame, sc.d_pd_mean, sc.d_pd_sum, sc.d_pd_sum2,
n_pd_samples, m_nSigma, nrows, ncols, sc.d_clusters,
sc.d_cluster_count, static_cast<uint32_t>(m_capacity));
CUDA_CHECK(cudaEventRecord(sc.kernel_stop, sc.stream));
CUDA_CHECK(cudaGetLastError());
// Read back cluster count into pinned buffer
CUDA_CHECK(cudaMemcpyAsync(sc.h_cluster_count, sc.d_cluster_count,
sizeof(uint32_t), cudaMemcpyDeviceToHost,
sc.stream));
// Synchronize to ensure count is available before the CPU reads
// clusters
CUDA_CHECK(cudaStreamSynchronize(sc.stream));
record_kernel_time(sc);
// Clamp to max in case of overflow
uint32_t n_found = *sc.h_cluster_count;
n_found = std::min(n_found, static_cast<uint32_t>(m_capacity));
// Read back clusters
m_clusters.set_frame_number(frame_number);
if (n_found > 0) {
append_device_clusters_to(m_clusters, sc, n_found);
}
}
/**
* @brief Batched cluster finding across multiple frames, using n_streams
* CUDA streams to overlap H2D transfer, kernel, and D2H transfer.
*
* Returns one ClusterVector per input frame (with frame_number set to
* first_frame + i).
*/
std::vector<ClusterVector<ClusterType>>
find_clusters_batched(NDView<FRAME_TYPE, 3> frames,
uint64_t first_frame = 0) {
if (m_pedestal_dirty) {
sync_pedestal_to_device();
m_pedestal_dirty = false;
}
const size_t n_frames = frames.shape(0);
const uint32_t n_pd_samples =
static_cast<uint32_t>(m_pedestal.n_samples());
std::vector<ClusterVector<ClusterType>> results;
results.reserve(n_frames);
for (size_t i = 0; i < n_frames; ++i) {
results.emplace_back();
results.back().set_frame_number(first_frame + i);
}
const size_t n_rounds = (n_frames + n_streams - 1) / n_streams;
for (size_t round = 0; round < n_rounds; ++round) {
// Launch phase: fan out kernels on all streams for this round
for (int k = 0; k < n_streams; ++k) {
// OOB guard
const size_t frame_idx = round * n_streams + k;
if (frame_idx >= n_frames)
continue;
auto &sc_k = v_sc[k];
const FRAME_TYPE *h_src =
frames.data() + frame_idx * m_image_size;
std::memcpy(sc_k.h_frame, h_src, m_image_bytes);
CUDA_CHECK(cudaMemsetAsync(sc_k.d_cluster_count, 0,
sizeof(uint32_t), sc_k.stream));
CUDA_CHECK(
cudaMemcpyAsync(sc_k.d_frame, sc_k.h_frame, m_image_bytes,
cudaMemcpyHostToDevice, sc_k.stream));
CUDA_CHECK(cudaEventRecord(sc_k.kernel_start, sc_k.stream));
device::find_clusters_in_single_frame<ClusterType, FRAME_TYPE,
PEDESTAL_TYPE>
<<<grid, block, shmem_bytes, sc_k.stream>>>(
sc_k.d_frame, sc_k.d_pd_mean, sc_k.d_pd_sum,
sc_k.d_pd_sum2, n_pd_samples, m_nSigma, nrows, ncols,
sc_k.d_clusters, sc_k.d_cluster_count,
static_cast<uint32_t>(m_capacity));
CUDA_CHECK(cudaEventRecord(sc_k.kernel_stop, sc_k.stream));
CUDA_CHECK(cudaGetLastError());
// Queue count D2H immediately after the kernel
CUDA_CHECK(cudaMemcpyAsync(
sc_k.h_cluster_count, sc_k.d_cluster_count,
sizeof(uint32_t), cudaMemcpyDeviceToHost, sc_k.stream));
}
// Drain phase: fan in results from all streams
for (int k = 0; k < n_streams; ++k) {
const size_t frame_idx = round * n_streams + k;
if (frame_idx >= n_frames)
continue;
auto &sc_k = v_sc[k];
// Wait for memset -> H2D -> kernel -> count D2H
CUDA_CHECK(cudaStreamSynchronize(sc_k.stream));
record_kernel_time(sc_k);
uint32_t n_found = *sc_k.h_cluster_count;
n_found = std::min<uint32_t>(n_found,
static_cast<uint32_t>(m_capacity));
if (n_found > 0) {
append_device_clusters_to(results[frame_idx], sc_k,
n_found);
}
}
}
return results;
}
float avg_kernel_time_ms() const {
return m_frames_processed > 0 ? m_total_kernel_ms / m_frames_processed
: 0.0f;
}
void reset_timers() {
m_total_kernel_ms = 0.0f;
m_frames_processed = 0;
}
private:
/**
* Upload the current host pedestal (mean, sum, sum2) to every stream's
* device buffers. Called lazily before a find_clusters call when the
* host pedestal has been updated.
*/
void sync_pedestal_to_device() {
// These return-by-value NDArrays must stay alive until the async
// copies complete, so we synchronise at the end before they go out
// of scope.
NDArray<PEDESTAL_TYPE, 2> h_mean = m_pedestal.mean();
NDArray<PEDESTAL_TYPE, 2> h_sum = m_pedestal.get_sum();
NDArray<PEDESTAL_TYPE, 2> h_sum2 = m_pedestal.get_sum2();
const size_t bytes = m_image_size * sizeof(PEDESTAL_TYPE);
for (auto &sc : v_sc) {
CUDA_CHECK(cudaMemcpyAsync(sc.d_pd_mean, h_mean.data(), bytes,
cudaMemcpyHostToDevice, sc.stream));
CUDA_CHECK(cudaMemcpyAsync(sc.d_pd_sum, h_sum.data(), bytes,
cudaMemcpyHostToDevice, sc.stream));
CUDA_CHECK(cudaMemcpyAsync(sc.d_pd_sum2, h_sum2.data(), bytes,
cudaMemcpyHostToDevice, sc.stream));
}
for (auto &sc : v_sc)
CUDA_CHECK(cudaStreamSynchronize(sc.stream));
}
/**
* Copy n_found clusters from sc.d_clusters into the given ClusterVector
* and block on the transfer.
*/
void append_device_clusters_to(ClusterVector<ClusterType> &cv, SC &sc,
uint32_t n_found) {
CUDA_CHECK(cudaMemcpyAsync(sc.h_clusters, sc.d_clusters,
n_found * sizeof(ClusterType),
cudaMemcpyDeviceToHost, sc.stream));
CUDA_CHECK(cudaStreamSynchronize(sc.stream));
for (uint32_t i = 0; i < n_found; ++i)
cv.push_back(sc.h_clusters[i]);
}
void record_kernel_time(SC &sc) {
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, sc.kernel_start, sc.kernel_stop));
m_total_kernel_ms += ms;
m_frames_processed++;
}
};
} // namespace aare
python/aare/ClusterFinder.py
+30 -6
View File
@@ -55,25 +55,49 @@ def _cuda_available():
def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
capacity=1024, n_streams=1):
max_clusters_per_frame=2048, n_streams=4):
"""
Factory function to create a ClusterFinderCUDA object. Provides a cleaner
syntax for the templated ClusterFinderCUDA in C++. API mirrors
ClusterFinder() plus CUDA-specific knobs (n_streams).
ClusterFinder() plus CUDA-specific knobs.
Parameters
----------
image_size : tuple of (int, int)
Detector shape as (nrows, ncols).
cluster_size : tuple of (int, int), optional
Cluster window size; default (3, 3).
n_sigma : float, optional
Threshold in units of per-pixel pedestal standard deviation.
dtype : numpy dtype, optional
Cluster value type (np.int32 or np.float32).
max_clusters_per_frame : int, optional
Tight upper bound on clusters per frame. Determines the fixed-size D2H
transfer per frame. Set this high enough to never truncate real frames
but as tight as possible to minimize PCIe traffic. Default 2048.
n_streams : int, optional
Number of CUDA streams for H2D/kernel/D2H pipelining. Default 4.
Example
-------
.. code-block:: python
from aare import ClusterFinderCUDA
cf = ClusterFinderCUDA(image_size=(512, 1024),
cf = ClusterFinderCUDA(image_size=(400, 400),
cluster_size=(3, 3),
n_sigma=5,
n_streams=4)
n_streams=5)
for frame in pedestal_frames:
cf.push_pedestal_frame(frame)
# Batched (recommended for throughput)
results = cf.find_clusters_batched(frames_3d, first_frame=0)
# Or single-frame (one launch per frame)
for i, frame in enumerate(data_frames):
cf.find_clusters(frame, frame_number=i)
clusters = cf.steal_clusters()
clusters = cf.steal_clusters()
"""
if not _cuda_available():
raise RuntimeError(
@@ -84,7 +108,7 @@ def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
cls = _get_class("ClusterFinderCUDA", cluster_size, dtype)
return cls(image_size,
n_sigma=n_sigma,
capacity=capacity,
max_clusters_per_frame=max_clusters_per_frame,
n_streams=n_streams)
def ClusterCollector(clusterfindermt, dtype=np.int32):
python/src/bind_ClusterFinderCUDA.hpp
+42 -20
View File
@@ -8,8 +8,8 @@
#include <cstdint>
#include <pybind11/pybind11.h>
// #include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
#include <pybind11/stl.h>
// #include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
@@ -28,18 +28,20 @@ void define_ClusterFinderCUDA(py::module &m, const std::string &typestr) {
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
using CF = ClusterFinderCUDA<ClusterType, uint16_t, pd_type>;
using ContigArr =
py::array_t<uint16_t, py::array::c_style | py::array::forcecast>;
py::class_<CF>(m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t, int>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000,
py::arg("n_streams") = 1)
.def(py::init<Shape<2>, float, size_t, int>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0f,
py::arg("max_clusters_per_frame") = 2048, py::arg("n_streams") = 4)
.def_property(
"nSigma", &CF::get_nSigma, &CF::set_nSigma,
R"(Number of sigma above the pedestal to consider a photon during cluster finding.)")
.def("push_pedestal_frame",
[](CF &self, py::array_t<uint16_t> frame) {
[](CF &self, ContigArr frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
@@ -63,39 +65,59 @@ void define_ClusterFinderCUDA(py::module &m, const std::string &typestr) {
.def(
"steal_clusters",
[](CF &self, bool realloc_same_capacity) {
ClusterVector<ClusterType> clusters =
self.steal_clusters(realloc_same_capacity);
return clusters;
return std::move(self.steal_clusters(realloc_same_capacity));
},
py::arg("realloc_same_capacity") = false)
py::arg("realloc_same_capacity") = true)
.def(
"find_clusters",
[](CF &self, py::array_t<uint16_t> frame, uint64_t frame_number) {
[](CF &self, ContigArr frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
},
py::arg("frame"), py::arg("frame_number") = 0)
py::arg("frame"), py::arg("frame_number") = 0,
py::call_guard<py::gil_scoped_release>())
.def(
"find_clusters_batched",
[](CF &self, py::array_t<uint16_t> frames, uint64_t first_frame) {
// frames is expected as a 3D numpy array (n_frames, nrows,
// ncols)
[](CF &self, ContigArr frames, uint64_t first_frame) {
auto view = make_view_3d(frames);
return self.find_clusters_batched(view, first_frame);
},
py::arg("frames"), py::arg("first_frame") = 0,
R"(Process a 3D array of frames (n_frames, nrows, ncols) in parallel
across the configured CUDA streams. Returns a list of ClusterVector, one per
input frame.)")
py::call_guard<py::gil_scoped_release>(),
R"(Process a 3D array of frames (n_frames, nrows, ncols) using
n_streams CUDA streams for H2D/kernel/D2H pipelining. Returns a
list of ClusterVector, one per input frame. The input array is
converted to C-contiguous uint16 if needed.)")
.def("avg_kernel_time_ms", &CF::avg_kernel_time_ms,
R"(Average kernel execution time per frame in milliseconds,
excluding PCIe transfers. Use wall_time - avg_kernel_time to estimate transfer overhead.)")
excluding PCIe transfers.)")
.def("reset_timers", &CF::reset_timers,
R"(Reset the internal kernel timing counters.)");
R"(Reset the internal kernel timing counters.)")
.def(
"register_input_buffer",
[](CF &self, py::array arr) {
auto info = arr.request();
self.register_input_buffer(
info.ptr, static_cast<size_t>(info.size) *
static_cast<size_t>(info.itemsize));
},
R"(Pin a numpy array as a locked host buffer so that
find_clusters_batched transfers it at full DMA bandwidth
(~22 GB/s) instead of going through the CUDA driver's
internal staging (~15 GB/s for pageable memory).
Call once before the processing loop with the full data
array. Slices of that array passed to find_clusters_batched
lie within the registered region and benefit automatically.
Call unregister_input_buffer() when done.)")
.def("unregister_input_buffer", &CF::unregister_input_buffer,
"Release the previously pinned input buffer.");
}
} // namespace aare
python/tests/ClusterFinderCUDA.ipynb
+52 -30
View File
File diff suppressed because one or more lines are too long
src/ClusterFinderCUDA.test.cu
+60 -36
View File
@@ -257,18 +257,20 @@ void feed_pedestal(
int main(int argc, char *argv[]) {
// Parse arguments
// const char *default_pedestal =
// "/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/Calibration_CZT/"
// "2025Sept_m694/Sn25300eV/500_us_voltage_40kV/"
// "250922_CZTonly_Pedestal_Tp15C_tint_500_master_0.json";
const char *default_pedestal =
"/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/Calibration_CZT/"
"2025Sept_m694/Sn25300eV/500_us_voltage_40kV/"
"250922_CZTonly_Pedestal_Tp15C_tint_500_master_0.json";
// const char* default_pedestal =
// "/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/November2025/sparse/dynamic/si_pedestal_200keV_DYNAMIC_tint_20us_master_0.json";
"/mnt/sls_det_storage/matterhorn_data/aare_test_data/Moench03new/"
"cu_half_speed_master_4.json";
// const char *default_data =
// "/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/Calibration_CZT/"
// "2025Sept_m694/Sn25300eV/500_us_voltage_40kV/"
// "250922_CZTonly_Xray_Tp15C_tint_500_master_0.json";
const char *default_data =
"/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/Calibration_CZT/"
"2025Sept_m694/Sn25300eV/500_us_voltage_40kV/"
"250922_CZTonly_Xray_Tp15C_tint_500_master_0.json";
// const char* default_data =
// "/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/November2025/sparse/dynamic/si_data_200keV_DYNAMIC_tint_20us_master_0.json";
"/mnt/sls_det_storage/matterhorn_data/aare_test_data/Moench03new/"
"cu_half_speed_master_4.json";
std::filesystem::path pedestal_path(argc > 1 ? argv[1] : default_pedestal);
std::filesystem::path data_path(argc > 2 ? argv[2] : default_data);
@@ -283,7 +285,7 @@ int main(int argc, char *argv[]) {
}
// Defaults: Adjust depending on the test dataset used
size_t n_pedestal_frames = 6000;
size_t n_pedestal_frames = 1000;
size_t n_data_frames = 10000;
double nSigma = 5.0;
@@ -381,6 +383,7 @@ int main(int argc, char *argv[]) {
// Pre-read all frames into memory to remove I/O from timing
std::vector<aare::NDArray<FRAME_TYPE, 2>> frames;
frames.reserve(use_data);
data_file.seek(n_pedestal_frames);
for (size_t i = 0; i < use_data; ++i) {
auto f = data_file.read_frame();
// Copy into NDArray for consistent access
@@ -440,6 +443,7 @@ int main(int argc, char *argv[]) {
size_t gpu_total_clusters = 0;
// --- Single frame on a single CUDA stream ---
/*
{
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> cuda_cf(
{ROWS, COLS}, nSigma);
@@ -465,56 +469,68 @@ int main(int argc, char *argv[]) {
printf("GPU: %zu clusters in %.1f ms (%.2f ms/frame)\n",
gpu_total_clusters, gpu_time, gpu_time / use_data);
}
*/
// --- Batched H2D + multi-stream (enable this block and disable the one
// above to benchmark the batched path against the CPU results) ---
/*
{
constexpr size_t BATCH_SIZE = 100;
constexpr int N_STREAMS = 2;
constexpr size_t BATCH_SIZE = 2000;
constexpr int N_STREAMS = 5;
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> cuda_cf(
{ROWS, COLS}, nSigma, 1'000'000, N_STREAMS);
{ROWS, COLS}, nSigma, 4096, N_STREAMS);
feed_pedestal(cuda_cf, pedestal_frames);
// Contiguous staging buffer reused across batches
// Contiguous staging buffer reused across batches — registered as
// pinned so that H2D transfers run at DMA bandwidth (~22 GB/s) instead
// of going through the CUDA driver's internal staging (~15 GB/s for
// pageable memory).
std::vector<FRAME_TYPE> batch_buffer(BATCH_SIZE * ROWS * COLS);
cuda_cf.register_input_buffer(batch_buffer.data(),
batch_buffer.size() * sizeof(FRAME_TYPE));
const size_t n_batches = (use_data + BATCH_SIZE - 1) / BATCH_SIZE;
double pack_ms = 0.0, gpu_ms = 0.0;
Timer t;
t.start();
for (size_t bi = 0; bi < n_batches; ++bi) {
const size_t offset = bi * BATCH_SIZE;
const size_t offset = bi * BATCH_SIZE;
const size_t actual_batch = std::min(BATCH_SIZE, use_data - offset);
t.start();
pack_frame_batch(frames, offset, actual_batch, batch_buffer);
pack_ms += t.elapsed_ms();
aare::NDView<FRAME_TYPE, 3> batch_view(
batch_buffer.data(),
{static_cast<ssize_t>(actual_batch), ROWS, COLS});
auto batch_results = cuda_cf.find_clusters_batched(batch_view,
offset);
t.start();
auto batch_results =
cuda_cf.find_clusters_batched(batch_view, offset);
gpu_ms += t.elapsed_ms();
for (size_t f = 0; f < actual_batch; ++f) {
auto& cv = batch_results[f];
auto& out = gpu_results[offset + f];
auto &cv = batch_results[f];
auto &out = gpu_results[offset + f];
out.clear();
out.reserve(cv.size());
for (size_t j = 0; j < cv.size(); ++j) out.push_back(cv[j]);
for (size_t j = 0; j < cv.size(); ++j)
out.push_back(cv[j]);
gpu_total_clusters += out.size();
}
}
double gpu_time = t.elapsed_ms();
printf("GPU(batched): %zu clusters in %.1f ms (%.2f ms/frame, batch=%zu,
streams=%d)\n", gpu_total_clusters, gpu_time, gpu_time / use_data,
BATCH_SIZE, N_STREAMS);
cuda_cf.unregister_input_buffer();
printf("GPU(batched): %zu clusters — pack=%.1f ms GPU=%.1f ms "
"total=%.1f ms"
" (%.2f µs/frame, batch=%zu, streams=%d)\n",
gpu_total_clusters, pack_ms, gpu_ms, pack_ms + gpu_ms,
1000.0 * (pack_ms + gpu_ms) / use_data, BATCH_SIZE, N_STREAMS);
}
*/
// =========================================================================
// Phase 5: Comparison
@@ -579,10 +595,10 @@ int main(int argc, char *argv[]) {
// =========================================================================
// Phase 6: Per-frame timing benchmark
// =========================================================================
printf("\n--- Phase 6: Detailed timing (%d iterations) ---\n", 1000);
printf("\n--- Phase 6: Detailed timing (%d iterations) ---\n", 10000);
if (use_data > 0) {
const int N_ITER = 1000;
const int N_ITER = 10000;
const auto &bench_frame = frames[0];
// CPU benchmark
@@ -614,7 +630,7 @@ int main(int argc, char *argv[]) {
// Warmup
cuda_cf.find_clusters(bench_frame.view(), 0);
cuda_cf.steal_clusters(true);
// cuda_cf.steal_clusters(true);
Timer t;
t.start();
@@ -629,18 +645,21 @@ int main(int argc, char *argv[]) {
// --- GPU benchmark (batched + multi-streamed) ---
{
constexpr size_t BATCH_SIZE = 500;
constexpr int N_STREAMS = 10;
constexpr size_t BATCH_SIZE = 2000;
constexpr int N_STREAMS = 5;
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>
cuda_cf({ROWS, COLS}, nSigma, 1'000'000, N_STREAMS);
cuda_cf({ROWS, COLS}, nSigma, 4096, N_STREAMS);
feed_pedestal(cuda_cf, pedestal_frames);
// Build one contiguous batch of BATCH_SIZE copies of bench_frame
// Build one contiguous batch of BATCH_SIZE copies of bench_frame.
// Register as pinned for DMA-speed H2D transfers.
std::vector<FRAME_TYPE> batch(BATCH_SIZE * ROWS * COLS);
for (size_t k = 0; k < BATCH_SIZE; ++k)
std::memcpy(batch.data() + k * ROWS * COLS, bench_frame.data(),
ROWS * COLS * sizeof(FRAME_TYPE));
cuda_cf.register_input_buffer(batch.data(),
batch.size() * sizeof(FRAME_TYPE));
aare::NDView<FRAME_TYPE, 3> batch_view(
batch.data(), {static_cast<ssize_t>(BATCH_SIZE), ROWS, COLS});
@@ -657,11 +676,16 @@ int main(int argc, char *argv[]) {
const size_t n_iter_batches =
(N_ITER + BATCH_SIZE - 1) / BATCH_SIZE;
std::vector<aare::ClusterVector<ClusterType>> batch_results;
Timer t;
t.start();
for (size_t b = 0; b < n_iter_batches; ++b) {
(void)cuda_cf.find_clusters_batched(batch_view, b * BATCH_SIZE);
batch_results =
cuda_cf.find_clusters_batched(batch_view, b * BATCH_SIZE);
}
cuda_cf.unregister_input_buffer();
printf("GPU(batched): %.3f ms/frame (H2D + kernel + D2H, "
"batch=%zu, streams=%d)\n",
t.elapsed_ms() / (n_iter_batches * BATCH_SIZE), BATCH_SIZE,
src/ClusterFinderCUDA_old.test.cu
+683
View File
@@ -0,0 +1,683 @@
// SPDX-License-Identifier: MPL-2.0
#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <filesystem>
#include <numeric>
#include <vector>
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderCUDA_old.hpp"
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/NDArray.hpp"
#include "aare/Pedestal.hpp"
#include "aare/defs.hpp"
#include "aare/utils/batch.hpp"
// _____________
//
// Timing helper
// _____________
struct Timer {
using clock = std::chrono::high_resolution_clock;
clock::time_point t0;
void start() { t0 = clock::now(); }
double elapsed_ms() const {
return std::chrono::duration<double, std::milli>(clock::now() - t0)
.count();
}
};
// __________________
//
// Cluster comparison
// __________________
template <typename ClusterType> struct ClusterComparison {
size_t cpu_count = 0;
size_t gpu_count = 0;
size_t matched = 0;
size_t position_mismatch = 0;
size_t data_mismatch = 0;
size_t cpu_only = 0;
size_t gpu_only = 0;
};
// Sort clusters by (y, x) for deterministic comparison
template <typename ClusterType>
void sort_clusters(std::vector<ClusterType> &clusters) {
std::sort(clusters.begin(), clusters.end(),
[](const ClusterType &a, const ClusterType &b) {
if (a.y != b.y)
return a.y < b.y;
return a.x < b.x;
});
}
// Compare two sorted cluster lists
template <typename ClusterType>
ClusterComparison<ClusterType>
compare_clusters(std::vector<ClusterType> &cpu_clusters,
std::vector<ClusterType> &gpu_clusters) {
sort_clusters(cpu_clusters);
sort_clusters(gpu_clusters);
ClusterComparison<ClusterType> result;
result.cpu_count = cpu_clusters.size();
result.gpu_count = gpu_clusters.size();
size_t ci = 0, gi = 0;
while (ci < cpu_clusters.size() && gi < gpu_clusters.size()) {
const auto &cc = cpu_clusters[ci];
const auto &gc = gpu_clusters[gi];
if (cc.y == gc.y && cc.x == gc.x) {
// Same position/check data
bool data_ok = true;
constexpr int N =
ClusterType::cluster_size_x * ClusterType::cluster_size_y;
for (int k = 0; k < N; ++k) {
// if (cc.data[k] != gc.data[k]) { // a bit too strict
// espacially that pedestal update is slightly different
if (std::abs(cc.data[k] - gc.data[k]) > 5) {
data_ok = false;
break;
}
}
if (data_ok)
result.matched++;
else
result.data_mismatch++;
ci++;
gi++;
} else if (cc.y < gc.y || (cc.y == gc.y && cc.x < gc.x)) {
result.cpu_only++;
ci++;
} else {
result.gpu_only++;
gi++;
}
}
result.cpu_only += (cpu_clusters.size() - ci);
result.gpu_only += (gpu_clusters.size() - gi);
return result;
}
// ________________
//
// Cluster printing
// ________________
template <typename ClusterType>
void print_cluster_comparison(
const std::vector<std::vector<ClusterType>> &cpu_results,
const std::vector<std::vector<ClusterType>> &gpu_results,
size_t max_per_frame = 10, size_t max_frames = 100) {
constexpr int NX = ClusterType::cluster_size_x;
constexpr int NY = ClusterType::cluster_size_y;
constexpr int N = NX * NY;
size_t frames_shown = 0;
for (size_t fi = 0; fi < cpu_results.size() && frames_shown < max_frames;
++fi) {
if (cpu_results[fi].empty() && gpu_results[fi].empty())
continue;
size_t n_cpu = cpu_results[fi].size();
size_t n_gpu = gpu_results[fi].size();
printf("\n Frame %zu: CPU=%zu clusters, GPU=%zu clusters\n", fi, n_cpu,
n_gpu);
// Merge-walk over sorted lists (assumes already sorted by y,x)
size_t ci = 0, gi = 0;
size_t shown = 0;
while ((ci < n_cpu || gi < n_gpu) && shown < max_per_frame) {
bool have_cpu = ci < n_cpu;
bool have_gpu = gi < n_gpu;
// Determine if current entries match in position
bool same_pos = have_cpu && have_gpu &&
cpu_results[fi][ci].x == gpu_results[fi][gi].x &&
cpu_results[fi][ci].y == gpu_results[fi][gi].y;
if (same_pos) {
const auto &cc = cpu_results[fi][ci];
const auto &gc = gpu_results[fi][gi];
// Check if data differs
bool differs = false;
for (int k = 0; k < N; ++k) {
if (cc.data[k] != gc.data[k]) {
differs = true;
break;
}
}
printf(" CPU and GPU clusters found at SAME position "
"(col=%3d, row=%3d)\n %s\n",
cc.x, cc.y, differs ? "DATA MISMATCH" : "MATCH");
printf(" CPU: [");
for (int k = 0; k < N; ++k) {
if (k)
printf(", ");
printf("%6d", static_cast<int>(cc.data[k]));
}
printf("]\n");
if (differs) {
printf(" GPU: [");
for (int k = 0; k < N; ++k) {
if (k)
printf(", ");
printf("%6d", static_cast<int>(gc.data[k]));
}
printf("]\n");
printf(" diff: [");
for (int k = 0; k < N; ++k) {
if (k)
printf(", ");
int d = static_cast<int>(gc.data[k]) -
static_cast<int>(cc.data[k]);
printf("%+6d", d);
}
printf("]\n");
}
ci++;
gi++;
shown++;
} else if (!have_gpu ||
(have_cpu &&
(cpu_results[fi][ci].y < gpu_results[fi][gi].y ||
(cpu_results[fi][ci].y == gpu_results[fi][gi].y &&
cpu_results[fi][ci].x < gpu_results[fi][gi].x)))) {
const auto &cc = cpu_results[fi][ci];
printf(" (%3d, %3d) CPU ONLY\n", cc.x, cc.y);
printf(" CPU: [");
for (int k = 0; k < N; ++k) {
if (k)
printf(", ");
printf("%6d", static_cast<int>(cc.data[k]));
}
printf("]\n");
ci++;
shown++;
} else {
const auto &gc = gpu_results[fi][gi];
printf(" (%3d, %3d) GPU ONLY\n", gc.x, gc.y);
printf(" GPU: [");
for (int k = 0; k < N; ++k) {
if (k)
printf(", ");
printf("%6d", static_cast<int>(gc.data[k]));
}
printf("]\n");
gi++;
shown++;
}
}
frames_shown++;
}
}
// _________________________________________
//
// Helpers for the updated (CPU-parity) API
// _________________________________________
// Copy a ClusterVector into a std::vector<ClusterType> for downstream
// comparison code that expects the latter.
template <typename Finder, typename ClusterType>
void drain_into(Finder &f, std::vector<ClusterType> &out) {
auto cv = f.steal_clusters(true);
out.clear();
out.reserve(cv.size());
for (size_t j = 0; j < cv.size(); ++j)
out.push_back(cv[j]);
}
// Feed a set of cached pedestal frames through any finder exposing the
// CPU-style push_pedestal_frame(NDView) method. Works for both
// ClusterFinder and ClusterFinderCUDA.
template <typename Finder, typename FRAME_TYPE>
void feed_pedestal(
Finder &f, const std::vector<aare::NDArray<FRAME_TYPE, 2>> &ped_frames) {
for (const auto &frame : ped_frames) {
f.push_pedestal_frame(frame.view());
}
}
// ____________
//
// Main
// ____________
int main(int argc, char *argv[]) {
// Parse arguments
// const char *default_pedestal =
// "/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/Calibration_CZT/"
// "2025Sept_m694/Sn25300eV/500_us_voltage_40kV/"
// "250922_CZTonly_Pedestal_Tp15C_tint_500_master_0.json";
const char *default_pedestal =
"/mnt/sls_det_storage/matterhorn_data/aare_test_data/Moench03new/"
"cu_half_speed_master_4.json";
// const char *default_data =
// "/mnt/sls_det_storage/highZ_data/CZT_Vienna/Khalil/Calibration_CZT/"
// "2025Sept_m694/Sn25300eV/500_us_voltage_40kV/"
// "250922_CZTonly_Xray_Tp15C_tint_500_master_0.json";
const char *default_data =
"/mnt/sls_det_storage/matterhorn_data/aare_test_data/Moench03new/"
"cu_half_speed_master_4.json";
std::filesystem::path pedestal_path(argc > 1 ? argv[1] : default_pedestal);
std::filesystem::path data_path(argc > 2 ? argv[2] : default_data);
if (!std::filesystem::exists(pedestal_path)) {
fprintf(stderr, "Pedestal file not found: %s\n", pedestal_path.c_str());
return 1;
}
if (!std::filesystem::exists(data_path)) {
fprintf(stderr, "Data file not found: %s\n", data_path.c_str());
return 1;
}
// Defaults: Adjust depending on the test dataset used
size_t n_pedestal_frames = 1000;
size_t n_data_frames = 10000;
double nSigma = 5.0;
if (argc > 3)
n_pedestal_frames = std::atol(argv[3]);
if (argc > 4)
n_data_frames = std::atol(argv[4]);
if (argc > 5)
nSigma = std::atof(argv[5]);
// Detector geometry from master file
constexpr uint8_t cs_x = 3;
constexpr uint8_t cs_y = 3;
using ClusterType = aare::Cluster<int32_t, cs_x, cs_y>;
using FRAME_TYPE = uint16_t;
using PEDESTAL_TYPE = double;
// Read actual frame dimensions from the pedestal file
ssize_t ROWS, COLS;
{
aare::File probe(pedestal_path, "r");
auto first_frame = probe.read_frame();
ROWS = first_frame.rows();
COLS = first_frame.cols();
}
printf("=== Cluster Finder: CPU vs CUDA (OLD) ===\n");
printf("Detector: %zu x %zu\n", ROWS, COLS);
printf("Cluster: %d x %d\n", ClusterType::cluster_size_x,
ClusterType::cluster_size_y);
printf("nSigma: %.1f\n", nSigma);
// =========================================================================
// Phase 1: Build pedestal from dark frames (sanity check only + frame
// cache)
// =========================================================================
//
// Neither ClusterFinder nor ClusterFinderCUDA needs an external Pedestal
// object; both build their own via push_pedestal_frame. We still read the
// pedestal file once, but cache the frames in memory so every subsequent
// finder can be fed without re-hitting disk. We also build a standalone
// Pedestal purely to print a sanity check.
// =========================================================================
printf("\n--- Phase 1: Pedestal accumulation (sanity check + cache) ---\n");
std::vector<aare::NDArray<FRAME_TYPE, 2>> pedestal_frames;
aare::Pedestal<PEDESTAL_TYPE> pedestal(ROWS, COLS, 1000);
{
aare::File ped_file(pedestal_path, "r");
size_t total_ped = ped_file.total_frames();
size_t use_ped =
(n_pedestal_frames == 0 || n_pedestal_frames > total_ped)
? total_ped
: n_pedestal_frames;
printf("Pedestal frames: %zu / %zu\n", use_ped, total_ped);
pedestal_frames.reserve(use_ped);
Timer t;
t.start();
for (size_t i = 0; i < use_ped; ++i) {
auto frame = ped_file.read_frame();
auto view = frame.view<FRAME_TYPE>();
// Copy into a standalone NDArray that we can reuse as many times
// as we have finders to feed.
aare::NDArray<FRAME_TYPE, 2> arr({ROWS, COLS});
for (ssize_t r = 0; r < ROWS; ++r)
for (ssize_t c = 0; c < COLS; ++c)
arr(r, c) = view(r, c);
pedestal_frames.push_back(std::move(arr));
pedestal.push_no_update(view);
}
pedestal.update_mean();
printf("Pedestal read+cached+built in %.1f ms\n", t.elapsed_ms());
}
printf("Pedestal mean[0,0] = %.2f, std[0,0] = %.4f\n", pedestal.mean(0, 0),
pedestal.std(0, 0));
// =========================================================================
// Phase 2: Read data frames
// =========================================================================
printf("\n--- Phase 2: Read data frames ---\n");
aare::File data_file(data_path, "r");
size_t total_data = data_file.total_frames();
size_t use_data = std::min(n_data_frames, total_data);
printf("Data frames: %zu / %zu\n", use_data, total_data);
// Pre-read all frames into memory to remove I/O from timing
std::vector<aare::NDArray<FRAME_TYPE, 2>> frames;
frames.reserve(use_data);
data_file.seek(n_pedestal_frames);
for (size_t i = 0; i < use_data; ++i) {
auto f = data_file.read_frame();
// Copy into NDArray for consistent access
aare::NDArray<FRAME_TYPE, 2> arr({ROWS, COLS});
auto view = f.view<FRAME_TYPE>();
for (size_t r = 0; r < ROWS; ++r)
for (size_t c = 0; c < COLS; ++c)
arr(r, c) = view(r, c);
frames.push_back(std::move(arr));
}
printf("Frames loaded into memory\n");
// =========================================================================
// Phase 3: Sequential (CPU) cluster finding
// =========================================================================
printf("\n--- Phase 3: CPU ClusterFinder ---\n");
std::vector<std::vector<ClusterType>> cpu_results(use_data);
size_t cpu_total_clusters = 0;
{
// Build a ClusterFinder with the same pedestal
aare::ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> cf(
{ROWS, COLS}, nSigma);
feed_pedestal(cf, pedestal_frames);
Timer t;
t.start();
for (size_t i = 0; i < use_data; ++i) {
cf.find_clusters(frames[i].view(), static_cast<uint64_t>(i));
drain_into(cf, cpu_results[i]);
cpu_total_clusters += cpu_results[i].size();
}
double cpu_time = t.elapsed_ms();
printf("CPU: %zu clusters in %.1f ms (%.2f ms/frame)\n",
cpu_total_clusters, cpu_time, cpu_time / use_data);
}
// =========================================================================
// Phase 4: CUDA cluster finding
// =========================================================================
//
// The API mirrors ClusterFinder: push_pedestal_frame to train, then
// find_clusters / steal_clusters for each frame. H2D transfer and kernel
// launch happen internally.
//
// Toggle between the single-stream and batched paths by swapping which
// block is enabled. They both write into gpu_results so only one at a
// time makes sense.
// =========================================================================
printf("\n--- Phase 4: CUDA ClusterFinder ---\n");
std::vector<std::vector<ClusterType>> gpu_results(use_data);
size_t gpu_total_clusters = 0;
// --- Single frame on a single CUDA stream ---
/*
{
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> cuda_cf(
{ROWS, COLS}, nSigma);
feed_pedestal(cuda_cf, pedestal_frames);
// Warmup: first CUDA call pays driver/context init overhead. The
// pedestal drifts slightly during this single frame, which is
// acceptable for timing purposes.
cuda_cf.find_clusters(frames[0].view(), 0);
cuda_cf.steal_clusters(true);
Timer t;
t.start();
for (size_t i = 0; i < use_data; ++i) {
cuda_cf.find_clusters(frames[i].view(), static_cast<uint64_t>(i));
drain_into(cuda_cf, gpu_results[i]);
gpu_total_clusters += gpu_results[i].size();
}
double gpu_time = t.elapsed_ms();
printf("GPU: %zu clusters in %.1f ms (%.2f ms/frame)\n",
gpu_total_clusters, gpu_time, gpu_time / use_data);
}
*/
// --- Batched H2D + multi-stream (enable this block and disable the one
// above to benchmark the batched path against the CPU results) ---
{
constexpr size_t BATCH_SIZE = 2000;
constexpr int N_STREAMS = 5;
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> cuda_cf(
{ROWS, COLS}, nSigma, 50'000, N_STREAMS);
feed_pedestal(cuda_cf, pedestal_frames);
// Contiguous staging buffer reused across batches
std::vector<FRAME_TYPE> batch_buffer(BATCH_SIZE * ROWS * COLS);
const size_t n_batches = (use_data + BATCH_SIZE - 1) / BATCH_SIZE;
double pack_ms = 0.0, gpu_ms = 0.0;
Timer t;
for (size_t bi = 0; bi < n_batches; ++bi) {
const size_t offset = bi * BATCH_SIZE;
const size_t actual_batch = std::min(BATCH_SIZE, use_data - offset);
t.start();
pack_frame_batch(frames, offset, actual_batch, batch_buffer);
pack_ms += t.elapsed_ms();
aare::NDView<FRAME_TYPE, 3> batch_view(
batch_buffer.data(),
{static_cast<ssize_t>(actual_batch), ROWS, COLS});
t.start();
auto batch_results =
cuda_cf.find_clusters_batched(batch_view, offset);
gpu_ms += t.elapsed_ms();
for (size_t f = 0; f < actual_batch; ++f) {
auto &cv = batch_results[f];
auto &out = gpu_results[offset + f];
out.clear();
out.reserve(cv.size());
for (size_t j = 0; j < cv.size(); ++j)
out.push_back(cv[j]);
gpu_total_clusters += out.size();
}
}
printf("GPU(batched): %zu clusters — pack=%.1f ms GPU=%.1f ms "
"total=%.1f ms"
" (%.2f µs/frame, batch=%zu, streams=%d)\n",
gpu_total_clusters, pack_ms, gpu_ms, pack_ms + gpu_ms,
1000.0 * (pack_ms + gpu_ms) / use_data, BATCH_SIZE, N_STREAMS);
}
// =========================================================================
// Phase 5: Comparison
// =========================================================================
printf("\n--- Phase 5: Comparison ---\n");
size_t total_matched = 0;
size_t total_data_mismatch = 0;
size_t total_cpu_only = 0;
size_t total_gpu_only = 0;
size_t frames_with_differences = 0;
for (size_t i = 0; i < use_data; ++i) {
auto result = compare_clusters(cpu_results[i], gpu_results[i]);
total_matched += result.matched;
total_data_mismatch += result.data_mismatch;
total_cpu_only += result.cpu_only;
total_gpu_only += result.gpu_only;
bool has_diff = (result.cpu_only > 0 || result.gpu_only > 0 ||
result.data_mismatch > 0);
if (has_diff) {
frames_with_differences++;
// Print details for first few mismatching frames
if (frames_with_differences <= 5) {
printf(" Frame %zu: CPU=%zu GPU=%zu matched=%zu "
"data_mismatch=%zu cpu_only=%zu gpu_only=%zu\n",
i, result.cpu_count, result.gpu_count, result.matched,
result.data_mismatch, result.cpu_only, result.gpu_only);
}
}
}
printf("\nSummary over %zu frames:\n", use_data);
printf(" CPU total clusters: %zu\n", cpu_total_clusters);
printf(" GPU total clusters: %zu\n", gpu_total_clusters);
printf(" Matched: %zu\n", total_matched);
printf(" Data mismatch: %zu\n", total_data_mismatch);
printf(" CPU only: %zu\n", total_cpu_only);
printf(" GPU only: %zu\n", total_gpu_only);
printf(" Frames with diffs: %zu / %zu\n", frames_with_differences,
use_data);
if (total_cpu_only == 0 && total_gpu_only == 0 &&
total_data_mismatch == 0) {
printf("\n*** PASS: CPU and GPU results match exactly ***\n");
} else {
printf("\n*** DIFFERENCES DETECTED ***\n");
}
// // Print detailed cluster comparison (side-by-side with diffs)
// if (cpu_total_clusters > 0 || gpu_total_clusters > 0) {
// size_t max_clusters_per_frame = 10;
// size_t max_frames = 100;
// printf("\n--- Cluster details (up to %zu frames, %zu clusters each)
// ---\n", max_frames, max_clusters_per_frame);
// print_cluster_comparison(cpu_results, gpu_results,
// max_clusters_per_frame, max_frames);
// }
// =========================================================================
// Phase 6: Per-frame timing benchmark
// =========================================================================
printf("\n--- Phase 6: Detailed timing (%d iterations) ---\n", 10000);
if (use_data > 0) {
const int N_ITER = 10000;
const auto &bench_frame = frames[0];
// CPU benchmark
{
aare::ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> cf(
{ROWS, COLS}, nSigma);
// Load pedestal
feed_pedestal(cf, pedestal_frames);
// Warmup
cf.find_clusters(bench_frame.view(), 0);
cf.steal_clusters(true);
Timer t;
t.start();
for (int iter = 0; iter < N_ITER; ++iter) {
cf.find_clusters(bench_frame.view(), 0);
cf.steal_clusters(true);
}
double cpu_per_frame = t.elapsed_ms() / N_ITER;
printf("CPU: %.3f ms/frame\n", cpu_per_frame);
}
// --- GPU benchmark (single frame, single stream) ---
{
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>
cuda_cf({ROWS, COLS}, nSigma);
feed_pedestal(cuda_cf, pedestal_frames);
// Warmup
cuda_cf.find_clusters(bench_frame.view(), 0);
cuda_cf.steal_clusters(true);
Timer t;
t.start();
for (int iter = 0; iter < N_ITER; ++iter) {
cuda_cf.find_clusters(bench_frame.view(), 0);
cuda_cf.steal_clusters(true);
}
printf("GPU: %.3f ms/frame (H2D + kernel + D2H, single "
"frame, single stream)\n",
t.elapsed_ms() / N_ITER);
}
// --- GPU benchmark (batched + multi-streamed) ---
{
constexpr size_t BATCH_SIZE = 2000;
constexpr int N_STREAMS = 5;
aare::ClusterFinderCUDA<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>
cuda_cf({ROWS, COLS}, nSigma, 50'000, N_STREAMS);
feed_pedestal(cuda_cf, pedestal_frames);
// Build one contiguous batch of BATCH_SIZE copies of bench_frame
std::vector<FRAME_TYPE> batch(BATCH_SIZE * ROWS * COLS);
for (size_t k = 0; k < BATCH_SIZE; ++k)
std::memcpy(batch.data() + k * ROWS * COLS, bench_frame.data(),
ROWS * COLS * sizeof(FRAME_TYPE));
aare::NDView<FRAME_TYPE, 3> batch_view(
batch.data(), {static_cast<ssize_t>(BATCH_SIZE), ROWS, COLS});
// Warmup. The kernel mutates the device-side pedestal for every
// non-photon pixel, so after a 500-frame warmup the pedestal
// state has drifted. Reset to a clean baseline by clearing the
// host pedestal and re-feeding the cached frames; this also
// re-arms the dirty flag so the next find_clusters re-uploads.
(void)cuda_cf.find_clusters_batched(batch_view, 0);
cuda_cf.clear_pedestal();
feed_pedestal(cuda_cf, pedestal_frames);
const size_t n_iter_batches =
(N_ITER + BATCH_SIZE - 1) / BATCH_SIZE;
Timer t;
t.start();
for (size_t b = 0; b < n_iter_batches; ++b) {
(void)cuda_cf.find_clusters_batched(batch_view, b * BATCH_SIZE);
}
printf("GPU(batched): %.3f ms/frame (H2D + kernel + D2H, "
"batch=%zu, streams=%d)\n",
t.elapsed_ms() / (n_iter_batches * BATCH_SIZE), BATCH_SIZE,
N_STREAMS);
}
}
printf("\nDone.\n");
return 0;
}
src/Makefile
+12 -5
View File
@@ -7,17 +7,24 @@ LDFLAGS := -L../build -L../build/_deps/fmt-build
LIBS := -laare_core -lfmt -lstdc++fs
DEFINES := -DAARE_LOG_LEVEL=logERROR
TARGET := test_cf_cuda
SRC := ClusterFinderCUDA.test.cu
DEP := $(SRC:.cu=.d)
TARGET_OLD := test_cf_cuda_old
TARGET := test_cf_cuda
all: $(TARGET)
SRC_OLD := ClusterFinderCUDA_old.test.cu
SRC := ClusterFinderCUDA.test.cu
DEP := $(SRC:.cu=.d) $(SRC_OLD:.cu=.d)
all: $(TARGET) $(TARGET_OLD)
$(TARGET): $(SRC) ../include/aare/clusterfinder_kernel.cuh
$(NVCC) -Xptxas=-v $(ARCH) $(CXXFLAGS) $(DEFINES) $(INCLUDES) $(LDFLAGS) $< -o $@ $(LIBS)
$(TARGET_OLD): $(SRC_OLD) ../include/aare/clusterfinder_kernel.cuh
$(NVCC) -Xptxas=-v $(ARCH) $(CXXFLAGS) $(DEFINES) $(INCLUDES) $(LDFLAGS) $< -o $@ $(LIBS)
clean:
rm -f $(TARGET) $(DEP)
rm -f $(TARGET) $(TARGET_OLD) $(DEP)
-include $(DEP)