Building PixelHistogramCUDA

Start here

Steps 1, 3 and 4 are done and verified. Steps 2 and 5 are written and compile, but are only verified by step 6. Details in Steps.

How the pipeline works

Two streams, two device slots, one histogram that stays on the GPU for the whole run.

• host framescaller-owned, pinned once with register_input_buffer
• copy streamH2D of batch i into slot i mod 2
• compute streamone fill kernel per batch, in order
• d_histogramuint32_t counts, bin-major
• values()transpose + cap, one D2H

Who waits for what

Expected case: copy-bound

If the kernel were slower: compute-bound

Throughput is one batch per max(copy time, kernel time). A third slot would only absorb jitter, not make the slower stage faster.

Kernel design

Fill kernel: bin-major, one thread per pixel

• 1D grid, one thread per pixel. Histogramming never looks at neighbouring pixels, so a pixel is just a flat index p. Consecutive threads take consecutive pixels, so each warp reads one contiguous chunk per frame.
• Each thread loops over the batch's frames and increments its own pixel's bins.
• No atomics. Only one thread ever writes pixel p's bins, and kernels never overlap because they share one in-order stream.

Bin-major is never worse, and much better when neighbouring pixels land in the same bin. What matters is how many distinct bins a warp hits, not how far apart they are.

No shared memory in the fill kernel

• It wouldn't be shared: each thread owns its pixel's bins, so it would just be a big private array per thread.
• 2 KB per pixel (1000 bins × 2 B) fits only ~50 pixels per SM, which collapses occupancy.
• Writing all bins back after each batch would move the whole histogram (~1.3 GB per batch).
• The 4090's 72 MB L2 already holds the hot bins: ~2 µs of L2 traffic per frame, against ~13 µs to copy one frame over PCIe.

Finalize kernel: transpose and cap Done

• bin_to_pixel_major turns the [nbins][npix] accumulator into [npix][nbins] StorageType, capping each count. One thread per output element: writes are coalesced, reads are strided but mostly served from cache.
• It writes into a temporary buffer and only reads the accumulator, so values() works mid-run and filling can continue afterwards.
• Cost: roughly 2–4 ms, against ~13 ms for the 320 MB D2H that follows. A tiled version is deferred (see Decisions).

Implementation reference

Exact behaviour of the code as of 13 Oct 2026. Everything below is implemented.

Layout conventions

• NDView / NDArray are always C-order. Strides are computed from the shape (c_strides in NDView.hpp); the last index moves fastest; neither class accepts custom strides. sub_view(first, last) slices axis 0 and stays contiguous.
• Input batch: NDView<AxisType, 3> of shape (n_frames, rows, cols). With p = row × ncols + col and npix = rows × cols, frame f, pixel p is at f × npix + p.
• Device accumulator d_histogram: npix × nbins uint32_t, bin-major: count of (p, b) at b × npix + p.
• API output: shape (rows, cols, nbins) of StorageType, pixel-major: (p, b) at p × nbins + b. Same as the CPU classes.

Class template and members

template <StorageType = uint16_t, DeviceStorageType = uint32_t, AxisType = float, int BLOCK_SIZE = 128> class PixelHistogramCUDA, in aare::cuda. Non-copyable, non-movable. NUM_SLOTS = 2 (compile-time).

Constructor

• check_args(rows, cols, n_bins, xmin, xmax, batch_size) throws std::invalid_argument unless rows, cols, bins ≥ 1, xmax > xmin (also rejects NaN) and batch_size ≥ 1; returns rows.
• Body: create events in a loop; cudaMemsetAsync(d_histogram, 0, d_histogram.bytes(), compute_stream) (on the compute stream so every kernel is ordered after it; plain cudaMemset would run on the legacy default stream, which non-blocking streams don't wait for); allocate the frame buffers in a loop.

Input registration

• register_input_buffer(ptr, bytes): synchronise the copy stream → unregister the previous buffer explicitly (re-registering overlapping memory would fail with cudaErrorHostMemoryAlreadyRegistered) → HostRegistration(ptr, bytes).
• unregister_input_buffer(): synchronise the copy stream → unregister.
• HostRegistration::unregister() releases the pointer from its unique_ptr before calling cudaHostUnregister, so a failure is never retried by the destructor.
• Contract: the caller's dataset must outlive the registration. Unregistered (pageable) input still works, but copies don't overlap.

submit_batch(frames): exact sequence

1. Validate: shape(0) ≤ m_batch_size, shape(1) == nrows, shape(2) == ncols.
2. slot = m_next_batch++ % NUM_SLOTS; n_frames = shape(0).
3. Copy stream: cudaStreamWaitEvent(kernel_done[slot]).
4. Copy stream: cudaMemcpyAsync H2D of npix × n_frames × sizeof(AxisType) into d_frames_batch[slot].
5. Copy stream: cudaEventRecord(copy_done[slot]).
6. Compute stream: cudaStreamWaitEvent(copy_done[slot]).
7. Compute stream: fill_bin_major<<<ceil(npix / BLOCK_SIZE), BLOCK_SIZE>>>(frames, d_histogram, npix, nbins, xmin, xmax, scale, n_frames), then cudaGetLastError.
8. Compute stream: cudaEventRecord(kernel_done[slot]).

Invariants: the host never waits. In the first round, waits on events never recorded are no-ops. All kernels run in order on one stream, which is what makes plain increments safe. n_frames is passed by value, so a partial last batch is safe. fill_histogram(frames) cuts frames into sub_views of at most m_batch_size along axis 0 and submits each. flush() const synchronises copy, then compute.

fill_bin_major

• p = (size_t)blockIdx.x × blockDim.x + threadIdx.x; return if p ≥ npix.
• For f in [0, n_frames): v = frames[f × npix + p]; skip (continue) unless v ≥ xmin && v < xmax, written this way so NaN is skipped too; bin = min((size_t)((v − xmin) × scale), nbins − 1); ++hist[bin × npix + p]. No saturation check: the uint32_t counter can't realistically overflow, and values() caps.
• Matches PixelHistogramImpl::fill_flat_unchecked operation for operation. No FMA risk: there's no addition after the multiply.
• No atomics, no shared memory, no __syncthreads(). All indices in size_t; pointers const / __restrict__.

values(out) const: exact sequence

1. Check out.shape() == (nrows, ncols, nbins); otherwise throw std::invalid_argument naming the expected shape. The shape fully determines the layout (C-order), so nothing else needs checking.
2. tmp = make_async_buffer<StorageType>(npix × nbins, compute_stream).
3. bin_to_pixel_major<<<ceil(total / BLOCK_SIZE), BLOCK_SIZE, 0, compute>>>(d_histogram, tmp, npix, nbins, count_cap()), then cudaGetLastError.
4. cudaMemcpyAsync(out.data(), tmp, total × sizeof(StorageType), D2H, compute_stream).
5. cudaStreamSynchronize(compute_stream). tmp is freed with cudaFreeAsync when it goes out of scope, also if a CUDA_CHECK throws.

• Semantics: a snapshot of every batch submitted so far; the accumulator is never modified, so filling can continue afterwards. The caller needs no sync before calling.
• values() allocates an uninitialised NDArray (shape-only constructor), calls values(view) and returns it (NRVO, else a move).
• count_cap() (private, static constexpr): for integer StorageType, the smaller of its max and DeviceStorageType's max (compared as uintmax_t); for floating-point storage, DeviceStorageType's max. Capping once here equals the CPU's per-increment saturation, because counts only increase.
• bin_to_pixel_major: one thread per output element t (64-bit index, return if t ≥ npix × nbins); pix = t / nbins, bin = t − pix × nbins; out[t] = (StorageType)min(in[bin × npix + pix], cap) with a ternary. Writes coalesced, reads strided (served mostly from cache).

values_pixel(row, col) const

• A debug/convenience call for inspecting one pixel; for many pixels, call values() once. Throws std::invalid_argument unless row < nrows and col < ncols.
• cudaMemcpy2DAsync(dst = host DeviceStorageType[nbins], dpitch = sizeof(DeviceStorageType), src = d_histogram + p, spitch = npix × sizeof(DeviceStorageType), width = sizeof(DeviceStorageType), height = nbins, D2H, compute_stream), then synchronise.
• Cap each value at count_cap() and convert into the returned NDArray<StorageType, 1> in a host loop. (NDArray's expression constructor wouldn't do this: it needs an rvalue and truncates instead of capping.)
• Cost: the copy engine moves one 4-byte piece per bin, so per byte it's slow (2D copies are always slower than 1D unless width = pitch). Fine for a single pixel.

bin_edges() / bin_centers()

Computed on demand, nothing stored: edges[i] = xmin + (i × (xmax − xmin)) / nbins for i in [0, nbins] (same formula as PixelHistogramImpl); centers[i] = (edges[i] + edges[i+1]) / 2.

RAII helpers (include/aare/RAII_device.hpp)

All deleters are noexcept and ignore errors; CUDA_CHECK (which throws) is only used outside destructors.

Decisions

Settled

Open

• Input type (D1): frames are float today. Raw uint16_t would halve PCIe traffic, the likely bottleneck. For the performance pass.
• Batch size (D6): 64–256 frames should be enough (launch overhead ~5 µs vs ~13 µs per frame copy). Choose with nsys.
• Streaming from file: today the dataset must already be in host memory. If needed: a blocking submit, or a pinned staging ring that MultiThreadedFileReader writes into directly.

Deferred optimisations

• Tiled transpose, only if the transpose shows up next to the D2H in nsys.
• A window of hot bins in shared memory (e.g. 32 × uint16_t per pixel, centred per pixel), only if the fill kernel turns out to be the bottleneck.

Steps

Each step ends with a passing test. Specs are in the reference.

1. RAII wrappers Done

   All CUDA resources are owned by wrappers in RAII_device.hpp; the class has no destructor.

2. Host-side pipeline Written, untested

   Constructor, input registration, submit_batch, fill_histogram, flush(), bin edges and centres. Verified in step 6.

3. Finalize kernel + values() Done, verified

   Checked with a scratch nvcc test on the RTX 4090: a known pattern count = bin × npix + p on 37 × 53 px and 45 bins (counts up to ~88 k), with uint16_t and uint32_t storage, plus a 1 × 1 case and a run through values(); count_cap() checked at compile time for six types. Zero mismatches. Becomes a real Catch2 test in step 6.

4. values_pixel() Done, verified

   cudaMemcpy2DAsync gather, then cap and convert on the host. Scratch test on the known pattern (37 × 53 px, 45 bins, counts up to 88 244): values_pixel(r, c) equals values()(r, c, :) for the four corners and two interior pixels, with uint16_t and uint32_t storage. Zero mismatches.

5. Fill kernel fill_bin_major Written, untested

   Implemented to the spec in the reference; the class compiles. Verified by step 6's comparison with the CPU class.

6. Tests and overlap Next

   • CMake wiring: reuse the AARE_CUDA option and the test pattern already in this branch's base (see ClusterFinderCUDA_testimpl.cu / ClusterFinderCUDA.test.cpp). The kernel launch means only nvcc can compile files that include the header, so a .cu test implementation is called from the Catch2 .test.cpp.
   • The tests listed under Testing and measuring.
   • nsys shows H2D i+1 overlapping kernel i.

7. Performance pass To do

   One change at a time, behind the same tests: uint16_t input, batch size, block size, several frames loaded ahead per thread. Aim for throughput within ~10–20 % of measured H2D bandwidth. Once copy-bound, kernel speedups show up as nothing.

8. Integration To do

   Python binding (with py::keep_alive on register_input_buffer), docs in docs/src/, a RELEASE.md entry, the explicit source/test lists in CMake (see AGENTS.md). Done when ctest and pytest pass on a CUDA build and a non-CUDA build is unaffected.

Loose ends

Small items in the current files.

• hpp L202
  tidyThe "TODO: Needs transpose" above values(out) is done.
• hpp L240
  tidyA bare // TODO under the Option B comment, with nothing to do.
• hpp L173
  tidyThe "let's not use shmem at first" comment inside the launch can go: the decision is final.
• cuda_check.cuh
  tidyNo newline at the end of the file, so compilers warn "backslash-newline at end of file". The file comes from the cluster-finder branch: fix it there.

Testing and measuring

Tests

• Known pattern written into d_histogram (via FRIEND_TEST from tests/friend_test.hpp): checks the transpose and capping without the fill kernel. Use odd sizes so partial blocks are covered.
• Exact match with PixelHistogramImpl on random frames, including values at xmin, just below xmax, at xmax and outside the range.
• Pipeline cases: a partial last batch; several fill_histogram calls in a row; values() mid-run, then more filling.
• Race test (a test of the tests): temporarily remove the kernel_done wait and slow the kernel down. Results must become wrong; if they don't, the tests can't catch a synchronisation bug.

Measuring

• Ceiling: measure pinned H2D bandwidth (nvbandwidth or bandwidthTest) and divide by the frame size to get the maximum frames/s.
• Timeline: wrap submit, values() and file reading in NVTX ranges, then:

nsys profile --trace=cuda,nvtx -o pixhist ./build/your_bench
nsys stats --report cuda_gpu_trace pixhist.nsys-rep

• Gaps on the copy row while kernels run: compute-bound, so optimise the kernel.
• Gaps on both rows: the host is the bottleneck (file reading, a hidden sync).
• Copies labelled "Pageable": the input isn't pinned, so nothing overlaps.
• A memcpy or memset on the default stream: something escaped your streams.

Pitfalls

• Default-stream callsPlain cudaMemcpy and cudaMemset run on the legacy default stream, which streams created with cudaStreamNonBlocking don't wait for. Keep every operation on your own streams.
• Pageable memorycudaMemcpyAsync from unpinned memory works but is effectively synchronous. Nothing errors; you just lose the overlap.
• Throwing in destructorsCUDA_CHECK throws, and a throw from a destructor calls std::terminate. Deleters ignore errors instead.
• 32-bit index overflowblockIdx.x × blockDim.x is computed in 32 bits; widen to size_t first. bin × npix and frame × npix also need 64 bits on large detectors.
• Host constexpr in device codestd::numeric_limits<T>::max() needs nvcc's --expt-relaxed-constexpr in a kernel. That's why the cap is computed on the host and passed in.
• Kernel arguments vs pointersArguments are copied at launch, so passing the frame count by value is safe. Anything read through a pointer is read when the kernel runs.

Self-check

Answers are folded.