From 8567d3f00eb993e62faacd5b6945387737d19fd3 Mon Sep 17 00:00:00 2001 From: kferjaoui Date: Wed, 2 Sep 2026 13:49:35 +0200 Subject: [PATCH] Delete detailed optimisation reports --- docs/cuda_optimization_recap.md | 257 -------------------- docs/optimization_summary.md | 202 --------------- docs/pedestal_precision_f32_cancellation.md | 237 ------------------ 3 files changed, 696 deletions(-) delete mode 100644 docs/cuda_optimization_recap.md delete mode 100644 docs/optimization_summary.md delete mode 100644 docs/pedestal_precision_f32_cancellation.md diff --git a/docs/cuda_optimization_recap.md b/docs/cuda_optimization_recap.md deleted file mode 100644 index 9986fbde..00000000 --- a/docs/cuda_optimization_recap.md +++ /dev/null @@ -1,257 +0,0 @@ -# ClusterFinderCUDA — Optimization Recap - -A step-by-step walkthrough of the optimizations applied to `ClusterFinderCUDA`, -in the order they were implemented. Each section = one slide: what changed, -why, and what it bought us. - -**Setup** (constant across all measurements below): -- GPU: RTX 4090 (Ada, sm_89), PCIe 4.0 x16 -- Detector: 400×400, `uint16` frames -- Cluster size: 3×3, `nSigma = 5` -- Reference: single-threaded CPU `ClusterFinder` - -| Reference | ms / frame | -|-----------------------------|-----------:| -| CPU (single-threaded) | ~1.7 – 1.8 | - ---- - -## Step 0 — First CUDA port: single-frame `find_clusters()` - -The very first GPU path: one frame in, H2D copy → kernel → D2H copy, fully -synchronous, no streams. - -```cpp -cuda_cf.find_clusters(frame, frame_number); -``` - -| Metric | ms / frame | -|-------------------------------------------|-----------:| -| GPU, single frame, single stream | **0.057 – 0.062** | - -**Already ~28–30× faster than CPU**, just from moving the stencil to the GPU — -before any pipelining work begins. This number stays roughly constant -throughout the rest of the journey: it's the "no batching, no overlap" floor -that the batched path is measured against. - ---- - -## Step 1 — Multi-stream batched cluster finder - -*Commit `3ed773e` — Add multi-stream ClusterFinderCUDA with batched processing* - -- Introduced `StreamContext`: each CUDA stream owns its own device buffers - and pedestal arrays -- `find_clusters_batched()`: N frames distributed round-robin across - `n_streams` streams -- Split `ClusterFinderCUDA.cuh` into the host wrapper (`ClusterFinderCUDA.hpp`) - and the device kernel (`clusterfinder_kernel.cuh`) - -```cpp -struct StreamContext { - cudaStream_t stream; - FRAME_TYPE *d_frame; - ClusterType *d_clusters; - uint32_t *d_cluster_count; - PEDESTAL_TYPE *d_pd_mean, *d_pd_sum, *d_pd_sum2; -}; -std::vector m_streams; -``` - -This is the scaffolding everything else builds on — no isolated benchmark -yet, but it's what makes "more streams = more overlap" possible in Step 4. - ---- - -## Step 2 — Mixed precision: FP32 stencil, FP64 pedestal - -*Commit `ac96d1f` — Implement mixed precision: f32 stencil, f64 pedestal* - -```cpp -// clusterfinder_kernel.cuh -using COMPUTE_TYPE = float; // stencil arithmetic + shared memory - -auto load_pixel = [&] __device__(ssize_t gr, ssize_t gc) -> COMPUTE_TYPE { - auto gid = gc + ncols * gr; - return static_cast(d_frame[gid]) - - static_cast(d_pd_mean[gid]); -}; -``` - -**Why:** -- FP32 throughput on Ada is ~64× FP64 — the per-pixel stencil sum/max/quadrant - reductions don't need double precision -- `float` shared memory (stride-18 tile) maps to distinct banks; the previous - `double` layout caused bank conflicts - -Pedestal accumulation (`mean`/`sum`/`sum2`) stays `double` at this stage — -that's addressed in Step 5. No isolated number here; this change is a -prerequisite for the kernel speedup measured later. - ---- - -## Step 3 — Per-frame kernel timing via CUDA events - -*Commit `34e69a8` — Add per-frame kernel timing via CUDA events* - -```cpp -cudaEventRecord(m_kernel_start_pool[slot], stream); -find_clusters_in_single_frame<<>>(...); -cudaEventRecord(m_kernel_stop_pool[slot], stream); -... -float ms = 0; -cudaEventElapsedTime(&ms, start, stop); -``` - -This is purely **instrumentation** — it adds `avg_kernel_time_ms()`, which -isolates *kernel-only* time from PCIe transfer time. Every "kernel: X µs" -number quoted from here on comes from this pool. No performance change, but -without it Step 5's "3.3× kernel speedup" would be invisible (it's masked by -PCIe in the wall-clock total). - ---- - -## Step 4 — Eliminate sync barriers + pinned transfers - -*Commits `88e0e8d` (transfer/kernel hot path) + `6a12e3d` (refactor)* - -**88e0e8d:** -- Per-stream **pinned host staging buffers** for truly async H2D/D2H -- Stop reserving full device capacity per output frame — delay cluster - payload construction until a candidate is confirmed -- Replace per-pixel `sqrtf()` threshold checks with squared comparisons - -**6a12e3d:** -- Replace **one `cudaStreamSynchronize` per frame** with **one per stream per - batch** — cuts sync calls from `O(n_frames × n_streams)` to `O(n_streams)` -- Unified D2H output layout `[uint32_t count | clusters[max]]` in one - lazily-allocated pinned pool -- New `register_input_buffer()` / `unregister_input_buffer()` — pin the - caller's batch buffer once via `cudaHostRegister`, so every - `find_clusters_batched()` slice transfers at DMA speed (~22 GB/s) instead of - ~15 GB/s for pageable memory - -```cpp -// continuous-queue launch loop (no per-frame sync) -for (auto &sc : m_streams) { - for (frame : frames_for_this_stream) { - cudaMemcpyAsync(sc.d_frame, ..., cudaMemcpyHostToDevice, sc.stream); - find_clusters_in_single_frame<<<..., sc.stream>>>(...); - cudaMemcpyAsync(h_output_pinned[slot], sc.d_clusters, ..., sc.stream); - } -} -for (auto &sc : m_streams) cudaStreamSynchronize(sc.stream); // once per stream -``` - -| Configuration | ms / frame | -|----------------------------------------|-----------:| -| Batched, 1 stream | 0.052 | -| Batched, 5 streams — **before** refactor | ~0.034 | -| Batched, 5 streams — **after** refactor | **0.028** | - -**−18%** at 5 streams, purely from removing host-side sync barriers and -making transfers DMA-speed. Going from 1→5 streams alone (52→28 µs) shows the -multi-stream scaffolding from Step 1 finally paying off once the host isn't -serializing on every frame. - ---- - -## Step 5 — FP32 device pedestal + bulk memcpy drain - -*Commit `4c66802` — FP32 device pedestal and bulk memcpy drain* - -```cpp -// device pedestal arrays: double -> float -float *__restrict__ d_pd_mean, *__restrict__ d_pd_sum, *__restrict__ d_pd_sum2; - -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); -``` - -```cpp -// D2H drain: per-cluster push_back loop -> single resize + memcpy -results[frame_idx].resize(n_found); -std::memcpy(results[frame_idx].data(), src, n_found * sizeof(ClusterType)); -``` - -| Metric | Before | After | Change | -|-------------------------------|-------:|-------:|-------:| -| Kernel-only time | 15 µs | 4.6 µs | **3.3×** | -| Batched, 5 streams (overall) | 28 µs | 26 µs | ~7% | - -**Why the kernel got 3.3× faster but the overall number barely moved:** -this is Amdahl's law in action. At 28 µs total, the kernel was already a -small slice — PCIe H2D/D2H dominates. Cutting the kernel to a quarter of its -size can only remove *the kernel's share* of the time. The big remaining -lever is PCIe, not compute — which motivates Step 6. - ---- - -## Step 6 — Async `submit_batch()` / `collect()` pipeline - -*Commit `5922c73` — async submit_batch/collect API* - -```cpp -auto tok = cf.submit_batch(buf_a, first_frame=0); // enqueue, don't wait -for (...) { - buf_b[:n] = data[start:start+n]; // CPU fills next buffer - auto next_tok = cf.submit_batch(buf_b, first_frame=start); - results += cf.collect(tok); // GPU runs buf_b meanwhile - tok = next_tok; - std::swap(buf_a, buf_b); -} -results += cf.collect(tok); -``` - -- `submit_batch()`: enqueues H2D + kernel + D2H, returns a `BatchToken`, - **never blocks** -- `collect()`: waits on a `cudaEvent` (not `cudaStreamSynchronize`) so a - second batch already queued behind the first keeps running - -| Metric | ms / frame | -|---------------------------------|-----------:| -| Batched, 5 streams (Step 5) | 0.026 | -| **Async pipeline** | **0.022** | - -**−15–18%.** The mechanism: for `BATCH_SIZE=2000` (640 MB/batch), the CPU -memcpy that fills the next buffer (~25 ms) is *longer* than the GPU's batch -execution time (~12 ms). The async pipeline hides the entire GPU batch inside -that CPU memcpy window — the GPU is no longer idle while the CPU prepares the -next chunk. - ---- - -## Summary — the full journey - -| Stage | ms / frame | Speedup vs CPU | -|------------------------------------------|-----------:|---------------:| -| CPU (reference) | 1.7 – 1.8 | 1× | -| GPU, single frame (Step 0) | 0.057 | ~30× | -| Batched, 1 stream (Step 1) | 0.052 | ~33× | -| Batched, 5 streams, sync-barrier removed (Step 4) | 0.028 | ~63× | -| + FP32 device pedestal (Step 5) | 0.026 | ~67× | -| + Async pipeline (Step 6) | **0.022** | **~78×** | - -``` -1.8 ms ──────────────────────────────────────────────────► CPU -0.057 ms ──────► GPU naive -0.052 ms ─────► batched, 1 stream -0.028 ms ───► sync-barrier removal + pinned transfers -0.026 ms ───► FP32 pedestal (kernel 3.3×, overall ~flat — Amdahl) -0.022 ms ──► async submit/collect pipeline -``` - ---- - -## What's next - -See `optimization_summary.md` for the original profiling-driven roadmap. -Remaining levers, roughly by expected impact: - -1. **GPU-side cluster compaction before D2H** — D2H currently moves - `max_clusters_per_frame` slots regardless of occupancy -2. **Larger batch size** — amortize per-batch launch overhead further -3. **Smaller `Cluster` struct** — reduce D2H payload per cluster -4. **GPUDirect RDMA** — eliminate H2D entirely if data arrives via - InfiniBand/RoCE diff --git a/docs/optimization_summary.md b/docs/optimization_summary.md deleted file mode 100644 index fc91ada4..00000000 --- a/docs/optimization_summary.md +++ /dev/null @@ -1,202 +0,0 @@ -# ClusterFinder CUDA — Profiling, Bottleneck Analysis, and Optimization Roadmap - -## 1. Setup - -- **Hardware**: RTX 4090 (Ada, sm_89), PCIe 4.0 x16 -- **Workload**: 50,000 frames of 400×400 uint16 detector data, 3×3 cluster finding -- **Tooling**: Nsight Systems 2024.5.1 (system-wide timeline), Nsight Compute 2024.1.1 (per-kernel deep dive) -- **Profiler counter access** unlocked via: - ``` - /etc/modprobe.d/nvidia-profiler.conf: - options nvidia NVreg_RestrictProfilingToAdminUsers=0 - ``` - Confirmed by `cat /proc/driver/nvidia/params | grep -i profil` → `RmProfilingAdminOnly: 0`. - ---- - -## 2. Bottlenecks Identified from Nsys Profiling - -### 2.1 Overall wall-time vs. ideal floor - -From the original Phase 4 run (50k frames, 5 streams, BATCH_SIZE=2000): - -| Resource | Total time | Phase 4 utilization | -|---|---|---| -| Kernel execution | 824 ms | 21.1% | -| H2D copy engine | 746 ms | 19.1% | -| D2H copy engine | 235 ms | 6.0% | -| **Phase 4 wall time** | **3903 ms** | — | - -These three engines are independent hardware resources. With perfect overlap, the wall-clock floor is `max(824, 746, 235) ≈ 824 ms`. Actual measured wall time: 3903 ms. **The GPU was idle ~80% of the time.** This is host-scheduling-limited, not GPU-compute-limited. - -### 2.2 Specific causes - -**Excessive host-device synchronization.** 110,029 `cudaStreamSynchronize` calls (~2.2 per frame) consuming 452 ms of host time. Each sync blocks the host scheduler thread and prevents queueing of subsequent stream commands. - -**Two D2H transfers per frame.** 110,002 D2H operations vs. 55,001 kernel launches — a count-then-data pattern where the second D2H size depends on the first D2H's result, forcing a per-frame sync. - -**Round-barrier pattern in the launch loop.** Code structure `for round { launch 5 streams; sync 5 streams }` forces stream 0 to wait for the slowest stream in each round. PCIe RX throughput visibly pulses with ~80 µs busy / ~40 µs idle in the timeline, confirming the H2D engine isn't continuously fed. - -**Redundant host memcpy.** The call site copies into a pageable `batch_buffer`, then `find_clusters_batched` copies *again* into per-stream pinned buffers. ~1.25 GB of host-DRAM memcpy per batch, all on one CPU core, blocking the launch loop. This shows up in the timeline as ~50–70 ms gaps between bursts where the GPU is fully idle. - -**Per-frame launch overhead at small kernel size.** Each frame triggers ~18 µs of fixed CUDA API overhead (launch + 3 memcpys + 2 syncs), vs. only ~15 µs of actual kernel execution. The host scheduling cost dominates the kernel cost. - -### 2.3 What was working - -PCIe pinned-memory transfers measured at 23.7 GB/s sustained (~95% of PCIe 4.0 x16 practical peak), confirming `cudaMallocHost` allocations are truly pinned. The kernel itself is regular and well-behaved (4% stddev across 55k launches, no thermal throttling, no spills per `ptxas` output: `0 bytes spill stores, 0 bytes spill loads`). - ---- - -## 3. Three-Tier Optimization Plan - -### Tier 1 — Recover the missing 2 seconds (host-side fixes) - -Largest expected impact. Targets the gap between 824 ms (ideal floor) and 3903 ms (measured) by fixing host scheduling and synchronization patterns. - -- **Drop the round-barrier**: continuous queueing across all streams, with a single sync per stream at end-of-batch -- **Collapse count-then-data D2H into a single fixed-size D2H**: per-frame device output as one contiguous `[count][clusters[max]]` block, eliminating the per-frame sync dependency -- **Eliminate the redundant host memcpy**: ClusterFinder owns one large pinned staging buffer, populated by a single bulk memcpy at the top of `find_clusters_batched` -- **Expected gain**: Phase 4 wall time from 3903 ms → ~1000–1500 ms - -### Tier 2 — Per-launch overhead reduction - -After Tier 1 closes the host-scheduling gap, the next bottleneck becomes the per-launch CUDA API cost. - -- **True work coalescing**: one kernel processes N frames via `blockIdx.z`, amortizing the 2.5 µs `cudaLaunchKernel` cost over N×. Requires kernel signature change and a decision on how to handle per-frame pedestal updates (likely: move pedestal update to a separate sparse kernel) -- **CUDA Graphs**: capture the per-batch sequence (memset → H2D → kernel → D2H) and replay as a single command-stream submission -- **Expected gain**: ~1.5–2× over Tier 1, approaching the H2D copy engine floor (~13.5 µs/frame, ~675 ms for 50k frames) - -### Tier 3 — Kernel internals - -Worth doing only after Tier 1 and 2, and only after `ncu --set full` confirms which kernel-internal change matters. - -- **Drop the per-pixel double-precision arithmetic** for variance/threshold (Ada has 1/64 FP64 throughput; precompute thresholds as float during pedestal updates) -- **Pedestal update write-back traffic**: most pixels are non-photons and write 3 × 8 bytes per pixel, dominating DRAM traffic. Either move to float, or defer pedestal updates to a separate kernel -- **Halo-loading divergence cleanup**: replace the role-based corner-handling code with a uniform cooperative tile load -- **Block geometry**: A/B test 32×8 vs 16×16 once everything else is locked in - ---- - -## 4. What Was Implemented (Tier 1) - -### Changes to `ClusterFinderCUDA.hpp` - -**New constructor signature**: -```cpp -ClusterFinderCUDA(Shape<2> shape, COMPUTE_TYPE nSigma, - size_t max_clusters_per_frame, // was: capacity (global) - int n_streams, - size_t staging_batch_capacity) // new -``` -The third argument's semantics changed from "global cluster buffer size per stream" to "tight upper bound per frame." This is now a real number used for the fixed-size D2H, not an opaque global capacity. - -**Per-frame device output as single block**: -``` -d_output layout: [uint32_t count][padding to ClusterType alignment][ClusterType clusters[max]] -``` -Offset computed as `(count_bytes + cluster_align - 1) & ~(cluster_align - 1)` for power-of-two alignment safety. - -**Internal pinned staging buffer**: one `cudaMallocHost` of `staging_batch_capacity × image_bytes`, owned for the lifetime of the ClusterFinder. The caller's `NDView` is bulk-copied into this buffer once per batch. - -**Continuous-queue launch loop**: frames distributed round-robin across streams; per-stream nested loop queues all of a stream's frames back-to-back without intermediate syncs. Single `cudaStreamSynchronize` per stream at end of batch. - -**Lazy-grown pinned output pool**: one slot per frame within a batch (`n_frames × output_bytes_per_frame`), so each frame's D2H targets its own pinned destination and frames don't overwrite each other before parsing. - -**Per-frame kernel event timing pool** (restored after initial removal): one event pair per frame slot, allocated lazily up to `staging_batch_capacity`, enables accurate `avg_kernel_time_ms()` reporting without the aliasing problem that one-event-pair-per-stream would have in continuous-queue mode. - -### Changes to bindings (`bind_ClusterFinderCUDA.hpp`) - -- Constructor binding updated to the new 5-argument form with renamed `max_clusters_per_frame` and new `staging_batch_capacity` -- `steal_clusters()` no longer takes the `realloc_same_capacity` bool argument -- `steal_clusters()` lambda uses `std::move` to transfer ownership out of the internal reference (required because `ClusterVector` is move-only) -- Re-added `#include ` for automatic `std::vector` return conversion -- Added `py::array::c_style | py::array::forcecast` to all `py::array_t` parameters to enforce contiguity at the binding boundary - -### Changes to Python factory (`ClusterFinderCUDA` function) - -- Kwarg `capacity` → `max_clusters_per_frame` -- New kwarg `staging_batch_capacity` (defaults to 2000) -- Default `n_streams` raised from 1 to 4 -- Tighter default `max_clusters_per_frame=2048` (was unbounded-feel with `capacity=1024` that would have silently truncated) - -### Changes to test/benchmark code - -- Constructor call sites updated to new signature -- `steal_clusters(true)` → `steal_clusters()` (no bool argument) -- Optional: out-parameter overload of `find_clusters_batched` to avoid per-batch result vector allocation in tight benchmark loops - -### Changes to Jupyter notebook - -- Two `steal_clusters(realloc_same_capacity=...)` calls in the non-batched branch updated to `steal_clusters()` -- BATCHED branch unchanged - ---- - -## 5. What Still Needs to Be Done - -### Immediate next step: validate Tier 1 - -Re-run the test binary and check: -- Phase 4 wall time should drop from 3903 ms → ~1000–1500 ms -- `cudaStreamSynchronize` count in `nsys` stats should drop from 110k → ~`n_batches × n_streams` (~125 for 25 batches × 5 streams) -- `avg_kernel_time_ms()` from the restored event pool should report ~0.015 ms (15 µs) -- Per-frame wall time should land between 20–30 µs (vs. old 78 µs) - -If Phase 4 is significantly above ~1500 ms, the next thing to investigate in `nsys` is whether the host bulk-memcpy at the top of `find_clusters_batched` is on the critical path — visible as a gap at the start of each batch in the timeline. - -### Tier 1.5 — overlap host-memcpy with GPU execution - -The current implementation runs the bulk host-memcpy on the calling thread, blocking the launch loop. To hide this cost: - -- **Worker-thread approach**: dedicate one host thread to filling staging buffer N+1 while the GPU consumes batch N. Requires a ping-pong (or larger) staging buffer. -- **Ping-pong staging**: two pinned staging buffers, alternating between them per batch -- **Cost**: adds threading complexity, but the bulk memcpy (~50 ms per 2000-frame batch) is currently on the critical path; hiding it could close the remaining gap to the kernel floor - -### Tier 2 — Work coalescing and CUDA Graphs - -Pending Tier 1 validation. The two pieces: - -**Multi-frame kernel via `blockIdx.z`**: -- New kernel signature accepting `[n_frames, nrows, ncols]` input and per-frame output blocks -- Resolve the pedestal-update race: simplest option is to disable the in-kernel pedestal update and run a separate sparse update kernel once per batch (small mathematical difference from current frame-by-frame behavior, likely negligible for slowly-drifting pedestals) -- Estimated gain: collapses 2000 launches per batch into ~5–10 launches; saves ~120 ms of host overhead across Phase 4 - -**CUDA Graph capture**: once the per-batch sequence is regular (same operations on same buffers each time), capture once and replay. ~1.5–2× on launch-bound workloads. Synergizes with the multi-frame kernel because both reduce per-frame host overhead. - -### Tier 3 — Kernel internal optimization - -Requires `ncu --set full` profiling run on a single kernel launch to prioritize: -- Source counters → confirm 0 register spills -- Memory Workload Analysis → L1 hit rate (current 3×3 stencil should benefit from shared memory; check whether it actually does) -- Warp State Statistics → identify dominant stall reason (Long Scoreboard = DRAM-bound; Short Scoreboard = shared-memory-bound; Wait = barrier-bound) -- GPU Speed of Light / Roofline → confirm kernel is memory-bound (expected) and how far from the DRAM ceiling - -Then prioritize the kernel-internal changes (FP64 reduction, pedestal update offload, halo-load cleanup, block geometry) based on which one the profiler indicates as binding. - -### Open correctness issue (independent of perf work) - -The original Phase 5 comparison reported **13.7M `data_mismatch` clusters and ~0.5M each of CPU-only / GPU-only** out of 69M total. This is a kernel-vs-CPU algorithmic difference, not a performance issue, but worth tracing before further kernel optimization. Likely candidates: - -- Tie-breaking rule between adjacent local maxima (CPU and GPU may disagree on which pixel "wins" when two neighbors have equal values) -- Off-by-one in cluster window indexing at frame boundaries -- Difference in whether the central pixel must be a strict max (`>`) vs. ≥ max (`>=`) -- The kernel's pedestal-update timing (CPU updates per-pixel sequentially as it walks the frame; GPU updates the entire frame's pedestal in parallel after the cluster decisions) - -Profile a correct kernel, not a fast wrong one. - ---- - -## 6. Reference Numbers - -For future "did we improve?" comparisons: - -| Metric | Original | Tier 1 target | Tier 2 target | Theoretical floor | -|-------------------------------|-----------|---------------|---------------|--------------------| -| Phase 4 wall time | 3903 ms | ~1200 ms | ~700 ms | 675 ms (H2D-bound) | -| Per-frame wall | 78 µs | 24 µs | 14 µs | 13.5 µs | -| Per-frame kernel | 15 µs | 15 µs | 15 µs | 15 µs | -| `cudaStreamSynchronize` count | 110,029 | ~125 | ~25 | — | -| `cudaLaunchKernel` count | 55,001 | 55,001 | ~125 | — | -| H2D bandwidth | 23.7 GB/s | 23.7 GB/s | ~25 GB/s | 26 GB/s (PCIe 4.0 x16 pinned) | - -The H2D engine is the hard floor on this hardware. To go below 13.5 µs/frame requires either reducing PCIe payload (uint8 input → 2× faster) or doing more compute per frame transferred (currently the simplest finder; nothing else to amortize over). diff --git a/docs/pedestal_precision_f32_cancellation.md b/docs/pedestal_precision_f32_cancellation.md deleted file mode 100644 index a33f43d8..00000000 --- a/docs/pedestal_precision_f32_cancellation.md +++ /dev/null @@ -1,237 +0,0 @@ -# CUDA ClusterFinder — why an f32 pedestal breaks (catastrophic cancellation) - -## TL;DR - -The CUDA kernel computes the per-pixel pedestal variance as -`var = E[X²] − E[X]²` in `DEVICE_PED_TYPE`. Pedestals sit at **~4655 ADU**, so it -recovers a **~2025** variance by subtracting two **~2.17×10⁷** numbers. In -**float32** (~7 significant digits) that subtraction carries a fixed **±3 ADU²** -absolute error: negligible for noisy pixels, but it drives low-noise pixels' -`rms²` to zero (via the clamp on the next line), turning ~1–2 % of pixels into -hot pixels that fire every frame → **+28 % spurious clusters and an unphysical -high-energy tail**. **float64** has the digits to spare, so it never happens -(CPU↔CUDA agree to 0.0001 %). - -**Fix:** either ship `DEVICE_PED_TYPE = double`, or keep f32 and remove the -cancellation (per-pixel offset accumulation, or Welford online variance). - ---- - -## 1. What was observed - -Two builds of `aare`, same data (MOENCH 400×400, Cu XRF, 1000 pedestal frames, -20k–50k data frames, `n_sigma = 5`): - -| build | `COMPUTE_TYPE` | `DEVICE_PED_TYPE` | CPU↔CUDA agreement | -|---|---|---|---| -| correctness | double (f64) | **double (f64)** | ~100 % (35 / 46,478,563 clusters = **0.0001 %**) | -| performance | double (f64) | **float (f32)** | CPU 116,010,113 vs CUDA 148,559,598 = **+28.06 %** | - -The f32-pedestal build also shows a large **unphysical tail** in the -cluster-energy spectrum. The *compute* precision (f64 in both) is not the -variable — the **pedestal** precision is. - ---- - -## 2. What the kernel computes - -Per pixel, in `include/aare/clusterfinder_kernel.cuh` (~lines 186–191): - -```cpp -DEVICE_PED_TYPE var_px = - d_pd_sum2[global_tid] / n - mean_px * mean_px; // variance = E[X²] − E[X]² -COMPUTE_TYPE rms_sq = (var_px > 0) ? var_px : 0; // clamp negatives to 0 -COMPUTE_TYPE nSig_sq_rms_sq = nSigma * nSigma * rms_sq; // squared threshold -``` - -Decision (Test 1): a pixel is a photon if `max_val² > nSig_sq_rms_sq`, i.e. -`max_val > nSigma · rms`. Everything hinges on `rms_sq` being right. - -The stored pedestal accumulators (in `DEVICE_PED_TYPE`), at equilibrium: - -- `d_pd_sum[tid] ≈ n · mean` -- `d_pd_sum2[tid] ≈ n · E[X²]` - ---- - -## 3. The magnitudes (with `mean = 4655`, `rms = 45`, `n = 1000`) - -| quantity | formula | value | -|---|---|---| -| `mean` | — | 4 655 | -| `mean²` | 4655² | **21 669 025** ≈ 2.17×10⁷ | -| `E[X²]` | `mean² + var` = 21 669 025 + 2025 | **21 671 050** ≈ 2.17×10⁷ | -| `var` | `E[X²] − mean²` | **2 025** | -| `sum2` | `n · E[X²]` = 1000 × 21 671 050 | **2.167×10¹⁰** | - -Critical ratio: - -``` -var / mean² = 2025 / 21 669 025 ≈ 9×10⁻⁵ -``` - -**We recover a number (2025) ~11 000× smaller than the two numbers subtracted to -get it.** That is the entire problem. - ---- - -## 4. float32 precision - -IEEE-754 single = 24 significant bits: - -- **≈ 7.2 significant decimal digits** (24 · log₁₀2). -- **ULP rule:** a value in `[2ᵉ, 2ᵉ⁺¹)` lies on a grid of spacing `ULP = 2ᵉ⁻²³`. - Absolute precision therefore *grows* with magnitude. - -Applied to each quantity (find `e` with `2ᵉ ≤ value < 2ᵉ⁺¹`): - -| quantity | value | range | e | ULP = 2ᵉ⁻²³ | -|---|---|---|---|---| -| `mean` | 4 655 | [2¹², 2¹³) | 12 | 2⁻¹¹ ≈ **0.0005** | -| `mean²` | 2.17×10⁷ | [2²⁴, 2²⁵) | 24 | 2¹ = **2** | -| `E[X²]` = `sum2/n` | 2.17×10⁷ | [2²⁴, 2²⁵) | 24 | **2** | -| `sum2` | 2.17×10¹⁰ | [2³⁴, 2³⁵) | 34 | 2¹¹ = **2048** | - -Key readings: - -- `mean` in f32 is essentially perfect (±0.0005 ADU) → **pedestal subtraction is - not the problem.** -- `mean²` and `E[X²]` are representable only to the **nearest 2 ADU²**. -- The raw accumulator `sum2` is good only to the **nearest 2048** — its bottom - ~4 decimal digits are noise (7-digit float holding an 11-digit number). - ---- - -## 5. Error propagation through the subtraction - -``` -var = E[X²] − mean² - = (21 671 050 ± 2) − (21 669 025 ± 2) ← each operand on a grid of 2 - = 2025 ± ~3 ← plus sum2/n quantization (±2) → ±3–4 ADU² -``` - -The result 2025 is finely representable, **but carries a ±3–4 ADU² error -inherited from the two giant operands**, and that error is a **fixed absolute -size — it does not shrink for quiet pixels.** - ---- - -## 6. Absolute error → relative error explodes for quiet pixels - -| pixel | true rms | true var | var in f32 | rms error | 5σ threshold effect | -|---|---|---|---|---|---| -| noisy (bulk) | 45 | 2025 | 2025 ± 3 | 0.07 % | 225.0 → 225.1 — harmless | -| moderate | 10 | 100 | 100 ± 3 | 1.5 % | 50 → 50.8 | -| quiet | 5 | 25 | 25 ± 3 | 6 % | 25 → 26.5 | -| very quiet | 3 | 9 | 9 ± 3 | 17 % | 15 → 12–17 | -| near-flat | 2 | 4 | **4 ± 3 → can be ≤ 0** | **∞** | **clamp → 0** | - -The bulk of the detector is fine (why the spectrum mostly looks right), but every -low-noise pixel has a corrupted threshold, and the flattest fall off the cliff. - ---- - -## 7. The clamp turns "wrong" into "catastrophic" - -`rms_sq = (var_px > 0) ? var_px : 0` - -When cancellation drives a quiet pixel's `var_px` to **≤ 0**, it clamps to -`rms_sq = 0` → `nSig_sq_rms_sq = 0` → Test 1 becomes `max_val² > 0`. The pixel is -declared a photon on **any** positive fluctuation: its 5σ gate (~225 ADU) has -become a 0σ gate, so it fires essentially every frame. (The clamp is not a bug — -it exists *because* negative variances were already occurring; it just converts -them into hot pixels rather than NaNs.) - ---- - -## 8. Order-of-magnitude count check (the +28 %) - -Measured excess: **+32.5M** clusters over **20 000** frames: - -``` -32.5×10⁶ / 20 000 ≈ 1600 extra clusters per frame -``` - -A stuck (threshold-0) pixel yields ~1 cluster/frame (whenever it is the local max -of its window). So this matches **~1600–3000 corrupted pixels out of 160 000 — -about 1–2 % of the array** having low enough noise to be underestimated/clamped -by the ±3 ADU² f32 error. That accounts for the whole 28 %. - ---- - -## 9. Why the tail is *high*-energy - -A correct 5σ threshold (~225 ADU) removes everything below it — the spectrum -starts at ~225. A stuck pixel with threshold ≈ 0 admits the **entire positive -side** of its distribution. Its cluster energy is the window sum -`Σ(raw − mean)`, and with the gate wide open it accepts ordinary positive noise -(low energy) **plus** the occasional large excursion or real charge drifting -through that a proper cut would remove. Instead of a clean edge at 225 you get a -population smeared from ~0 upward — an unphysical **tail**. Same corrupted pixels -produce both the count excess and the tail. - ---- - -## 10. Compounding effect: the running update - -The per-frame EMA update (~lines 293–294) runs in the same f32: - -```cpp -sum2 += raw*raw − sum2/n; -``` - -`raw² ≈ 21.67M` (ULP 2); the increment fluctuates by `±2·mean·δ ≈ ±2·4655·45 ≈ -±419 000`; it is added to `sum2 ≈ 2.17×10¹⁰` (ULP 2048). Each update therefore -rounds by up to ±1024 → injects ~±1 ADU² of jitter into `E[X²]` every frame, on -top of the static cancellation error. Over thousands of frames the EMA stays -bounded but permanently noisy, nudging borderline pixels across the clamp. - -Also relevant: `sync_pedestal_to_device` (in `ClusterFinderCUDA.hpp`) casts the -host's **double** `sum2` **down to f32 on upload**, so in the f32 build the -precision is destroyed *before the first data frame*. - ---- - -## 11. Why double precision makes it vanish - -Same formula, `DEVICE_PED_TYPE = double` → 52 mantissa bits ≈ 15–16 decimal -digits. ULP at `2.17×10⁷` is `2²⁴⁻⁵² = 2⁻²⁸ ≈ 4×10⁻⁹ ADU²`: - -``` -var = (21 671 050 ± 4e-9) − (21 669 025 ± 4e-9) = 2025 ± 1e-8 -``` - -Variance is exact for all practical purposes, no pixel is ever clamped, and the -host→device cast is a no-op → CPU and CUDA agree to 0.0001 %. - ---- - -## 12. Fix options - -| option | what | pro | con | -|---|---|---|---| -| **A. f64 pedestal** | `DEVICE_PED_TYPE = double` | trivially correct (proven) | 2× pedestal memory + bandwidth; loses the f32 perf win | -| **B1. per-pixel offset** | accumulate `X − X0` with `X0 ≈ round(mean)` per pixel; keep `sum`/`sum2` of the *centered* values | keeps f32 bandwidth **and** correctness; small kernel + host-sync change | needs a per-pixel offset array and one-time rebase | -| **B2. Welford** | store running `mean` + `M2` directly; never form `E[X²]` | numerically ideal; no offset bookkeeping | more work per update; larger change to the update step | - -**Why B works:** with an offset `X0 ≈ 4655`, the accumulated quantities become -O(45) and O(2025) instead of O(4655) and O(2×10⁷), so `var = sum2'/n − mean'²` -is a difference of ~2000-scale numbers → f32 resolves it to ~1e-4 relative. The -cancellation is gone and the variance no longer depends on the absolute pedestal -level. - -**Recommendation:** **B1 (per-pixel offset)** — recovers the f32 performance -without the errors. B2 (Welford) is the cleanest if a larger refactor of the -update path is acceptable. - ---- - -## References - -- Kernel: `include/aare/clusterfinder_kernel.cuh` — variance ~L186–191, clamp - L190, running update ~L293–298. -- Host sync / f32 cast: `include/aare/ClusterFinderCUDA.hpp` — - `sync_pedestal_to_device`. -- Validation: `python/tests/ClusterFinderFrozen_vs_CUDA.ipynb` (f64 build, - ~100 % concordance) and `python/tests/ClusterFinderCUDA_perf.ipynb` (f32 - pedestal build, +28 % + tail).