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

Delete detailed optimisation reports
Build on RHEL8 / build (push) Successful in 3m28s
Details
Build on RHEL9 / build (push) Successful in 3m39s
Details
Run tests using data on local RHEL8 / build (push) Successful in 4m48s
Details

Browse Source
This commit is contained in:
kferjaoui
2026-09-02 13:49:35 +02:00
parent 9921944772
commit 8567d3f00e
3 changed files with 0 additions and 696 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/cuda_optimization_recap.md
-257
View File
@@ -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<StreamContext> 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<COMPUTE_TYPE>(d_frame[gid])
- static_cast<COMPUTE_TYPE>(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<<<grid, block, smem, stream>>>(...);
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<float>(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
docs/optimization_summary.md
-202
View File
@@ -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 <pybind11/stl.h>` 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).
docs/pedestal_precision_f32_cancellation.md
-237
View File
@@ -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).