mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2026-09-04 23:30:55 +02:00
Benchmarking of CUDA cluster finder
This commit is contained in:
kferjaoui
@@ -0,0 +1 @@
|
||||
,ferjao_k,pc-moench-04.psi.ch,11.08.2026 11:46,file:///home/ferjao_k/.config/libreoffice/4;
|
||||
Binary file not shown.
@@ -0,0 +1,467 @@
|
||||
# ClusterFinderCUDA — Benchmark Results, opt1 → opt6
|
||||
|
||||
Consolidated, verified performance numbers for the CUDA cluster-finder optimization
|
||||
ladder. Every number below is traced to the code, notebook cell, or profiler report
|
||||
that produced it, and each is tagged **quotable** or **not quotable** with the reason.
|
||||
|
||||
Companion documents:
|
||||
- `docs/pedestal_precision_f32_cancellation.md` — why the naive f32 pedestal failed and how B1 fixes it
|
||||
- `docs/cuda_optimization_recap.md`, `docs/optimization_summary.md` — earlier narrative/roadmap notes
|
||||
- `docs/ClusterFinderCUDA_optimizations.pptx` — the deck these numbers feed
|
||||
|
||||
---
|
||||
|
||||
## 1. Environment
|
||||
|
||||
| item | value |
|
||||
|---|---|
|
||||
| GPU | NVIDIA GeForce RTX 4090 (Ada, sm_89), 24 GB, driver 595.71.05 |
|
||||
| GPU clocks | idle 210 MHz → boost 3120 MHz; **persistence mode disabled** |
|
||||
| FP64 rate | 1/64 of FP32 on this part (relevant to opt6) |
|
||||
| CPU | AMD Ryzen 9 7950X, 16 cores / 32 threads |
|
||||
| RAM | 125 GiB, **no swap** |
|
||||
| CUDA | 12.4 (nvcc V12.4.131) |
|
||||
| Profiler | Nsight Systems 2024.5.1 (`/opt/nvidia/nsight-systems/2024.5.1`) |
|
||||
| Host | `pc-moench-04` |
|
||||
| Branch | `bench/opt2-pipeline` (off `feature/cuda_clusterfinder` @ `ce256dd`) |
|
||||
|
||||
**Dataset** — MOENCH, MAX IV beamtime, Cu fluorescence:
|
||||
|
||||
```
|
||||
/mnt/sls_det_storage/moench_data/2603_MaxIVBeamtime/2026032408/process/xrf/
|
||||
Cu_factor_10_data_master_0.json (100 000 frames, 400×400 uint16)
|
||||
Cu_factor_10_pedestal_master_0.json (1 000 frames used for pedestal training)
|
||||
```
|
||||
|
||||
Frame size 400×400×2 B = 320 000 B (312.5 KiB). All finders trained on the
|
||||
**same 1 000 pedestal frames**;
|
||||
`n_sigma = 5` throughout. Data pre-loaded into RAM with `read_n()` so file I/O is
|
||||
outside every timing loop.
|
||||
|
||||
---
|
||||
|
||||
## 2. The optimization ladder
|
||||
|
||||
| step | what changed | how it is measured |
|
||||
|---|---|---|
|
||||
| **baseline** | `ClusterFinderMT`, 48 threads | `ClusterFinderMT(..., n_threads=48)` + `ClusterCollector` |
|
||||
| **opt1** | first CUDA port: 1 stream, one launch per frame, no batching | `ClusterFinderCUDAOpt2(..., n_streams=1)` + `find_clusters()` per frame |
|
||||
| **opt2** | multi-stream scaffolding + host-side batching (bulk memcpy) | `ClusterFinderCUDAOpt2(..., n_streams=4)` + `find_clusters_batched()`, batch 2000 |
|
||||
| **opt3** | pipeline rework: remove per-round sync barriers, fixed-size D2H | `ClusterFinderCUDA(..., n_streams=4)` + `find_clusters_batched()`, **no pinning** |
|
||||
| **opt4** | DMA-speed transfers via pinned host input | opt3 + `register_input_buffer(data)` |
|
||||
| **opt5** | CUDA Graphs (pre-recorded H2D→kernel→D2H per stream) | `ClusterFinderCUDAGraph(...)` + pinned input |
|
||||
| **opt6** | **first kernel optimization**: f32 device pedestal + variance rewrite (B1) | rebuild with `DEVICE_PED_TYPE = float` |
|
||||
|
||||
opt1–opt5 are **pipeline/host-side**; opt6 is the first change to the **kernel itself**.
|
||||
|
||||
### Code behind each step
|
||||
|
||||
| step | primary source |
|
||||
|---|---|
|
||||
| opt1, opt2 | [`include/aare/ClusterFinderCUDAOpt2.hpp`](../include/aare/ClusterFinderCUDAOpt2.hpp), [`include/aare/clusterfinder_kernel_opt2.cuh`](../include/aare/clusterfinder_kernel_opt2.cuh) — snapshot of commit `88e0e8d` (pre-refactor pipeline), namespace `aare::device_opt2` |
|
||||
| opt3, opt4, opt6 | [`include/aare/ClusterFinderCUDA.hpp`](../include/aare/ClusterFinderCUDA.hpp), [`include/aare/clusterfinder_kernel.cuh`](../include/aare/clusterfinder_kernel.cuh) |
|
||||
| opt5 | [`include/aare/ClusterFinderCUDA_graph.hpp`](../include/aare/ClusterFinderCUDA_graph.hpp) |
|
||||
| bindings | [`python/src/bind_ClusterFinderCUDAOpt2.hpp`](../python/src/bind_ClusterFinderCUDAOpt2.hpp), [`python/src/cuda_bindings.cu`](../python/src/cuda_bindings.cu) |
|
||||
| factories | [`python/aare/ClusterFinder.py`](../python/aare/ClusterFinder.py) |
|
||||
|
||||
Relevant history: `3ed773e` (multi-stream+batched) → `ac96d1f` (mixed precision) →
|
||||
`88e0e8d` (**opt1/opt2 snapshot**) → `6a12e3d` (pipeline refactor = opt3) →
|
||||
`4c66802` (FP32 pedestal, introduced the tail) → `5922c73` (async API) →
|
||||
`1bf317f` (local-max fix) → `a42d71c` (graphs = opt5).
|
||||
|
||||
### Precision configuration
|
||||
|
||||
Two type aliases in [`clusterfinder_kernel.cuh:16-17`](../include/aare/clusterfinder_kernel.cuh#L16-L17):
|
||||
|
||||
```cpp
|
||||
using COMPUTE_TYPE = float; // stencil arithmetic — float in ALL builds below
|
||||
using DEVICE_PED_TYPE = double; // device pedestal — the opt6 knob (double → float)
|
||||
```
|
||||
|
||||
- **"f64 build"** in this document = `COMPUTE_TYPE=float`, `DEVICE_PED_TYPE=double` (mixed precision).
|
||||
- **"f32 build"** (opt6) = `COMPUTE_TYPE=float`, `DEVICE_PED_TYPE=float` (100 % f32).
|
||||
|
||||
> ⚠️ `ClusterFinderCUDAOpt2` is templated on `PEDESTAL_TYPE`, which its binding pins
|
||||
> to `double` ([`bind_ClusterFinderCUDAOpt2.hpp:14`](../python/src/bind_ClusterFinderCUDAOpt2.hpp#L14)).
|
||||
> **opt1 and opt2 are therefore unaffected by the opt6 flip** — their numbers are
|
||||
> identical in both campaigns by construction. Only opt3/opt4/opt5 respond to opt6.
|
||||
|
||||
---
|
||||
|
||||
## 3. Methodology — what is measurable and what is not
|
||||
|
||||
Three measurement artifacts were identified and controlled. **All three matter for
|
||||
how numbers may be quoted on a slide.**
|
||||
|
||||
### 3.1 First-run soft page faults (dominant, ~1.5–4 s per cell)
|
||||
|
||||
Each timed cell materializes ~233 M clusters ≈ **10 GB of host result heap**. On the
|
||||
first execution in a process, every 4 kB page must be faulted in and zeroed by the OS
|
||||
(millions of minor faults). Measured cost on this machine: **0.7 µs/fault**.
|
||||
|
||||
Instrumentation added to all seven timed notebook cells:
|
||||
|
||||
```python
|
||||
import resource
|
||||
def _faults():
|
||||
r = resource.getrusage(resource.RUSAGE_SELF)
|
||||
return r.ru_minflt, r.ru_majflt
|
||||
```
|
||||
bracketing the timed region, printing `minor faults: N (~N×0.7µs est.)`.
|
||||
|
||||
**Protocol: quote the run where the fault counter has plateaued** (< ~200 k). Cold
|
||||
numbers are inflated by 25–45 %. Verified by direct correlation (§9).
|
||||
|
||||
### 3.2 CUDA-event `kernel_ms` inflates under multi-stream saturation
|
||||
|
||||
`avg_kernel_time_ms()` uses a CUDA event pair on the kernel's own stream. It measures
|
||||
**elapsed time on that stream's timeline**, which includes queue-wait when other
|
||||
streams are competing for SMs. Consequences:
|
||||
|
||||
- In transfer-paced regimes: `event ≈ true kernel + ~5–7 µs` launch gap → usable.
|
||||
- In kernel-saturated regimes: **inflated up to 3.5×** → **not quotable**.
|
||||
- Symptom of saturation: the derived `PCIe + overhead = wall/N − kernel_ms` **goes
|
||||
negative** (kernels overlap, so wall/frame < kernel/frame).
|
||||
|
||||
Ground truth requires Nsight Systems (§8).
|
||||
|
||||
### 3.3 Profiler distorts wall clock
|
||||
|
||||
Under `nsys`, wall time per frame is ~4× the unprofiled value (API tracing overhead).
|
||||
**Take per-operation GPU times from nsys; take wall times from unprofiled runs.** The
|
||||
standalone probe script also pays the full first-touch fault tax inside its single
|
||||
timed call (fresh process, no warm pass), so **its wall times are not throughput
|
||||
numbers** either.
|
||||
|
||||
### 3.4 Other controls
|
||||
|
||||
- GPU clock ramp (210 MHz → 3.1 GHz) is real but < 0.1 % of a multi-second run; the
|
||||
invariance of `kernel_ms` across loaded/idle runs confirms it is not a factor.
|
||||
- The `cf_cuda_v1` notebook cell builds a histogram **inside** its timed loop
|
||||
(~4.5 s over 100 k frames). It is a reference row, **not part of the arc**.
|
||||
- `ClusterFinderMT` cannot restart after `stop()`; the CPU baseline is therefore
|
||||
always a first-pass number and carries its own ~1.8 s of allocator faults.
|
||||
|
||||
---
|
||||
|
||||
## 4. Campaign A — 3×3 pipeline arc (f64 pedestal) ★ headline
|
||||
|
||||
**Config**: `cluster_size=(3,3)`, `N=100 000`, `n_frames_pd=1000`, `n_sigma=5`,
|
||||
`BATCH_SIZE=2000`, `n_streams=4`, `max_clusters_per_frame=3000`.
|
||||
**Source**: [`python/tests/ClusterFinderCUDA_perf.ipynb`](../python/tests/ClusterFinderCUDA_perf.ipynb),
|
||||
warm pass (steady state, faults plateaued).
|
||||
|
||||
| step | variant | wall [s] | FPS | µs/frame | kernel [µs/fr] | host ovhd [µs/fr] | **vs CPU** | step gain | minor faults |
|
||||
|---|---|--:|--:|--:|--:|--:|--:|--:|--:|
|
||||
| baseline | CPU MT (48 threads) | 21.00 | 4 761 | 210.0 | — | — | 1.00× | — | 2.63 M |
|
||||
| **opt1** | 1 stream, per-frame | 6.68 | 14 968 | 66.8 | 23 | 44 | **3.14×** | 3.14× | 1 |
|
||||
| **opt2** | 4 streams + batching | 4.32 | 23 134 | 43.2 | 22 | 21 | **4.86×** | 1.55× | 64 931 |
|
||||
| **opt3** | pipeline rework, no pin | 3.76 | 26 588 | 37.6 | 24 | 14 | **5.58×** | 1.15× | 59 814 |
|
||||
| **opt4** | + pinned input (DMA) | 2.72 | 36 810 | 27.2 | 24 | 3 | **7.73×** | 1.38× | 1 027 |
|
||||
| **opt5** | + CUDA Graph | 2.53 | 39 472 | 25.3 | 25 * | — | **8.29×** | 1.07× | 1 032 |
|
||||
|
||||
\* graph reports kernel+PCIe+overhead combined; the kernel is not timed separately.
|
||||
|
||||
Reference row (not in the arc): current finder, per-frame, histogram inside the timed
|
||||
loop — 11.35 s / 8 813 FPS / 1.85×.
|
||||
|
||||
**Status: quotable.** All CUDA rows at steady state; CPU baseline is first-pass by
|
||||
necessity (see §3.4) — if a fault-corrected CPU (≈19.2 s) is preferred, all speedups
|
||||
shrink by ~9 % (e.g. opt5 8.29× → 7.57×). Be consistent and state the choice.
|
||||
|
||||
### The story in one row — host overhead
|
||||
|
||||
```
|
||||
opt1 → opt2 → opt3 → opt4 host overhead per frame
|
||||
44 → 21 → 14 → 3 µs (kernel fixed at ~23 µs)
|
||||
```
|
||||
|
||||
opt4's 3 µs means H2D/kernel/D2H are essentially fully overlapped; opt5 then trims
|
||||
residual launch overhead. **After opt4 the 3×3 pipeline is host/PCIe-bound, not
|
||||
kernel-bound** — which is exactly why opt6 does nothing here (§5).
|
||||
|
||||
---
|
||||
|
||||
## 5. Campaign B — opt6 at 3×3 (100 % f32)
|
||||
|
||||
Same config as Campaign A; rebuilt with `DEVICE_PED_TYPE = float`.
|
||||
**Source**: same notebook, warm pass.
|
||||
|
||||
| step | wall f64 → f32 [s] | FPS f64 → f32 | kernel f64 → f32 [µs] | Δ wall |
|
||||
|---|--:|--:|--:|--:|
|
||||
| CPU baseline | 21.00 → 21.45 | 4 761 → 4 662 | — | (run-to-run) |
|
||||
| v1 per-frame (serialized) | 11.35 → 10.57 | 8 813 → 9 465 | **24 → 13** | **−7 %** ✓ |
|
||||
| opt3 (no pin) | 3.76 → 3.88 | 26 588 → 25 798 | **24 → 14** | +3 % (noise) |
|
||||
| opt4 (pinned) | 2.72 → 2.83 | 36 810 → 35 372 | **24 → 16** | +4 % (noise) |
|
||||
| opt5 (graph) | 2.53 → 2.53 | 39 472 → 39 503 | — | **0 %** |
|
||||
| opt1 / opt2 | 6.68 → 6.73 / 4.32 → 4.88 | — | 23 / 22 unchanged | n/a — f64-pinned class |
|
||||
|
||||
**Result: at 3×3, opt6 halves the kernel and buys nothing end-to-end.** The kernel
|
||||
(13–16 µs) sits *below* the ~25–27 µs/frame transfer+host floor, so it is entirely
|
||||
hidden by stream overlap. The control is the serialized `v1` path, where the kernel
|
||||
*cannot* hide: there the saving does appear (−7 % wall ≈ the 8–10 µs/frame kernel gain).
|
||||
|
||||
The f32 opt2 row (4.88 s) still carried 822 k faults and is **not quotable**; it is
|
||||
also irrelevant, being the f64-pedestal class.
|
||||
|
||||
**Status: quotable** (opt3/opt4/opt5/v1 rows).
|
||||
|
||||
---
|
||||
|
||||
## 6. Campaign C — opt6 at 9×9 (where the kernel *is* the bottleneck) ★
|
||||
|
||||
**Config change**: `cluster_size=(9,9)`, `N=20 000`, `n_streams=8`,
|
||||
`max_clusters_per_frame=1500`, everything else unchanged.
|
||||
Cap 1500 chosen so the CUDA finders record **all** clusters the CPU finds.
|
||||
**Source**: same notebook, warm runs (faults ≤ 102).
|
||||
|
||||
| build | path | wall [s] | FPS | µs/frame | kernel [µs/fr] | derived ovhd | vs CPU | clusters/frame |
|
||||
|---|---|--:|--:|--:|--:|--:|--:|--:|
|
||||
| **f64 ped** | opt4 batched+pin | 1.540 | 12 988 | 77.0 | 171 † | **−94** † | 6.01× | 1422.13 |
|
||||
| **f64 ped** | opt5 graph | 1.482 | 13 491 | 74.1 | 74 * | — | 6.24× | 1422.13 |
|
||||
| **f32 ped** | opt4 batched+pin | 1.423 | 14 057 | 71.2 | 33 | +38 | 6.61× | 1422.33 |
|
||||
| **f32 ped** | opt5 graph | 1.372 | 14 580 | 68.6 | 69 * | — | 6.86× | 1422.33 |
|
||||
|
||||
† **Not quotable** — event-timer inflation under 8-stream saturation; the negative
|
||||
overhead is the tell-tale. True kernel time is 43 µs (§8). \* combined metric.
|
||||
|
||||
CPU baselines (derived from the printed speedup ratios): f64 run ≈ 9.25 s (2 162 FPS),
|
||||
f32 run ≈ 9.41 s (2 126 FPS).
|
||||
|
||||
**Result: at 9×9, opt6 is worth ~8 % end-to-end** (1.540 → 1.423 s batched;
|
||||
1.482 → 1.372 s graph) **and flips the regime**: the f64 build is kernel-bound
|
||||
(negative derived overhead), the f32 build is transfer-bound (+38 µs).
|
||||
|
||||
**Status: wall/FPS/counts quotable; kernel column must come from §8.**
|
||||
|
||||
---
|
||||
|
||||
## 7. The rule that unifies Campaigns B and C
|
||||
|
||||
Per-frame GPU operation profile at 9×9 (nsys, §8):
|
||||
|
||||
```
|
||||
f64: kernel 43 µs | D2H 19 µs | H2D 13 µs → kernel is the tallest bar
|
||||
f32: kernel 26 µs | D2H 20 µs | H2D 13 µs → near-balanced
|
||||
3×3: kernel 13–24 µs, below the ~27 µs host/PCIe floor → always hidden
|
||||
```
|
||||
|
||||
> **f32 buys exactly the distance between the kernel bar and the next-tallest bar.**
|
||||
> Nothing at 3×3 (kernel already hidden); ~8 % at 9×9/f64 (kernel on the critical
|
||||
> path); more for larger windows, faster links, or fewer streams.
|
||||
|
||||
---
|
||||
|
||||
## 8. Nsight Systems ground truth (9×9, cap 1500)
|
||||
|
||||
**Script**: [`python/tests/nsys_kernel_probe.py`](../python/tests/nsys_kernel_probe.py)
|
||||
(also in the session scratchpad). Trains 1 000 pedestal frames, runs one batched pass
|
||||
over 2 000 frames, prints wall + event `kernel_ms`.
|
||||
|
||||
```bash
|
||||
nsys profile --trace=cuda --sample=none --cpuctxsw=none -o probe_s1 \
|
||||
python nsys_kernel_probe.py 1 2000 # 1 stream → kernels serialized
|
||||
nsys profile --trace=cuda --sample=none --cpuctxsw=none -o probe_s8 \
|
||||
python nsys_kernel_probe.py 8 2000 # 8 streams → deck config
|
||||
nsys stats --report cuda_gpu_kern_sum --report cuda_gpu_mem_time_sum probe_s1.nsys-rep
|
||||
```
|
||||
|
||||
### Kernel — `aare::device::find_clusters_in_single_frame<Cluster<int,9,9,uint16_t>>`
|
||||
|
||||
| build | streams | instances | **avg** | median | min | max | σ |
|
||||
|---|--:|--:|--:|--:|--:|--:|--:|
|
||||
| f64 ped | 1 | 2 000 | **43 011 ns** | 42 945 | 40 929 | 45 632 | 580 |
|
||||
| f64 ped | 8 | 2 000 | 46 714 ns | 46 880 | 41 632 | 61 441 | 1 988 |
|
||||
| f32 ped | 1 | 2 000 | **25 602 ns** | 25 569 | 24 737 | 27 073 | 411 |
|
||||
| f32 ped | 8 | 2 000 | 25 840 ns | 25 728 | 24 864 | 30 624 | 578 |
|
||||
|
||||
> ### **opt6 kernel result: 43.0 µs → 25.6 µs = −40 %** (9×9, exclusive time)
|
||||
|
||||
### Memory operations (1-stream, uncontended)
|
||||
|
||||
| op | f64 | f32 | payload | effective BW (decimal) |
|
||||
|---|--:|--:|---|--:|
|
||||
| D2H (clusters) | 19 424 ns | 19 845 ns | 1500 × 328 B = 492 000 B (480.5 KiB) | **25.3 GB/s** (23.6 GiB/s) |
|
||||
| H2D (frame) | 13 210 ns | 13 451 ns | 400×400×2 B = 320 000 B (312.5 KiB) | **24.2 GB/s** (22.6 GiB/s) |
|
||||
| memset | 359 ns | 365 ns | — | — |
|
||||
|
||||
The H2D payload is **one frame**; the D2H payload is one frame's cluster buffer
|
||||
(`max_clusters_per_frame` × `sizeof(Cluster<int,9,9,uint16_t>)`, transferred at
|
||||
fixed size regardless of how many clusters were actually found).
|
||||
|
||||
PCIe 4.0 ×16 theoretical is 31.5 GB/s (16 GT/s × 16 lanes × 128b/130b), so H2D
|
||||
reaches **77% of theoretical** — the signature of a real DMA path. Pageable
|
||||
transfers, which the driver stages through a hidden pinned buffer, run ~15 GB/s.
|
||||
|
||||
Transfers are identical between builds and at full DMA speed — confirming pinning
|
||||
(opt4) is doing its job and that only the kernel changed.
|
||||
|
||||
### Cross-validation of the event timer
|
||||
|
||||
| build | streams | nsys avg | event `kernel_ms` | offset |
|
||||
|---|--:|--:|--:|--:|
|
||||
| f64 | 1 | 43.0 µs | 0.050 ms | +7.0 µs (launch gap) |
|
||||
| f64 | 8 | 46.7 µs | 0.149 ms | **+102 µs (queue-wait)** |
|
||||
| f32 | 1 | 25.6 µs | 0.033 ms | +7.4 µs |
|
||||
| f32 | 8 | 25.8 µs | 0.033 ms | +7.2 µs |
|
||||
|
||||
Independently reproduced unprofiled by the user (f32 build, `N=20000`): event
|
||||
`kernel_ms` = **0.030–0.032** warm (= 25.6 µs + ~5–6 µs gap); 0.047 and 0.184 on
|
||||
cold-clock first invocations.
|
||||
|
||||
Note the 8-stream instance time stretches only +9 % (f64) / +1 % (f32): **one 9×9
|
||||
kernel nearly fills the GPU, so streams queue rather than co-execute.** Multi-streaming
|
||||
at 9×9 buys transfer overlap, not kernel concurrency.
|
||||
|
||||
**Reports retained** (openable in `nsys-ui` for timeline figures):
|
||||
`probe_s1.nsys-rep`, `probe_s8.nsys-rep` (f64), `probe_f32_s1.nsys-rep`,
|
||||
`probe_f32_s8.nsys-rep` (f32).
|
||||
|
||||
---
|
||||
|
||||
## 9. Supporting study — the page-fault artifact
|
||||
|
||||
### Synthetic isolation
|
||||
|
||||
Allocating ~8 GB in ClusterVector-sized chunks (90 000 × 93 kB), touching every page,
|
||||
freeing, repeating in one process:
|
||||
|
||||
| run | wall | minor faults |
|
||||
|---|--:|--:|
|
||||
| 1 (cold heap) | 2.05 s | 2 046 594 |
|
||||
| 2 (warm heap) | 0.07 s | 2 232 |
|
||||
| 3 (warm heap) | 0.07 s | 2 016 |
|
||||
|
||||
**30× faster, 1000× fewer faults**, same allocations — glibc retains the arenas.
|
||||
|
||||
### In situ (opt2 cell, three consecutive executions)
|
||||
|
||||
| run | wall [s] | FPS | minor faults |
|
||||
|---|--:|--:|--:|
|
||||
| 1 | 6.110 | 16 366 | 2 625 948 |
|
||||
| 2 | 4.872 | 20 526 | 729 866 |
|
||||
| 3 | 4.452 | 22 459 | 185 885 |
|
||||
|
||||
Correlation Δwall vs Δfaults:
|
||||
|
||||
| interval | Δwall | Δfaults | implied cost |
|
||||
|---|--:|--:|--:|
|
||||
| 1 → 2 | 1.238 s | 1 895 667 | **0.65 µs/fault** |
|
||||
| 2 → 3 | 0.420 s | 543 981 | **0.77 µs/fault** |
|
||||
| 1 → 3 | 1.658 s | 2 439 648 | **0.68 µs/fault** |
|
||||
|
||||
Reconstruction of run 1 from run 3:
|
||||
`4.452 s + 2 439 648 × 0.68 µs = 6.111 s` vs **measured 6.110 s** (1 ms error over 6 s).
|
||||
|
||||
Kernel time was constant (0.022 ms) across all three — the GPU is not involved.
|
||||
**Conclusion: the entire first-run penalty is OS page population of the host result
|
||||
heap.** Persistence mode / clock ramp are not responsible.
|
||||
|
||||
---
|
||||
|
||||
## 10. Correctness (held constant across the whole arc)
|
||||
|
||||
### 3×3, N = 100 000
|
||||
|
||||
| finder | clusters | /frame | diff vs CPU |
|
||||
|---|--:|--:|--:|
|
||||
| CPU MT | 233 085 343 | 2330.85 | — |
|
||||
| opt1 | 233 094 770 | 2330.95 | 0.0040 % |
|
||||
| opt2 | 233 093 553 | 2330.94 | 0.0035 % |
|
||||
| opt3 / opt4 / opt5 (f64) | 233 093 484 | 2330.93 | 0.0035 % |
|
||||
| opt3 / opt4 / opt5 (**f32**) | 233 093 554 – 233 094 465 | 2330.94 | 0.0039 % |
|
||||
|
||||
The residual ~0.004 % is the known per-frame vs per-pixel pedestal-update difference
|
||||
(analysed in `ClusterFinderFrozen_vs_CUDA.ipynb`), **not** a precision effect.
|
||||
|
||||
**f64 vs f32 builds differ by ~70 clusters out of 233 M (3 × 10⁻⁷).**
|
||||
|
||||
### 9×9, N = 20 000
|
||||
|
||||
| build | clusters | /frame |
|
||||
|---|--:|--:|
|
||||
| f64 | 28 442 582 | 1422.13 |
|
||||
| f32 | 28 446 667 | 1422.33 |
|
||||
|
||||
### Two fixes that made this possible
|
||||
|
||||
**(a) B1 — per-pixel offset accumulation** (`docs/pedestal_precision_f32_cancellation.md`).
|
||||
Before: naive f32 pedestal produced **+28.06 % clusters and an unphysical high-energy
|
||||
tail** — catastrophic cancellation in `var = sum2/n − mean²` (both terms ≈ 2.17 × 10⁷,
|
||||
variance ≈ 2025) drove quiet pixels to `rms=0`, so they fired every frame. Fix: freeze
|
||||
a per-pixel baseline `X0 ≈ round(mean)` at t=0 and accumulate centered `Y = X − X0`.
|
||||
After: full-f32 matches f64 to **3 × 10⁻⁷**. **opt6 is only shippable because of B1.**
|
||||
|
||||
**(b) Test3 local-max gate backported** into
|
||||
[`clusterfinder_kernel_opt2.cuh`](../include/aare/clusterfinder_kernel_opt2.cuh)
|
||||
so opt1/opt2 stop over-counting extended charge-shared events:
|
||||
|
||||
| | before | after |
|
||||
|---|--:|--:|
|
||||
| opt1 | 233 940 268 (2339.40/fr) | 233 094 770 (2330.95/fr) |
|
||||
| opt2 | 233 931 015 (2339.31/fr) | 233 093 553 (2330.94/fr) |
|
||||
|
||||
−0.36 %, now matching CPU. Kernel time unchanged (the gate is an early `return`).
|
||||
|
||||
---
|
||||
|
||||
## 11. Numbers that must NOT be used
|
||||
|
||||
| number | why |
|
||||
|---|---|
|
||||
| 9×9 **cap = 1000** runs (f64: 1.277 s / 15 657 FPS / 7.37×; f32: 1.087 s / 18 395 FPS / 8.53×) | cluster cap saturated → **exactly** 1000.00/frame, ~30 % of clusters truncated; CPU found ~1422/frame. Timing story survives but counts and speedups are dishonest. Superseded by cap 1500 (§6). |
|
||||
| `kernel_ms` under multi-stream saturation (f64 9×9: 0.171 / 0.192; 8-stream first-launch 0.184) | event-timer queue-wait inflation up to 3.5× (§3.2). Use nsys §8. |
|
||||
| Negative "PCIe + overhead" (−0.094, −0.128, −0.007) | arithmetic artifact of `wall/N − inflated kernel`. Useful only as a *kernel-bound indicator*, never as a transfer cost. |
|
||||
| Any wall time from `nsys_kernel_probe.py` (5.3–6.6 s for 20 k frames) | fresh process → full first-touch fault tax inside the timed call; single mega-batch; profiler overhead when traced. |
|
||||
| Cold-pass notebook numbers (opt2 6.18 s, opt3 5.29 s, batched 4.45 s, graph 4.02 s) | 2.1–2.6 M unresolved page faults each (§9). |
|
||||
| `cf_cuda_v1` cell as an arc data point | builds a histogram inside the timed loop (~4.5 s / 100 k frames). Valid only as the *serialized-path* control in §5. |
|
||||
|
||||
---
|
||||
|
||||
## 12. Reproduction index
|
||||
|
||||
| artifact | path | role |
|
||||
|---|---|---|
|
||||
| main benchmark notebook | `python/tests/ClusterFinderCUDA_perf.ipynb` | all wall/FPS/count numbers; 7 instrumented timed cells |
|
||||
| nsys probe | `python/tests/nsys_kernel_probe.py` | exclusive kernel + memcpy times |
|
||||
| correctness notebook | `python/tests/ClusterFinderFrozen_vs_CUDA.ipynb` | CPU↔CUDA agreement analysis |
|
||||
| precision study | `docs/pedestal_precision_f32_cancellation.md` | B1 derivation |
|
||||
| opt1/opt2 class | `include/aare/ClusterFinderCUDAOpt2.hpp` | pre-refactor pipeline snapshot |
|
||||
| opt1/opt2 kernel | `include/aare/clusterfinder_kernel_opt2.cuh` | `88e0e8d` kernel + backported Test3 gate |
|
||||
| opt3/opt4/opt6 | `include/aare/ClusterFinderCUDA.hpp`, `include/aare/clusterfinder_kernel.cuh` | current pipeline; precision knob at lines 16–17 |
|
||||
| opt5 | `include/aare/ClusterFinderCUDA_graph.hpp` | graph-based finder |
|
||||
| deck | `docs/ClusterFinderCUDA_optimizations.pptx` | target |
|
||||
|
||||
### Protocol to reproduce a clean number
|
||||
|
||||
1. Idle machine. Optionally `sudo nvidia-smi -pm 1` (kills ~1–3 s of per-process
|
||||
driver init; does **not** affect the fault artifact).
|
||||
2. Restart the Jupyter kernel; run the notebook top to bottom (cold pass).
|
||||
3. **Re-run each timed CUDA cell** until its `minor faults` line plateaus
|
||||
(< ~200 k; usually 2–3 executions). Quote that run.
|
||||
4. The CPU baseline cannot be re-run (`ClusterFinderMT.stop()` is terminal) —
|
||||
re-run the *Build finders* + pedestal cells first if a warm CPU number is needed.
|
||||
5. For kernel times: use `nsys` at `n_streams=1`, never the event timer under load.
|
||||
|
||||
---
|
||||
|
||||
## 13. Slide-ready takeaways
|
||||
|
||||
1. **The pipeline arc (opt1 → opt5) is monotonic**: 3.14× → 4.86× → 5.58× → 7.73× →
|
||||
**8.29×** over a 48-thread CPU baseline, at identical correctness (0.004 %).
|
||||
2. **The optimization story is host overhead collapsing**: 44 → 21 → 14 → **3 µs/frame**
|
||||
against a fixed ~23 µs kernel. By opt4 the GPU is fed almost perfectly.
|
||||
3. **opt6 (f32 pedestal) is a kernel win, not always a throughput win**: kernel
|
||||
**43.0 → 25.6 µs (−40 %, nsys-verified)**, but end-to-end **0 % at 3×3** and
|
||||
**−8 % wall at 9×9**. The bottleneck is workload-dependent.
|
||||
4. **The unifying rule**: f32 buys the distance between the kernel bar and the
|
||||
next-tallest bar in the per-frame GPU profile.
|
||||
5. **opt6 is only correct because of the B1 variance rewrite** — the naive f32 pedestal
|
||||
gave +28 % clusters and an unphysical tail; the centered accumulation restores
|
||||
agreement with f64 to 3 × 10⁻⁷.
|
||||
6. **Measurement discipline was necessary to get here**: soft page faults (0.7 µs each,
|
||||
up to 4 s per run) and CUDA-event queue-wait inflation (up to 3.5×) both had to be
|
||||
identified and controlled before any number was trustworthy.
|
||||
7. **Where the next win is**: after opt4 the 3×3 pipeline is PCIe/host-bound. Further
|
||||
kernel work has no payoff at that window size; the next target is transfer volume
|
||||
(on-GPU reduction, keeping results on device) rather than kernel speed.
|
||||
@@ -0,0 +1,744 @@
|
||||
"""Rebuild docs/ClusterFinderCUDA_optimizations.pptx — opt1..opt6, in the deck's
|
||||
own design language (extracted from the original file)."""
|
||||
from pptx import Presentation
|
||||
from pptx.util import Inches as In, Pt, Emu
|
||||
from pptx.dml.color import RGBColor
|
||||
from pptx.enum.text import PP_ALIGN, MSO_ANCHOR
|
||||
from pptx.enum.shapes import MSO_SHAPE
|
||||
from lxml import etree
|
||||
from pathlib import Path
|
||||
from PIL import Image
|
||||
|
||||
FIGS = Path(__file__).parent / "figs"
|
||||
OUT = Path("/home/ferjao_k/aare/docs/ClusterFinderCUDA_optimizations.pptx")
|
||||
|
||||
# ---------------------------------------------------------------- design tokens
|
||||
BG = RGBColor(0x0B, 0x10, 0x18)
|
||||
PANEL = RGBColor(0x12, 0x1A, 0x28)
|
||||
CODEBG = RGBColor(0x0E, 0x14, 0x20)
|
||||
RULE = RGBColor(0x1E, 0x28, 0x36)
|
||||
ACCENT = RGBColor(0x1E, 0x90, 0xC2)
|
||||
AMBER = RGBColor(0xE8, 0xB2, 0x5C)
|
||||
PALE = RGBColor(0xE7, 0xED, 0xF4)
|
||||
TEXT2 = RGBColor(0xA5, 0xB2, 0xC4)
|
||||
MUTED = RGBColor(0x6B, 0x7A, 0x90)
|
||||
|
||||
UI, MONO = "Segoe UI", "Consolas"
|
||||
W, H = 13.333, 7.5
|
||||
M = 0.7 # left margin
|
||||
COL = 7.9 # left column width
|
||||
RAIL_X, RAIL_W = 9.2, 3.5 # right rail
|
||||
|
||||
prs = Presentation()
|
||||
prs.slide_width, prs.slide_height = In(W), In(H)
|
||||
BLANK = prs.slide_layouts[6]
|
||||
N_SLIDES = 19
|
||||
|
||||
|
||||
# ------------------------------------------------------------------- helpers
|
||||
def new_slide():
|
||||
s = prs.slides.add_slide(BLANK)
|
||||
bg = etree.SubElement(s._element, "{http://schemas.openxmlformats.org/presentationml/2006/main}bg")
|
||||
pr = etree.SubElement(bg, "{http://schemas.openxmlformats.org/presentationml/2006/main}bgPr")
|
||||
fill = etree.SubElement(pr, "{http://schemas.openxmlformats.org/drawingml/2006/main}solidFill")
|
||||
clr = etree.SubElement(fill, "{http://schemas.openxmlformats.org/drawingml/2006/main}srgbClr")
|
||||
clr.set("val", "0B1018")
|
||||
etree.SubElement(pr, "{http://schemas.openxmlformats.org/drawingml/2006/main}effectLst")
|
||||
s._element.insert(0, bg)
|
||||
return s
|
||||
|
||||
|
||||
def rect(s, x, y, w, h, color, shape=MSO_SHAPE.RECTANGLE):
|
||||
sh = s.shapes.add_shape(shape, In(x), In(y), In(w), In(h))
|
||||
sh.fill.solid(); sh.fill.fore_color.rgb = color
|
||||
sh.line.fill.background(); sh.shadow.inherit = False
|
||||
return sh
|
||||
|
||||
|
||||
def tb(s, x, y, w, h, anchor=MSO_ANCHOR.TOP):
|
||||
box = s.shapes.add_textbox(In(x), In(y), In(w), In(h))
|
||||
tf = box.text_frame
|
||||
tf.word_wrap = True
|
||||
tf.margin_left = tf.margin_right = tf.margin_top = tf.margin_bottom = 0
|
||||
tf.vertical_anchor = anchor
|
||||
return tf
|
||||
|
||||
|
||||
def para(tf, first=False, space_after=0, space_before=0, line=None, align=None):
|
||||
p = tf.paragraphs[0] if first else tf.add_paragraph()
|
||||
p.space_after = Pt(space_after); p.space_before = Pt(space_before)
|
||||
if line: p.line_spacing = line
|
||||
if align: p.alignment = align
|
||||
return p
|
||||
|
||||
|
||||
def run(p, text, size=11, color=TEXT2, font=UI, bold=False, italic=False, spc=None):
|
||||
r = p.add_run(); r.text = text
|
||||
f = r.font
|
||||
f.name, f.size, f.bold, f.italic = font, Pt(size), bold, italic
|
||||
f.color.rgb = color
|
||||
if spc is not None:
|
||||
r.font._rPr.set("spc", str(int(spc * 100)))
|
||||
return r
|
||||
|
||||
|
||||
# ------------------------------------------------------------------ chrome
|
||||
def chrome(s, idx, eyebrow, title, title_size=27):
|
||||
rect(s, M, 0.60, 0.35, 0.035, ACCENT)
|
||||
tf = tb(s, 1.17, 0.50, 10.33, 0.32)
|
||||
run(para(tf, True), eyebrow.upper(), 9, MUTED, bold=True, spc=1.6)
|
||||
|
||||
tf = tb(s, M, 0.86, 11.9, 1.0)
|
||||
run(para(tf, True, line=1.05), title, title_size, PALE, bold=True)
|
||||
|
||||
# progress bar
|
||||
span, n = 11.0, N_SLIDES
|
||||
pitch = span / n; wseg = pitch * 0.90
|
||||
for i in range(n):
|
||||
rect(s, M + i * pitch, 7.28, wseg, 0.045, ACCENT if i <= idx - 1 else RULE)
|
||||
tf = tb(s, 12.0, 7.14, 0.9, 0.3)
|
||||
run(para(tf, True, align=PP_ALIGN.RIGHT), f"{idx} / {n}", 8.5, MUTED)
|
||||
|
||||
|
||||
def bullets(s, x, y, w, items, size=11, gap=7):
|
||||
tf = tb(s, x, y, w, 0.3)
|
||||
for i, it in enumerate(items):
|
||||
color, txt = (it if isinstance(it, tuple) else (TEXT2, it))
|
||||
p = para(tf, i == 0, space_after=gap, line=1.25)
|
||||
run(p, "• ", size, MUTED)
|
||||
# inline emphasis with **...**
|
||||
for j, part in enumerate(txt.split("**")):
|
||||
if part:
|
||||
run(p, part, size, PALE if j % 2 else color, bold=bool(j % 2))
|
||||
return tf
|
||||
|
||||
|
||||
def code(s, x, y, w, lines, size=8.5, title=None):
|
||||
lh = 0.148
|
||||
h = 0.24 + len(lines) * lh + (0.22 if title else 0)
|
||||
rect(s, x, y, w, h, CODEBG, MSO_SHAPE.ROUNDED_RECTANGLE)
|
||||
ty = y + 0.12
|
||||
if title:
|
||||
tf = tb(s, x + 0.18, ty, w - 0.36, 0.2)
|
||||
run(para(tf, True), title, 7.5, MUTED, bold=True, spc=1.2)
|
||||
ty += 0.22
|
||||
tf = tb(s, x + 0.18, ty, w - 0.36, h - 0.24)
|
||||
for i, ln in enumerate(lines):
|
||||
p = para(tf, i == 0, line=1.12)
|
||||
if ln.strip().startswith(("//", "#")):
|
||||
run(p, ln, size, MUTED, MONO)
|
||||
continue
|
||||
for j, part in enumerate(ln.split("«")):
|
||||
for k, seg in enumerate(part.split("»")):
|
||||
if not seg: continue
|
||||
hi = (j > 0 and k == 0)
|
||||
run(p, seg, size, ACCENT if hi else TEXT2, MONO, bold=hi)
|
||||
return h
|
||||
|
||||
|
||||
def callout(s, x, y, w, text, h=0.78, color=ACCENT, size=10.5):
|
||||
rect(s, x + 0.045, y, w - 0.045, h, PANEL)
|
||||
rect(s, x, y, 0.045, h, color)
|
||||
tf = tb(s, x + 0.28, y + 0.10, w - 0.5, h - 0.2, MSO_ANCHOR.MIDDLE)
|
||||
p = para(tf, True, line=1.2)
|
||||
for j, part in enumerate(text.split("**")):
|
||||
if part:
|
||||
run(p, part, size, PALE if j % 2 else TEXT2, bold=bool(j % 2))
|
||||
|
||||
|
||||
def rail(s, items, y0=2.0, divider=True):
|
||||
if divider:
|
||||
rect(s, 8.95, 2.0, 0.012, 4.55, RULE)
|
||||
y = y0
|
||||
for it in items:
|
||||
kind = it[0]
|
||||
if kind == "label":
|
||||
tf = tb(s, RAIL_X, y, RAIL_W, 0.26)
|
||||
run(para(tf, True), it[1].upper(), 8.5, MUTED, bold=True, spc=1.4)
|
||||
y += 0.28
|
||||
elif kind == "stat":
|
||||
_, lab, val, col = it
|
||||
tf = tb(s, RAIL_X, y, RAIL_W, 0.24)
|
||||
run(para(tf, True), lab.upper(), 8.5, MUTED, spc=1.2)
|
||||
tf = tb(s, RAIL_X, y + 0.24, RAIL_W, 0.6)
|
||||
run(para(tf, True), val, 26, col, bold=True)
|
||||
y += 0.98
|
||||
elif kind == "row":
|
||||
_, lab, val, col = it
|
||||
tf = tb(s, RAIL_X, y, RAIL_W, 0.24)
|
||||
run(para(tf, True), lab.upper(), 8.5, MUTED, spc=1.2)
|
||||
tf = tb(s, RAIL_X, y + 0.22, RAIL_W, 0.3)
|
||||
run(para(tf, True), val, 13, col, bold=True)
|
||||
y += 0.66
|
||||
elif kind == "note":
|
||||
tf = tb(s, RAIL_X, y, RAIL_W, 0.9)
|
||||
run(para(tf, True, line=1.25), it[1], 9, TEXT2)
|
||||
y += 0.30 + 0.17 * (len(it[1]) // 42 + 1)
|
||||
elif kind == "gap":
|
||||
y += it[1]
|
||||
return y
|
||||
|
||||
|
||||
def figure(s, name, x, y, w):
|
||||
p = FIGS / f"{name}.png"
|
||||
iw, ih = Image.open(p).size
|
||||
h = w * ih / iw
|
||||
s.shapes.add_picture(str(p), In(x), In(y), In(w), In(h))
|
||||
return h
|
||||
|
||||
|
||||
def caption(s, x, y, w, text, size=9):
|
||||
tf = tb(s, x, y, w, 0.3)
|
||||
run(para(tf, True, line=1.25), text, size, MUTED)
|
||||
|
||||
|
||||
# =========================================================== 1 · TITLE
|
||||
s = new_slide()
|
||||
rect(s, 0, 0, 0.16, H, ACCENT)
|
||||
tf = tb(s, M + 0.3, 0.85, 11, 0.3)
|
||||
run(para(tf, True), "AARE · PSI HYBRID PIXEL DETECTORS · CUDA CLUSTERFINDER",
|
||||
9.5, MUTED, bold=True, spc=1.8)
|
||||
|
||||
tf = tb(s, M + 0.3, 1.35, 11.4, 1.7)
|
||||
run(para(tf, True, line=1.02), "Feeding the GPU", 46, PALE, bold=True)
|
||||
tf = tb(s, M + 0.3, 2.30, 11.4, 0.8)
|
||||
run(para(tf, True, line=1.05), "Six optimization steps of the CUDA ClusterFinder",
|
||||
22, ACCENT)
|
||||
|
||||
tf = tb(s, M + 0.3, 3.25, 9.6, 0.8)
|
||||
run(para(tf, True, line=1.3),
|
||||
"Five of the six steps never touch the arithmetic. They are about keeping "
|
||||
"a 24 µs kernel supplied with data — and about learning to measure honestly.",
|
||||
12.5, TEXT2)
|
||||
|
||||
stats = [("×8.29", "VS 48-THREAD CPU", ACCENT), ("39,472", "FRAMES / SECOND", PALE),
|
||||
("25.3 µs", "PER FRAME, END TO END", PALE), ("0.004%", "CLUSTER-COUNT DRIFT", AMBER)]
|
||||
for i, (v, l, c) in enumerate(stats):
|
||||
x = M + 0.3 + i * 2.85
|
||||
rect(s, x, 4.35, 0.035, 0.95, c)
|
||||
tf = tb(s, x + 0.22, 4.35, 2.5, 0.55)
|
||||
run(para(tf, True), v, 30, c, bold=True)
|
||||
tf = tb(s, x + 0.22, 4.98, 2.5, 0.3)
|
||||
run(para(tf, True), l, 8.5, MUTED, spc=1.2)
|
||||
|
||||
rect(s, M + 0.3, 6.05, 11.0, 0.012, RULE)
|
||||
tf = tb(s, M + 0.3, 6.25, 11.4, 0.6)
|
||||
run(para(tf, True, line=1.35),
|
||||
"RTX 4090 (Ada, sm_89) · PCIe 4.0 ×16 · Mönch 400×400 uint16 · 3×3 clusters · "
|
||||
"100 000 frames · Cu fluorescence, MAX IV", 10, MUTED)
|
||||
tf = tb(s, M + 0.3, 6.62, 11.4, 0.3)
|
||||
run(para(tf, True), "Khalil Ferjaoui · Paul Scherrer Institut", 10, TEXT2)
|
||||
|
||||
# =========================================================== 2 · PROBLEM
|
||||
s = new_slide()
|
||||
chrome(s, 2, "The problem & the baseline", "What has to happen to every frame")
|
||||
bullets(s, M, 1.95, COL, [
|
||||
"Per pixel: subtract a **running pedestal** (mean ± rms), keep pixels above "
|
||||
"**nσ · rms**, cut a 3×3 cluster around each local maximum.",
|
||||
"400×400 = 160 k pixels, **312.5 kB per frame**; Cu data yields ~2 330 clusters "
|
||||
"per frame at 3×3.",
|
||||
"The pedestal is **updated by every non-photon pixel**, every frame — so the "
|
||||
"arithmetic and the data movement are coupled.",
|
||||
])
|
||||
code(s, M, 3.55, COL, [
|
||||
"// the whole algorithm, per pixel",
|
||||
"v = frame[i] - pedestal_mean[i]",
|
||||
"rms = sqrt(pedestal_sum2[i]/n - pedestal_mean[i]^2)",
|
||||
"if (v > «nSigma» * rms && v == max(3x3 window)) -> emit cluster",
|
||||
"else -> update pedestal",
|
||||
], title="THE KERNEL IN FIVE LINES")
|
||||
callout(s, M, 5.55, COL,
|
||||
"**Thesis of this talk:** the compute was fast almost immediately. "
|
||||
"Five of six steps are about feeding it.")
|
||||
rail(s, [
|
||||
("label", "Baseline · same data, same threshold"),
|
||||
("gap", 0.15),
|
||||
("stat", "CPU, 1 thread", "1.75 ms", MUTED),
|
||||
("stat", "CPU MT, 48 threads", "210 µs", PALE),
|
||||
("gap", 0.1),
|
||||
("row", "That is the bar", "4 761 frames / s", TEXT2),
|
||||
("gap", 0.25),
|
||||
("note", "Every CUDA number in this deck is measured against the 48-thread "
|
||||
"CPU on the same 100 000 frames."),
|
||||
])
|
||||
|
||||
# =========================================================== 3 · THE LADDER
|
||||
s = new_slide()
|
||||
chrome(s, 3, "Roadmap", "Two acts: feed the GPU, then speed up the kernel")
|
||||
rows = [
|
||||
("opt1", "First CUDA port", "1 stream, one launch per frame", "×3.14", ACCENT),
|
||||
("opt2", "Streams + batching", "4 streams, 2 000-frame batches", "×4.86", ACCENT),
|
||||
("opt3", "Pipeline rework", "sync barriers removed", "×5.58", ACCENT),
|
||||
("opt4", "Pinned memory", "DMA-speed host transfers", "×7.73", ACCENT),
|
||||
("opt5", "CUDA Graphs", "one launch replaces six", "×8.29", ACCENT),
|
||||
("opt6", "FP32 pedestal + variance rewrite", "the first kernel change", "kernel −40%", AMBER),
|
||||
]
|
||||
y = 2.05
|
||||
for i, (tag, name, sub, gain, col) in enumerate(rows):
|
||||
rect(s, M, y, 11.9, 0.72, PANEL if i % 2 == 0 else BG)
|
||||
rect(s, M, y, 0.035, 0.72, col)
|
||||
tf = tb(s, M + 0.28, y + 0.13, 1.0, 0.4)
|
||||
run(para(tf, True), tag, 15, col, bold=True, font=MONO)
|
||||
tf = tb(s, M + 1.45, y + 0.10, 5.0, 0.3)
|
||||
run(para(tf, True), name, 13, PALE, bold=True)
|
||||
tf = tb(s, M + 1.45, y + 0.38, 5.6, 0.3)
|
||||
run(para(tf, True), sub, 10, MUTED)
|
||||
tf = tb(s, 9.4, y + 0.18, 3.1, 0.4)
|
||||
run(para(tf, True, align=PP_ALIGN.RIGHT), gain, 15, col, bold=True)
|
||||
y += 0.78
|
||||
rect(s, M, 2.05, 0.012, 4.68, RULE)
|
||||
caption(s, M, 6.72, 11.9,
|
||||
"opt1–opt5 change only how work is scheduled and moved — the arithmetic is "
|
||||
"byte-identical. opt6 is the first step that changes the kernel itself.")
|
||||
|
||||
# =========================================================== 4 · METHODOLOGY
|
||||
s = new_slide()
|
||||
chrome(s, 4, "Before any number is believed", "Three ways a GPU benchmark lies")
|
||||
items = [
|
||||
("First-touch page faults", AMBER,
|
||||
"Each run materialises ~10 GB of clusters. The first pass faults in ~2.6 M "
|
||||
"pages at 0.7 µs each — up to 4 s of pure OS work inside the timer.",
|
||||
"Fix: re-run until getrusage() minor faults plateau (< 200 k)."),
|
||||
("CUDA-event kernel timing", AMBER,
|
||||
"avg_kernel_time_ms() measures elapsed time on a stream — including waiting "
|
||||
"for other streams. Under 8-stream load it over-reads by up to 3.5×.",
|
||||
"Fix: Nsight Systems per-instance times; 1 stream for exclusive numbers."),
|
||||
("The profiler itself", AMBER,
|
||||
"Under nsys, wall time per frame inflates ~4× from API tracing.",
|
||||
"Fix: GPU op times from nsys, wall times from unprofiled runs."),
|
||||
]
|
||||
x = M
|
||||
for title, col, body, fix in items:
|
||||
rect(s, x, 2.0, 3.83, 3.15, PANEL)
|
||||
rect(s, x, 2.0, 3.83, 0.035, col)
|
||||
tf = tb(s, x + 0.26, 2.28, 3.3, 0.6)
|
||||
run(para(tf, True, line=1.15), title, 13, PALE, bold=True)
|
||||
tf = tb(s, x + 0.26, 3.02, 3.3, 1.5)
|
||||
run(para(tf, True, line=1.3), body, 10, TEXT2)
|
||||
tf = tb(s, x + 0.26, 4.42, 3.3, 0.65)
|
||||
run(para(tf, True, line=1.25), fix, 9.5, ACCENT)
|
||||
x += 4.03
|
||||
code(s, M, 5.4, 11.9, [
|
||||
"# every timed cell in the benchmark notebook is bracketed with:",
|
||||
"mf0 = resource.getrusage(resource.RUSAGE_SELF).ru_minflt",
|
||||
"... t = time.perf_counter() - t0 ...",
|
||||
"print(f'minor faults: {mf1-mf0:,}') # quote the run where this plateaus",
|
||||
], title="THE FAULT PROTOCOL — python/tests/ClusterFinderCUDA_perf.ipynb")
|
||||
callout(s, M, 6.52, 11.2,
|
||||
"Validated: **wall = steady-state + faults × 0.68 µs** reproduced a 6.110 s "
|
||||
"run to within **1 ms**. Kernel time stayed constant throughout — the GPU was never the variable.",
|
||||
h=0.70, size=9.5)
|
||||
|
||||
# =========================================================== 5 · OPT1
|
||||
s = new_slide()
|
||||
chrome(s, 5, "opt1 · the first CUDA port", "One frame, one stream, fully synchronous")
|
||||
bullets(s, M, 1.95, COL, [
|
||||
"Shared-memory tiling with **halo loading** for any cluster size; pedestal "
|
||||
"subtraction fused into the tile load.",
|
||||
"Cluster geometry is a **compile-time template parameter** → the 3×3 stencil "
|
||||
"is fully unrolled.",
|
||||
"One cudaMemcpy in, one kernel, one cudaMemcpy out — **the host blocks "
|
||||
"on every frame**.",
|
||||
])
|
||||
code(s, M, 3.62, COL, [
|
||||
"// one frame at a time — the host waits at every step",
|
||||
"cudaMemcpy(d_frame, h_frame, bytes, cudaMemcpyHostToDevice);",
|
||||
"find_clusters_in_single_frame<ClusterType, FRAME_TYPE>",
|
||||
" <<<grid, block, shmem>>>(d_frame, d_pd_mean, ...);",
|
||||
"cudaMemcpy(h_out, d_out, out_bytes, cudaMemcpyDeviceToHost);",
|
||||
], title="ClusterFinderCUDAOpt2.hpp · find_clusters()")
|
||||
callout(s, M, 5.62, COL,
|
||||
"The stencil was **already fast**: 23 µs of kernel inside a 67 µs frame. "
|
||||
"The other 44 µs is the host standing still.")
|
||||
rail(s, [
|
||||
("label", "opt1 · 3×3 · 100 k frames"),
|
||||
("gap", 0.15),
|
||||
("stat", "End to end", "66.8 µs", PALE),
|
||||
("stat", "vs 48-thread CPU", "×3.14", ACCENT),
|
||||
("gap", 0.05),
|
||||
("row", "Kernel (GPU)", "23 µs", TEXT2),
|
||||
("row", "Host + PCIe", "44 µs", AMBER),
|
||||
("gap", 0.2),
|
||||
("note", "Two thirds of the frame is spent not computing."),
|
||||
])
|
||||
|
||||
# =========================================================== 6 · OPT2
|
||||
s = new_slide()
|
||||
chrome(s, 6, "opt2 · streams and batching", "A CUDA stream is a queue the GPU can overlap")
|
||||
bullets(s, M, 1.95, COL, [
|
||||
"A **stream** is an ordered queue of GPU work. Work in **different** streams may "
|
||||
"overlap — so a copy can run while another stream computes.",
|
||||
"Each stream gets its own **StreamContext**: device frame buffer, output buffer "
|
||||
"and pedestal. Frames are handed out **round-robin**.",
|
||||
"The host now submits **2 000 frames per call** instead of one.",
|
||||
])
|
||||
code(s, M, 3.62, COL, [
|
||||
"struct StreamContext {",
|
||||
" cudaStream_t stream;",
|
||||
" FRAME_TYPE *d_frame; ClusterType *d_clusters;",
|
||||
" PEDESTAL_TYPE *d_pd_mean, *d_pd_sum, *d_pd_sum2;",
|
||||
"};",
|
||||
"auto &sc = v_sc[frame_idx % «n_streams»]; // round-robin",
|
||||
], title="ClusterFinderCUDA.hpp · per-stream state")
|
||||
callout(s, M, 5.72, COL,
|
||||
"**Scaffolding, not yet the payoff.** The streams exist, but the host still "
|
||||
"synchronises after every round — see opt3.")
|
||||
rail(s, [
|
||||
("label", "opt2 · 4 streams · batch 2 000"),
|
||||
("gap", 0.15),
|
||||
("stat", "End to end", "43.2 µs", PALE),
|
||||
("stat", "vs CPU", "×4.86", ACCENT),
|
||||
("gap", 0.05),
|
||||
("row", "Step gain over opt1", "×1.55", ACCENT),
|
||||
("row", "Host + PCIe", "21 µs (was 44)", AMBER),
|
||||
])
|
||||
|
||||
# =========================================================== 7 · OPT3
|
||||
s = new_slide()
|
||||
chrome(s, 7, "opt3 · remove the sync barriers", "Stop draining the GPU between rounds")
|
||||
bullets(s, M, 1.95, 7.4, [
|
||||
"opt2 synchronised **all streams after every round** of n_streams frames. "
|
||||
"The GPU drained to empty each time.",
|
||||
"opt3 submits every frame's H2D → kernel → D2H **asynchronously**, then "
|
||||
"synchronises **once at the end of the batch**.",
|
||||
], size=10.5)
|
||||
figure(s, "fig_streams", M, 3.05, 7.55)
|
||||
code(s, 8.35, 1.95, 4.25, [
|
||||
"// opt2: barrier after every round",
|
||||
"for (round) {",
|
||||
" submit(n_streams frames);",
|
||||
" «cudaDeviceSynchronize»();",
|
||||
"}",
|
||||
"",
|
||||
"// opt3: submit everything, sync once",
|
||||
"for (frame : batch) {",
|
||||
" cudaMemcpyAsync(..., sc.stream);",
|
||||
" kernel<<<..., sc.stream>>>(...);",
|
||||
" cudaMemcpyAsync(..., sc.stream);",
|
||||
"}",
|
||||
"for (sc : streams)",
|
||||
" «cudaStreamSynchronize»(sc.stream);",
|
||||
], size=8, title="THE ONE-LINE IDEA")
|
||||
callout(s, 8.35, 5.05, 4.25,
|
||||
"**37.6 µs/frame · ×5.58**\nHost overhead 21 → **14 µs**", h=0.86, size=11)
|
||||
caption(s, 8.35, 6.15, 4.25,
|
||||
"Each lane is one stream. Removing the barrier lets a stream start its next "
|
||||
"frame while its neighbours are still copying.")
|
||||
|
||||
# =========================================================== 8 · OPT4
|
||||
s = new_slide()
|
||||
chrome(s, 8, "opt4 · pinned (page-locked) memory", "What pinning is, and why the GPU cares")
|
||||
bullets(s, M, 1.95, 12.0, [
|
||||
"Normal host memory is **pageable** — the OS may move or swap it. A DMA engine "
|
||||
"cannot safely read that, so the driver first copies your data into a **hidden "
|
||||
"pinned staging buffer**. Every transfer is copied twice.",
|
||||
"**Pinning** locks the pages in physical RAM. The GPU's DMA engine then reads "
|
||||
"host memory **directly** — no staging copy, and the transfer can be truly asynchronous.",
|
||||
], size=10.5)
|
||||
figure(s, "fig_pinning", M, 3.15, 7.6)
|
||||
code(s, 8.5, 3.15, 4.1, [
|
||||
"// pin the whole dataset once",
|
||||
"«cudaHostRegister»(ptr, bytes,",
|
||||
" cudaHostRegisterDefault);",
|
||||
"",
|
||||
"// ... run the whole campaign ...",
|
||||
"",
|
||||
"«cudaHostUnregister»(ptr);",
|
||||
], size=8, title="ClusterFinderCUDA.hpp")
|
||||
callout(s, 8.5, 4.72, 4.1,
|
||||
"**27.2 µs/frame · ×7.73**\nHost overhead 14 → **3 µs**", h=0.86, size=11)
|
||||
caption(s, 8.5, 5.82, 4.1,
|
||||
"Measured H2D: one 400×400 uint16 frame (312.5 KiB = 320 000 B) in 13.2 µs "
|
||||
"= 24.2 GB/s — 77% of PCIe 4.0 ×16 theoretical, i.e. true DMA speed. "
|
||||
"Pageable staging runs ~15 GB/s.")
|
||||
caption(s, M, 6.62, 7.6,
|
||||
"Caveat: pinned memory is a finite system resource — it cannot be swapped. "
|
||||
"aare exposes a budget helper (print_pinning_budget) before you pin 31 GB.")
|
||||
|
||||
# =========================================================== 9 · OPT5
|
||||
s = new_slide()
|
||||
chrome(s, 9, "opt5 · CUDA Graphs", "Record the pipeline once, replay it with one call")
|
||||
bullets(s, M, 1.95, 12.0, [
|
||||
"Every cudaMemcpyAsync / kernel launch costs the **CPU** a few microseconds of "
|
||||
"driver work — per frame, per operation. At 39 k frames/s that is the budget.",
|
||||
"A **CUDA Graph** captures the whole dependency DAG once. Replaying it is a "
|
||||
"**single** cudaGraphLaunch — the driver already knows every node and edge.",
|
||||
], size=10.5)
|
||||
figure(s, "fig_graphs", M, 3.15, 7.6)
|
||||
code(s, 8.5, 3.15, 4.1, [
|
||||
"// record once, at setup",
|
||||
"cudaStreamBeginCapture(sc.stream, ...);",
|
||||
" submit_h2d_kernel_d2h(sc);",
|
||||
"cudaStreamEndCapture(sc.stream, &sc.graph);",
|
||||
"«cudaGraphInstantiate»(&sc.graphExec, ...);",
|
||||
"",
|
||||
"// per batch — one call",
|
||||
"«cudaGraphLaunch»(sc.graphExec, sc.stream);",
|
||||
], size=8, title="ClusterFinderCUDA_graph.hpp")
|
||||
callout(s, 8.5, 4.88, 4.1,
|
||||
"**25.3 µs/frame · ×8.29**\nBest end-to-end result", h=0.86, size=11)
|
||||
caption(s, M, 6.62, 12.0,
|
||||
"Worth 7% here because opt4 already removed the transfer cost — what remains "
|
||||
"is CPU launch overhead, which is exactly what graphs eliminate. "
|
||||
"Trade-off: shapes are frozen at record time, so the batch geometry must be fixed.")
|
||||
|
||||
# =========================================================== 10 · OPT6 why
|
||||
s = new_slide()
|
||||
chrome(s, 10, "opt6 · FP32 device pedestal", "The first change to the kernel itself")
|
||||
bullets(s, M, 1.95, 7.5, [
|
||||
"~99.9% of pixels take the **pedestal-update** branch, which reads and writes "
|
||||
"mean, sum and sum². In FP64 that is **32 bytes per pixel**; in FP32, 16.",
|
||||
"The kernel is **bandwidth-bound**, so halving that traffic nearly halves the time.",
|
||||
"Second effect: on a GeForce part, **FP64 arithmetic runs at 1/64 of FP32**. "
|
||||
"The pedestal update was paying that tax on every pixel.",
|
||||
], size=10.5)
|
||||
figure(s, "fig_f32_kernel", M, 4.05, 7.5)
|
||||
code(s, 8.5, 1.95, 4.1, [
|
||||
"// clusterfinder_kernel.cuh",
|
||||
"using COMPUTE_TYPE = float;",
|
||||
"using DEVICE_PED_TYPE = «float»;",
|
||||
"// was: double",
|
||||
], size=8.5, title="ONE TYPEDEF")
|
||||
callout(s, 8.5, 3.05, 4.1,
|
||||
"Kernel, 9×9, measured with Nsight Systems\n**43.0 µs → 25.6 µs (−40%)**",
|
||||
h=0.92, size=11)
|
||||
caption(s, 8.5, 4.20, 4.1,
|
||||
"Exclusive per-instance kernel time, 2 000 instances, 1 stream. "
|
||||
"σ = 0.4 µs. Transfers are unchanged, as they must be.")
|
||||
callout(s, 8.5, 5.30, 4.1,
|
||||
"But a faster kernel is **not automatically a faster frame** — and naive FP32 "
|
||||
"is **wrong**. Both on the next slides.", h=1.05, size=10, color=AMBER)
|
||||
|
||||
# =========================================================== 11 · OPT6 trap
|
||||
s = new_slide()
|
||||
chrome(s, 11, "opt6 · the correctness trap", "Why the obvious FP32 pedestal is broken")
|
||||
bullets(s, M, 1.95, 7.9, [
|
||||
"The running variance was computed as **var = E[X²] − mean²**. With a pedestal "
|
||||
"mean of ~4 655 ADU, both terms are ≈ 2.17 × 10⁷ while the answer is ≈ 2 000.",
|
||||
"In FP32 the spacing between representable numbers at 2.17 × 10⁷ is **2 048** — "
|
||||
"larger than the variance itself. This is **catastrophic cancellation**.",
|
||||
], size=10.5)
|
||||
figure(s, "fig_cancellation", M, 3.25, 7.5)
|
||||
rail(s, [
|
||||
("label", "What it looked like"),
|
||||
("gap", 0.15),
|
||||
("stat", "Extra clusters", "+28.06%", AMBER),
|
||||
("gap", 0.05),
|
||||
("note", "Quiet pixels got rms → 0, so their threshold became 0 and they fired "
|
||||
"on every single frame — producing a large unphysical high-energy tail "
|
||||
"in the spectrum."),
|
||||
("gap", 0.35),
|
||||
("row", "Affected pixels", "~1–2% of the sensor", TEXT2),
|
||||
("row", "Written up in", "docs/pedestal_precision_…", MUTED),
|
||||
])
|
||||
|
||||
# =========================================================== 12 · OPT6 fix
|
||||
s = new_slide()
|
||||
chrome(s, 12, "opt6 · the variance rewrite", "Accumulate what is small, not what is large")
|
||||
bullets(s, M, 1.95, COL, [
|
||||
"Freeze a per-pixel baseline **X₀ = round(mean)** once, at the end of pedestal "
|
||||
"training, and never move it again.",
|
||||
"Accumulate the **centred** value Y = X − X₀ instead of X. Now both sums are "
|
||||
"O(rms)-sized — the huge common term is gone before the subtraction.",
|
||||
"The reported pedestal mean is still the full value, **X₀ + sum/n**, so nothing "
|
||||
"downstream changes.",
|
||||
])
|
||||
code(s, M, 4.0, COL, [
|
||||
"// before — both terms ~2.17e7, answer ~2000",
|
||||
"var = sum2/n - mean*mean;",
|
||||
"",
|
||||
"// after — centred on a frozen per-pixel offset X0",
|
||||
"DEVICE_PED_TYPE resid = mean - «d_pd_off»[i]; // ~O(1)",
|
||||
"DEVICE_PED_TYPE var_px = sum2[i]/n - resid*resid; // no cancellation",
|
||||
], title="clusterfinder_kernel.cuh")
|
||||
callout(s, M, 6.05, COL,
|
||||
"Result: the 100% FP32 build now matches the FP64 build to "
|
||||
"**3 × 10⁻⁷** — 70 clusters out of 233 million.")
|
||||
rail(s, [
|
||||
("label", "Why it works"),
|
||||
("gap", 0.15),
|
||||
("note", "Precision is relative. Floats resolve small numbers finely and large "
|
||||
"numbers coarsely — so never let a small answer be the difference of "
|
||||
"two large numbers."),
|
||||
("gap", 0.5),
|
||||
("row", "f32 vs f64 counts", "3 × 10⁻⁷", ACCENT),
|
||||
("row", "vs CPU", "0.0039%", ACCENT),
|
||||
("gap", 0.3),
|
||||
("note", "X₀ must never be updated — the accumulators are defined relative to it."),
|
||||
])
|
||||
|
||||
# =========================================================== 13 · OPT6 when
|
||||
s = new_slide()
|
||||
chrome(s, 13, "opt6 · when does a faster kernel help?", "Only if the kernel was the tallest bar")
|
||||
figure(s, "fig_bottleneck", M, 2.05, 11.9)
|
||||
callout(s, M, 5.35, 11.9,
|
||||
"**f32 buys exactly the distance between the kernel and the next-tallest bar.** "
|
||||
"At 3×3 the kernel already hides inside the transfers, so end-to-end throughput does not move at all. "
|
||||
"At 9×9 the kernel is on the critical path, and the same change is worth 8% of the frame.",
|
||||
h=0.95, size=11)
|
||||
caption(s, M, 6.55, 11.9,
|
||||
"9×9, 20 000 frames, 8 streams, cap 1 500 (all CPU clusters recorded): "
|
||||
"1.540 s → 1.423 s batched, 1.482 s → 1.372 s with graphs. "
|
||||
"The FP64 build shows the classic kernel-bound signature — per-frame wall time shorter than per-frame kernel time, because kernels from different streams queue.")
|
||||
|
||||
# =========================================================== 14 · RESULTS
|
||||
s = new_slide()
|
||||
chrome(s, 14, "Results", "The whole ladder, one dataset, one baseline")
|
||||
figure(s, "fig_arc", M, 1.95, 11.9)
|
||||
callout(s, M, 5.55, 5.85,
|
||||
"**×8.29 over 48 CPU threads** — 21.0 s → 2.53 s for 100 000 frames.", h=0.8)
|
||||
callout(s, 6.75, 5.55, 5.85,
|
||||
"Every step is **monotonic**, and correctness is held constant at **0.004%** throughout.",
|
||||
h=0.8, color=AMBER)
|
||||
caption(s, M, 6.55, 11.9,
|
||||
"3×3 clusters · nσ = 5 · 100 000 frames · batch 2 000 · 4 streams · warm run "
|
||||
"(page faults plateaued) · CPU baseline = ClusterFinderMT with 48 threads.")
|
||||
|
||||
# =========================================================== 15 · WHERE TIME GOES
|
||||
s = new_slide()
|
||||
chrome(s, 15, "Where the time actually went", "The kernel never changed — the overhead collapsed")
|
||||
figure(s, "fig_overhead", M, 2.05, 5.9)
|
||||
bullets(s, 7.0, 2.15, 5.6, [
|
||||
"The GPU kernel is a **flat ~23 µs** across opt1–opt5. Not one of those steps "
|
||||
"made the arithmetic faster.",
|
||||
"What changed is everything around it: **44 → 21 → 14 → 3 µs** of host and PCIe "
|
||||
"time per frame.",
|
||||
"By opt4 the pipeline is **fed almost perfectly** — which is precisely why opt5 "
|
||||
"(launch overhead) is the only lever left, and why opt6 shows nothing at 3×3.",
|
||||
], size=10.5)
|
||||
callout(s, 7.0, 5.15, 5.6,
|
||||
"The bottleneck moved from **the host**, to **PCIe**, and finally — only for "
|
||||
"large cluster windows — to **the kernel**.", h=0.95, size=10.5)
|
||||
caption(s, M, 6.4, 11.9,
|
||||
"Per-frame GPU operation profile at 9×9 (Nsight Systems): kernel 43 µs · "
|
||||
"D2H 19.4 µs · H2D 13.2 µs. At 3×3 the kernel is 13–24 µs against a ~25 µs "
|
||||
"transfer-and-host floor — which is the entire story of this deck in three numbers.")
|
||||
|
||||
# =========================================================== 16 · CORRECTNESS
|
||||
s = new_slide()
|
||||
chrome(s, 16, "Validation", "Same physics out of every variant")
|
||||
figure(s, "fig_correctness", M, 2.0, 7.6)
|
||||
bullets(s, 8.6, 2.05, 4.1, [
|
||||
"233 million clusters over 100 000 frames.",
|
||||
"All CUDA variants agree with the CPU to **0.004%**.",
|
||||
"The residual is **not** precision: it is the CUDA finder updating the pedestal "
|
||||
"**once per frame** vs the CPU's per-pixel update.",
|
||||
], size=10)
|
||||
code(s, M, 4.85, 7.6, [
|
||||
"CPU (ClusterFinderMT) 233 085 343 2330.85 / frame reference",
|
||||
"opt1 .. opt5 (f64 ped) 233 093 484 2330.93 / frame 0.0035 %",
|
||||
"opt6 (f32 ped) 233 094 465 2330.94 / frame 0.0039 %",
|
||||
], size=8.5, title="CLUSTER COUNTS · 3×3 · 100 000 FRAMES")
|
||||
callout(s, 8.6, 4.85, 4.1,
|
||||
"A local-maximum gate had to be **back-ported** into the opt1/opt2 snapshot so "
|
||||
"the whole ladder is compared at identical correctness.", h=1.0, size=10, color=AMBER)
|
||||
caption(s, M, 6.35, 11.9,
|
||||
"Cross-checks: ClusterFinderFrozen (a CPU finder with the CUDA pedestal-update "
|
||||
"timing) isolates that residual; energy spectra overlay within statistics.")
|
||||
|
||||
# =========================================================== 17 · API 1
|
||||
s = new_slide()
|
||||
chrome(s, 17, "For users · Python API", "The fast path in eight lines")
|
||||
code(s, M, 1.95, 7.6, [
|
||||
"from aare import File, ClusterFinderCUDA",
|
||||
"",
|
||||
"cf = ClusterFinderCUDA(image_size=(400, 400), cluster_size=(3, 3),",
|
||||
" n_sigma=5, «n_streams»=4,",
|
||||
" «max_clusters_per_frame»=3000)",
|
||||
"",
|
||||
"for _ in range(1000): # 1. train the pedestal",
|
||||
" cf.push_pedestal_frame(pd.read_frame())",
|
||||
"",
|
||||
"data = f.read_n(100_000) # 2. one contiguous array",
|
||||
"cf.«register_input_buffer»(data) # 3. pin it once",
|
||||
"",
|
||||
"for s in range(0, N, 2000): # 4. batch through it",
|
||||
" clusters = cf.«find_clusters_batched»(data[s:s+2000], first_frame=s)",
|
||||
"",
|
||||
"cf.unregister_input_buffer() # 5. release the pages",
|
||||
], size=9, title="THE RECOMMENDED PATTERN")
|
||||
bullets(s, 8.6, 2.0, 4.1, [
|
||||
"find_clusters_batched returns **one ClusterVector per frame**, in order.",
|
||||
"register_input_buffer is what turns opt3 into opt4 — **one call**.",
|
||||
"Swap in ClusterFinderCUDAGraph for opt5; the API is identical.",
|
||||
], size=10)
|
||||
callout(s, 8.6, 4.55, 4.1,
|
||||
"GIL is released for both find_clusters and find_clusters_batched, so "
|
||||
"reading the next file can overlap with the GPU.", h=1.05, size=10)
|
||||
callout(s, 8.6, 5.80, 4.1,
|
||||
"Pin **once**, outside the loop. Slices of a registered array inherit the pinning.",
|
||||
h=0.85, size=10, color=AMBER)
|
||||
|
||||
# =========================================================== 18 · API 2
|
||||
s = new_slide()
|
||||
chrome(s, 18, "For users · choosing the knobs", "What to set, and what it costs you")
|
||||
hdr = [("Parameter", 1.05), ("What it does", 3.6), ("Guidance", 5.2)]
|
||||
y = 2.0
|
||||
rect(s, M, y, 11.9, 0.4, PANEL)
|
||||
for lab, dx in hdr:
|
||||
tf = tb(s, M + dx - 0.85 if dx > 1.05 else M + 0.28, y + 0.09, 5.0, 0.3)
|
||||
run(para(tf, True), lab.upper(), 9, MUTED, bold=True, spc=1.3)
|
||||
y += 0.44
|
||||
params = [
|
||||
("n_streams", "How many frames are in flight at once.",
|
||||
"4 is right for 3×3. Larger windows saturate the GPU — 8 helps at 9×9."),
|
||||
("max_clusters_per_frame", "Fixed size of the per-frame D2H transfer.",
|
||||
"Must exceed the real maximum or clusters are silently dropped. Too high wastes PCIe."),
|
||||
("batch size", "Frames per find_clusters_batched call.",
|
||||
"2 000 amortises launch overhead without a large pinned footprint."),
|
||||
("cluster_size", "Compile-time stencil geometry.",
|
||||
"3×3 and 9×9 are registered; 9×9 shifts the bottleneck onto the kernel."),
|
||||
("register_input_buffer", "Page-locks the host array for DMA.",
|
||||
"Always, if the data is already in RAM. Check the pinning budget first."),
|
||||
]
|
||||
for i, (p_, what, guide) in enumerate(params):
|
||||
if i % 2 == 0:
|
||||
rect(s, M, y, 11.9, 0.82, PANEL)
|
||||
tf = tb(s, M + 0.28, y + 0.14, 2.6, 0.5)
|
||||
run(para(tf, True, line=1.1), p_, 9.5, ACCENT, font=MONO, bold=True)
|
||||
tf = tb(s, M + 3.0, y + 0.14, 2.9, 0.6)
|
||||
run(para(tf, True, line=1.2), what, 9.5, PALE)
|
||||
tf = tb(s, M + 6.15, y + 0.14, 5.4, 0.6)
|
||||
run(para(tf, True, line=1.2), guide, 9.5, TEXT2)
|
||||
y += 0.80
|
||||
callout(s, M, 6.55, 11.2,
|
||||
"The single most common mistake: leaving **max_clusters_per_frame** too low. "
|
||||
"It does not error — it truncates, and every frame quietly returns the same count.",
|
||||
h=0.66, size=10, color=AMBER)
|
||||
|
||||
# =========================================================== 19 · NEXT
|
||||
s = new_slide()
|
||||
chrome(s, 19, "Where this leaves us", "The bottleneck has moved — twice")
|
||||
cards = [
|
||||
("DONE", ACCENT, "×8.29 over 48 CPU threads",
|
||||
"25.3 µs/frame end to end, at 0.004% cluster agreement. Five pipeline steps "
|
||||
"and one kernel step."),
|
||||
("DONE", ACCENT, "FP32 pedestal, safely",
|
||||
"−40% kernel, and correct — because the variance is now accumulated on a frozen "
|
||||
"per-pixel offset instead of a raw second moment."),
|
||||
("NEXT", AMBER, "Attack the transfers, not the kernel",
|
||||
"At 3×3 the kernel is already hidden. The remaining per-frame cost is PCIe: "
|
||||
"492 kB of clusters out, 312 kB of frame in."),
|
||||
("NEXT", AMBER, "Keep results on the device",
|
||||
"On-GPU reduction, eta/interpolation on device, or compressed cluster formats — "
|
||||
"so the D2H bar stops setting the floor."),
|
||||
]
|
||||
x, y = M, 2.05
|
||||
for i, (tag, col, title, body) in enumerate(cards):
|
||||
cx = M + (i % 2) * 6.05
|
||||
cy = 2.05 + (i // 2) * 2.35
|
||||
rect(s, cx, cy, 5.85, 2.05, PANEL)
|
||||
rect(s, cx, cy, 5.85, 0.035, col)
|
||||
tf = tb(s, cx + 0.3, cy + 0.26, 1.4, 0.26)
|
||||
run(para(tf, True), tag, 8.5, col, bold=True, spc=1.5)
|
||||
tf = tb(s, cx + 0.3, cy + 0.60, 5.2, 0.4)
|
||||
run(para(tf, True, line=1.1), title, 14, PALE, bold=True)
|
||||
tf = tb(s, cx + 0.3, cy + 1.12, 5.2, 0.85)
|
||||
run(para(tf, True, line=1.3), body, 10, TEXT2)
|
||||
callout(s, M, 6.58, 11.2,
|
||||
"Full numbers, methodology and reproduction steps: **docs/benchmark_opt1_opt6_results.md** · "
|
||||
"notebook **python/tests/ClusterFinderCUDA_perf.ipynb** · profiler probe **python/tests/nsys_kernel_probe.py**",
|
||||
h=0.66, size=9.5)
|
||||
|
||||
prs.save(OUT)
|
||||
print(f"saved {OUT} ({len(prs.slides.__iter__.__self__._sldIdLst)} slides)")
|
||||
@@ -0,0 +1,363 @@
|
||||
"""Figures for ClusterFinderCUDA_optimizations.pptx — deck palette, dark, transparent."""
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from matplotlib.patches import FancyArrowPatch, Rectangle
|
||||
from pathlib import Path
|
||||
|
||||
OUT = Path(__file__).parent / "figs"
|
||||
OUT.mkdir(exist_ok=True)
|
||||
|
||||
BG = "#0B1018"
|
||||
PANEL = "#121A28"
|
||||
RULE = "#1E2836"
|
||||
ACCENT = "#1E90C2" # data 1
|
||||
AMBER = "#E8B25C" # data 2
|
||||
PALE = "#E7EDF4" # data 3 / primary text
|
||||
TEXT2 = "#A5B2C4"
|
||||
MUTED = "#6B7A90" # non-data only: grid, axes, annotation
|
||||
|
||||
plt.rcParams.update({
|
||||
"font.family": "DejaVu Sans", "font.size": 9,
|
||||
"text.color": PALE, "axes.labelcolor": TEXT2,
|
||||
"xtick.color": TEXT2, "ytick.color": TEXT2,
|
||||
"axes.edgecolor": RULE, "axes.facecolor": "none",
|
||||
"figure.facecolor": BG, "savefig.facecolor": BG,
|
||||
"axes.grid": False, "svg.fonttype": "none",
|
||||
})
|
||||
|
||||
|
||||
def save(fig, name):
|
||||
fig.savefig(OUT / f"{name}.png", dpi=220, transparent=False,
|
||||
bbox_inches="tight", pad_inches=0.08)
|
||||
plt.close(fig)
|
||||
print("wrote", name)
|
||||
|
||||
|
||||
def bare(ax, keep=("left", "bottom")):
|
||||
for s in ("top", "right", "left", "bottom"):
|
||||
ax.spines[s].set_visible(s in keep)
|
||||
|
||||
|
||||
# ---------------------------------------------------------------- 1. the arc
|
||||
def fig_arc():
|
||||
steps = ["CPU MT\n48 threads", "opt1\n1 stream", "opt2\nstreams+batch",
|
||||
"opt3\npipeline", "opt4\npinned", "opt5\ngraphs"]
|
||||
fps = [4761, 14968, 23134, 26588, 36810, 39472]
|
||||
spd = [1.0, 3.14, 4.86, 5.58, 7.73, 8.29]
|
||||
colors = [MUTED] + [ACCENT] * 4 + [AMBER]
|
||||
|
||||
fig, ax = plt.subplots(figsize=(11.4, 3.5))
|
||||
x = np.arange(len(steps))
|
||||
bars = ax.bar(x, fps, width=0.62, color=colors, zorder=3)
|
||||
for b in bars:
|
||||
b.set_linewidth(0)
|
||||
for xi, (f, s) in enumerate(zip(fps, spd)):
|
||||
ax.text(xi, f + 900, f"{f:,}", ha="center", va="bottom",
|
||||
color=PALE, fontsize=11, fontweight="bold")
|
||||
ax.text(xi, f + 3100, ("baseline" if s == 1.0 else f"×{s:.2f}"),
|
||||
ha="center", va="bottom", color=AMBER if s > 8 else TEXT2, fontsize=9)
|
||||
ax.set_xticks(x)
|
||||
ax.set_xticklabels(steps, fontsize=9, color=TEXT2)
|
||||
ax.set_ylim(0, 47000)
|
||||
ax.set_yticks([])
|
||||
bare(ax, keep=("bottom",))
|
||||
ax.spines["bottom"].set_color(RULE)
|
||||
ax.set_ylabel("")
|
||||
ax.text(0, 45500, "frames / second · 3×3 clusters, 100 000 frames, warm run",
|
||||
color=MUTED, fontsize=9, ha="left")
|
||||
save(fig, "fig_arc")
|
||||
|
||||
|
||||
# ------------------------------------------------- 2. host overhead collapse
|
||||
def fig_overhead():
|
||||
steps = ["opt1", "opt2", "opt3", "opt4"]
|
||||
ovhd = [44, 21, 14, 3]
|
||||
kern = [23, 22, 24, 24]
|
||||
|
||||
fig, ax = plt.subplots(figsize=(5.6, 3.0))
|
||||
x = np.arange(len(steps))
|
||||
ax.bar(x, kern, width=0.55, color=ACCENT, zorder=3, label="kernel (GPU)")
|
||||
ax.bar(x, ovhd, width=0.55, bottom=kern, color=AMBER, zorder=3,
|
||||
label="host + PCIe overhead")
|
||||
for xi, (k, o) in enumerate(zip(kern, ovhd)):
|
||||
ax.text(xi, k + o + 1.6, f"{o} µs", ha="center", color=AMBER,
|
||||
fontsize=10, fontweight="bold")
|
||||
ax.set_xticks(x); ax.set_xticklabels(steps, color=TEXT2)
|
||||
ax.set_ylabel("µs / frame", color=TEXT2)
|
||||
ax.set_ylim(0, 78)
|
||||
bare(ax)
|
||||
ax.legend(frameon=False, fontsize=8.5, labelcolor=TEXT2, loc="upper right")
|
||||
save(fig, "fig_overhead")
|
||||
|
||||
|
||||
# ------------------------------------------------------ 3. streams timeline
|
||||
def fig_streams():
|
||||
fig, axes = plt.subplots(3, 1, figsize=(7.7, 3.9))
|
||||
H, K, D = 12, 22, 12
|
||||
FR = H + K + D
|
||||
LANE = 0.68
|
||||
|
||||
def frame(ax, lane_y, t0):
|
||||
ax.broken_barh([(t0, H)], (lane_y, LANE), facecolors=AMBER, zorder=3)
|
||||
ax.broken_barh([(t0 + H, K)], (lane_y, LANE), facecolors=ACCENT, zorder=3)
|
||||
ax.broken_barh([(t0 + H + K, D)], (lane_y, LANE), facecolors=PALE, zorder=3)
|
||||
|
||||
# --- opt1: one stream, strictly serial
|
||||
ax = axes[0]
|
||||
for i in range(3):
|
||||
frame(ax, 1.0, i * FR)
|
||||
ax.set_ylim(0.4, 2.3)
|
||||
ax.text(3 * FR + 6, 1.34, "GPU idle between every stage", color=MUTED, fontsize=7.5,
|
||||
va="center")
|
||||
|
||||
# --- opt2: 4 streams, barrier after each round
|
||||
ax = axes[1]
|
||||
ROUND = FR + 3 * 8
|
||||
for r in range(2):
|
||||
for st in range(4):
|
||||
frame(ax, 3 - st * 1.0, r * (ROUND + 26) + st * 8)
|
||||
ax.axvspan(ROUND, ROUND + 26, color=AMBER, alpha=0.13, zorder=1)
|
||||
ax.text(ROUND + 13, 4.15, "barrier — GPU drains", color=AMBER, fontsize=7.5,
|
||||
ha="center", va="bottom")
|
||||
ax.set_ylim(-0.4, 4.9)
|
||||
|
||||
# --- opt3: no barriers, continuous
|
||||
ax = axes[2]
|
||||
for i in range(11):
|
||||
frame(ax, 3 - (i % 4) * 1.0, i * 11)
|
||||
ax.set_ylim(-1.5, 4.5)
|
||||
ax.text(0, -0.25, "streams never wait on each other — the GPU is continuously busy",
|
||||
color=ACCENT, fontsize=7.5, va="top")
|
||||
|
||||
titles = ["opt1 · 1 stream, synchronous",
|
||||
"opt2 · 4 streams, sync barrier per round",
|
||||
"opt3 · 4 streams, barriers removed"]
|
||||
for ax, t in zip(axes, titles):
|
||||
ax.set_xlim(-2, 190)
|
||||
ax.set_yticks([]); ax.set_xticks([])
|
||||
bare(ax, keep=())
|
||||
ax.set_title(t, color=TEXT2, fontsize=9, loc="left", pad=4)
|
||||
|
||||
handles = [Rectangle((0, 0), 1, 1, color=c) for c in (AMBER, ACCENT, PALE)]
|
||||
axes[0].legend(handles, ["H2D copy", "kernel", "D2H copy"], frameon=False,
|
||||
fontsize=8, labelcolor=TEXT2, ncol=3, loc="lower right",
|
||||
bbox_to_anchor=(1.02, 0.98), handlelength=1.1)
|
||||
axes[2].set_xlabel("time →", color=MUTED, fontsize=8.5, loc="left")
|
||||
fig.subplots_adjust(hspace=0.75)
|
||||
save(fig, "fig_streams")
|
||||
|
||||
|
||||
# ------------------------------------------------------------- 4. pinning
|
||||
def fig_pinning():
|
||||
fig = plt.figure(figsize=(7.7, 3.0))
|
||||
ax = fig.add_axes([0, 0.05, 0.63, 0.95]); ax.axis("off")
|
||||
ax.set_xlim(0, 10.4); ax.set_ylim(0, 6.4)
|
||||
|
||||
def box(x, y, w, h, label, sub=""):
|
||||
ax.add_patch(Rectangle((x, y), w, h, facecolor=PANEL, edgecolor=RULE, lw=1))
|
||||
ax.text(x + w / 2, y + h / 2 + 0.26, label, ha="center", va="center",
|
||||
color=PALE, fontsize=8.5, fontweight="bold")
|
||||
ax.text(x + w / 2, y + h / 2 - 0.34, sub, ha="center", va="center",
|
||||
color=MUTED, fontsize=7)
|
||||
|
||||
def arrow(x0, x1, y, color, label):
|
||||
ax.add_patch(FancyArrowPatch((x0, y), (x1, y), arrowstyle="-|>",
|
||||
mutation_scale=10, color=color, lw=1.6))
|
||||
ax.text((x0 + x1) / 2, y + 0.22, label, ha="center", va="bottom",
|
||||
color=color, fontsize=7)
|
||||
|
||||
ax.text(0, 5.95, "PAGEABLE · before opt4", color=AMBER, fontsize=8.5,
|
||||
fontweight="bold")
|
||||
box(0, 4.05, 2.5, 1.1, "numpy array", "pageable")
|
||||
box(4.0, 4.05, 2.4, 1.1, "driver staging", "hidden pinned buf")
|
||||
box(7.9, 4.05, 2.5, 1.1, "GPU", "device memory")
|
||||
arrow(2.5, 4.0, 4.60, AMBER, "memcpy")
|
||||
arrow(6.4, 7.9, 4.60, AMBER, "DMA")
|
||||
ax.text(0, 3.62, "every transfer is copied twice", color=MUTED, fontsize=7)
|
||||
|
||||
ax.text(0, 2.75, "PINNED · opt4", color=ACCENT, fontsize=8.5, fontweight="bold")
|
||||
box(0, 0.85, 2.5, 1.1, "numpy array", "page-locked")
|
||||
box(7.9, 0.85, 2.5, 1.1, "GPU", "device memory")
|
||||
arrow(2.5, 7.9, 1.40, ACCENT, "DMA — engine reads host RAM directly")
|
||||
ax.text(0, 0.42, "no staging copy, no page faults, fully async",
|
||||
color=MUTED, fontsize=7)
|
||||
|
||||
ax2 = fig.add_axes([0.75, 0.16, 0.25, 0.66])
|
||||
v = [14, 3]
|
||||
ax2.bar([0, 1], v, width=0.55, color=[AMBER, ACCENT], zorder=3)
|
||||
for i, val in enumerate(v):
|
||||
ax2.text(i, val + 0.5, f"{val} µs", ha="center", color=PALE,
|
||||
fontsize=10, fontweight="bold")
|
||||
ax2.set_xticks([0, 1])
|
||||
ax2.set_xticklabels(["opt3\npageable", "opt4\npinned"], color=TEXT2, fontsize=8)
|
||||
ax2.set_ylim(0, 18); ax2.set_yticks([]); bare(ax2, keep=("bottom",))
|
||||
ax2.set_title("host overhead / frame", color=MUTED, fontsize=7.5, pad=6)
|
||||
save(fig, "fig_pinning")
|
||||
|
||||
|
||||
# -------------------------------------------------------------- 5. graphs
|
||||
def fig_graphs():
|
||||
fig, ax = plt.subplots(figsize=(7.7, 2.6))
|
||||
ax.axis("off"); ax.set_xlim(0, 12.6); ax.set_ylim(0, 4.6)
|
||||
|
||||
def node(x, y, w, h, t, fc):
|
||||
ax.add_patch(Rectangle((x, y), w, h, facecolor=fc, edgecolor="none"))
|
||||
ax.text(x + w / 2, y + h / 2, t, ha="center", va="center",
|
||||
color=BG, fontsize=7.5, fontweight="bold")
|
||||
|
||||
ops = [("H2D", AMBER), ("kernel", ACCENT), ("D2H", PALE)] * 2
|
||||
|
||||
ax.text(0, 4.15, "WITHOUT GRAPHS · one driver call per operation, every frame",
|
||||
color=AMBER, fontsize=8.5, fontweight="bold")
|
||||
for i, (t, c) in enumerate(ops):
|
||||
x = 0.1 + i * 1.62
|
||||
node(x, 2.85, 1.4, 0.6, t, c)
|
||||
ax.add_patch(FancyArrowPatch((x + 0.7, 3.72), (x + 0.7, 3.52),
|
||||
arrowstyle="-|>", mutation_scale=7,
|
||||
color=MUTED, lw=0.9))
|
||||
ax.text(12.5, 3.15, "CPU cost\n≈ 6 launches", ha="right", va="center",
|
||||
color=MUTED, fontsize=7.5)
|
||||
|
||||
ax.text(0, 2.18, "WITH GRAPHS · opt5 · record once, replay with one launch",
|
||||
color=ACCENT, fontsize=8.5, fontweight="bold")
|
||||
ax.add_patch(Rectangle((0.1, 0.72), 9.85, 1.15, facecolor=PANEL,
|
||||
edgecolor=ACCENT, lw=1.2))
|
||||
for i, (t, c) in enumerate(ops):
|
||||
node(0.38 + i * 1.58, 0.98, 1.34, 0.6, t, c)
|
||||
ax.add_patch(FancyArrowPatch((0.8, 2.02), (0.8, 1.90), arrowstyle="-|>",
|
||||
mutation_scale=8, color=ACCENT, lw=1.3))
|
||||
ax.text(12.5, 1.30, "CPU cost\n≈ 1 launch", ha="right", va="center",
|
||||
color=ACCENT, fontsize=7.5, fontweight="bold")
|
||||
ax.text(0.1, 0.32, "cudaGraphLaunch() — the whole DAG is submitted as one unit; "
|
||||
"the driver already knows every dependency",
|
||||
color=MUTED, fontsize=7)
|
||||
save(fig, "fig_graphs")
|
||||
|
||||
|
||||
# ------------------------------------------------ 6. f32 kernel (nsys truth)
|
||||
def fig_f32_kernel():
|
||||
fig, (ax, ax2) = plt.subplots(1, 2, figsize=(7.7, 2.6),
|
||||
gridspec_kw={"width_ratios": [1, 1.35]})
|
||||
v = [43.0, 25.6]
|
||||
ax.bar([0, 1], v, width=0.5, color=[AMBER, ACCENT], zorder=3)
|
||||
for i, val in enumerate(v):
|
||||
ax.text(i, val + 1.2, f"{val} µs", ha="center", color=PALE,
|
||||
fontsize=11, fontweight="bold")
|
||||
ax.annotate("", xy=(1, 27.5), xytext=(0, 44.5),
|
||||
arrowprops=dict(arrowstyle="-|>", color=MUTED, lw=1.2,
|
||||
connectionstyle="arc3,rad=-0.25"))
|
||||
ax.text(0.5, 37, "−40%", ha="center", color=PALE, fontsize=10,
|
||||
fontweight="bold")
|
||||
ax.set_xticks([0, 1]); ax.set_xticklabels(["f64 pedestal", "f32 pedestal"],
|
||||
color=TEXT2, fontsize=8.5)
|
||||
ax.set_ylim(0, 52); ax.set_yticks([]); bare(ax, keep=("bottom",))
|
||||
ax.set_title("kernel, exclusive (nsys, 9×9)", color=MUTED, fontsize=8, pad=8)
|
||||
|
||||
labels = ["kernel", "D2H", "H2D"]
|
||||
f64 = [43.0, 19.4, 13.2]
|
||||
f32 = [25.6, 19.8, 13.5]
|
||||
y = np.arange(3); h = 0.35
|
||||
ax2.barh(y + h / 2, f64, height=h, color=AMBER, zorder=3, label="f64 ped")
|
||||
ax2.barh(y - h / 2, f32, height=h, color=ACCENT, zorder=3, label="f32 ped")
|
||||
for yi, (a, b) in enumerate(zip(f64, f32)):
|
||||
ax2.text(a + 1, yi + h / 2, f"{a:.1f}", va="center", color=TEXT2, fontsize=8)
|
||||
ax2.text(b + 1, yi - h / 2, f"{b:.1f}", va="center", color=TEXT2, fontsize=8)
|
||||
ax2.set_yticks(y); ax2.set_yticklabels(labels, color=TEXT2, fontsize=8.5)
|
||||
ax2.invert_yaxis(); ax2.set_xlim(0, 56); ax2.set_xticks([])
|
||||
bare(ax2, keep=("left",))
|
||||
ax2.legend(frameon=False, fontsize=8, labelcolor=TEXT2, loc="lower right")
|
||||
ax2.set_title("per-frame GPU operations (µs)", color=MUTED, fontsize=8, pad=8)
|
||||
save(fig, "fig_f32_kernel")
|
||||
|
||||
|
||||
# ------------------------------------------------------- 7. cancellation
|
||||
def fig_cancellation():
|
||||
fig, (ax, ax2) = plt.subplots(1, 2, figsize=(7.7, 2.7),
|
||||
gridspec_kw={"width_ratios": [1.25, 1]})
|
||||
names = ["E[X²]\n2.17e7", "mean²\n2.17e7", "variance\n2025"]
|
||||
vals = [2.17e7, 2.17e7, 2025]
|
||||
ax.bar([0, 1], vals[:2], width=0.5, color=[PALE, PALE], zorder=3)
|
||||
ax.bar([2], [2025], width=0.5, color=AMBER, zorder=3)
|
||||
ax.set_yscale("log"); ax.set_ylim(1e2, 2e8)
|
||||
ax.set_xticks([0, 1, 2]); ax.set_xticklabels(names, color=TEXT2, fontsize=8)
|
||||
ax.set_yticks([1e3, 1e5, 1e7])
|
||||
ax.axhline(2048, color=ACCENT, lw=1.3, ls="--", zorder=4)
|
||||
ax.text(2.42, 3000, "f32 rounding step\nat 2.17e7 = 2048", color=ACCENT,
|
||||
fontsize=7.5, ha="right", va="bottom")
|
||||
bare(ax)
|
||||
ax.set_title("var = E[X²] − mean² (f32, mean ≈ 4655 ADU)",
|
||||
color=MUTED, fontsize=8, pad=8)
|
||||
|
||||
rms = np.linspace(0, 12, 200)
|
||||
ax2.fill_between(rms, 0, np.where(rms < 6.5, 1, 0), color=AMBER, alpha=0.16,
|
||||
step="pre")
|
||||
ax2.plot(rms, rms**2, color=PALE, lw=1.8, label="true variance")
|
||||
ax2.axhline(42, color=ACCENT, lw=1.4, ls="--", label="f32 error floor")
|
||||
ax2.set_xlabel("pixel rms (ADU)", color=TEXT2, fontsize=8.5)
|
||||
ax2.set_ylabel("variance", color=TEXT2, fontsize=8.5)
|
||||
ax2.set_ylim(0, 150); ax2.set_xlim(0, 12)
|
||||
ax2.set_yticks([]); ax2.tick_params(labelsize=8)
|
||||
bare(ax2)
|
||||
ax2.text(1.0, 108, "quiet pixels:\nerror > variance\n→ rms clamped to 0\n→ fires every frame",
|
||||
color=AMBER, fontsize=7.5, va="top")
|
||||
ax2.legend(frameon=False, fontsize=7.5, labelcolor=TEXT2, loc="lower right")
|
||||
save(fig, "fig_cancellation")
|
||||
|
||||
|
||||
# ------------------------------------------------- 8. where f32 pays or not
|
||||
def fig_bottleneck():
|
||||
fig, (a1, a2) = plt.subplots(1, 2, figsize=(11.4, 3.0))
|
||||
|
||||
for ax, title, kern, floor, gain in [
|
||||
(a1, "3×3 clusters — pipeline-bound", (24, 13), 25,
|
||||
"kernel already hidden → 0% end-to-end"),
|
||||
(a2, "9×9 clusters — kernel-bound (f64)", (43, 26), 32,
|
||||
"kernel on the critical path → −8% wall"),
|
||||
]:
|
||||
x = [0, 1]
|
||||
ax.bar(x, kern, width=0.5, color=[AMBER, ACCENT], zorder=3)
|
||||
ax.axhline(floor, color=PALE, lw=1.4, ls="--", zorder=4)
|
||||
ax.text(1.62, floor + 1.2, "transfer + host floor", color=PALE,
|
||||
fontsize=8, ha="right")
|
||||
for i, v in enumerate(kern):
|
||||
ax.text(i, v + 1.2, f"{v} µs", ha="center", color=PALE,
|
||||
fontsize=10, fontweight="bold")
|
||||
ax.set_xticks(x); ax.set_xticklabels(["f64 pedestal", "f32 pedestal"],
|
||||
color=TEXT2, fontsize=9)
|
||||
ax.set_xlim(-0.6, 1.7); ax.set_ylim(0, 55); ax.set_yticks([])
|
||||
bare(ax, keep=("bottom",))
|
||||
ax.set_title(title, color=PALE, fontsize=9.5, pad=10)
|
||||
ax.text(-0.55, -9, gain, color=AMBER if "0%" in gain else ACCENT,
|
||||
fontsize=8.5, fontweight="bold")
|
||||
fig.subplots_adjust(bottom=0.22)
|
||||
save(fig, "fig_bottleneck")
|
||||
|
||||
|
||||
# ----------------------------------------------------------- 9. correctness
|
||||
def fig_correctness():
|
||||
fig, ax = plt.subplots(figsize=(7.4, 2.4))
|
||||
names = ["CPU MT", "opt1", "opt2", "opt3", "opt4", "opt5", "opt6 (f32)"]
|
||||
diff = [0.0, 0.0040, 0.0035, 0.0035, 0.0035, 0.0035, 0.0039]
|
||||
colors = [MUTED] + [ACCENT] * 5 + [AMBER]
|
||||
x = np.arange(len(names))
|
||||
ax.bar(x, diff, width=0.55, color=colors, zorder=3)
|
||||
for xi, d in enumerate(diff):
|
||||
ax.text(xi, d + 0.00022, ("reference" if d == 0 else f"{d:.4f}%"),
|
||||
ha="center", color=PALE if d else MUTED, fontsize=8.5)
|
||||
ax.axhline(0.01, color=PALE, lw=1.2, ls="--")
|
||||
ax.text(6.4, 0.0104, "0.01% — well inside statistical noise", color=PALE,
|
||||
fontsize=8, ha="right")
|
||||
ax.set_xticks(x); ax.set_xticklabels(names, color=TEXT2, fontsize=8.5)
|
||||
ax.set_ylim(0, 0.0125); ax.set_yticks([])
|
||||
bare(ax, keep=("bottom",))
|
||||
ax.set_title("cluster-count difference vs CPU · 233 million clusters, 3×3",
|
||||
color=MUTED, fontsize=8.5, pad=8)
|
||||
save(fig, "fig_correctness")
|
||||
|
||||
|
||||
for f in (fig_arc, fig_overhead, fig_streams, fig_pinning, fig_graphs,
|
||||
fig_f32_kernel, fig_cancellation, fig_bottleneck, fig_correctness):
|
||||
f()
|
||||
print("done ->", OUT)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 66 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 56 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 73 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 39 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 48 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 55 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 33 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 72 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 54 KiB |
@@ -0,0 +1,435 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#pragma once
|
||||
#include "aare/ClusterFinder.hpp"
|
||||
#include "aare/clusterfinder_kernel_opt2.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 StreamContextOpt2 {
|
||||
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 ClusterFinderCUDAOpt2 {
|
||||
using COMPUTE_TYPE =
|
||||
device_opt2::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 = StreamContextOpt2<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 ClusterFinderCUDAOpt2
|
||||
*
|
||||
* @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
|
||||
*/
|
||||
ClusterFinderCUDAOpt2(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(
|
||||
"ClusterFinderCUDAOpt2: n_streams must be > 0");
|
||||
}
|
||||
|
||||
if (capacity >
|
||||
static_cast<size_t>(std::numeric_limits<uint32_t>::max())) {
|
||||
throw std::invalid_argument(
|
||||
"ClusterFinderCUDAOpt2: capacity must fit in uint32_t");
|
||||
}
|
||||
|
||||
if (capacity == 0) {
|
||||
throw std::invalid_argument(
|
||||
"ClusterFinderCUDAOpt2: 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));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
~ClusterFinderCUDAOpt2() {
|
||||
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
|
||||
ClusterFinderCUDAOpt2(const ClusterFinderCUDAOpt2 &) = delete;
|
||||
ClusterFinderCUDAOpt2 &operator=(const ClusterFinderCUDAOpt2 &) = delete;
|
||||
ClusterFinderCUDAOpt2(ClusterFinderCUDAOpt2 &&) = delete;
|
||||
ClusterFinderCUDAOpt2 &operator=(ClusterFinderCUDAOpt2 &&) = 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_opt2::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_opt2::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
|
||||
@@ -0,0 +1,346 @@
|
||||
#pragma once
|
||||
#include "aare/Cluster.hpp"
|
||||
#include "aare/ClusterFinder.hpp"
|
||||
#include <cuda_runtime.h>
|
||||
#include <type_traits>
|
||||
|
||||
// OPT2 SNAPSHOT (commit 88e0e8d) — the pre-refactor kernel, kept for
|
||||
// benchmarking the opt2 pipeline. Lives in namespace aare::device_opt2 so it
|
||||
// coexists with the modern kernel. opt2 pipeline behaviour: f32 stencil /
|
||||
// f64 pedestal, raw E[X^2]-E[X]^2 variance. The Test3 local-max gate has been
|
||||
// backported from the current kernel so opt1/opt2 cluster counts match the CPU
|
||||
// and current finders (holding correctness constant across the opt arc).
|
||||
namespace aare::device_opt2 {
|
||||
|
||||
// Implementing mixed precision for shared memory and stencil arithmetic
|
||||
using COMPUTE_TYPE = float;
|
||||
|
||||
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>>
|
||||
__global__ void find_clusters_in_single_frame(
|
||||
const FRAME_TYPE *__restrict__ d_frame,
|
||||
PEDESTAL_TYPE *__restrict__ d_pd_mean, PEDESTAL_TYPE *__restrict__ d_pd_sum,
|
||||
PEDESTAL_TYPE *__restrict__ d_pd_sum2, const uint32_t n_pd_samples,
|
||||
const COMPUTE_TYPE m_nSigma, const size_t nrows, const size_t ncols,
|
||||
// const uint64_t frame_number,
|
||||
ClusterType *d_clusters, uint32_t *d_cluster_count,
|
||||
const uint32_t max_clusters) {
|
||||
using CT = typename ClusterType::value_type;
|
||||
|
||||
// Compile-time cluster geometry useful for unrolling loops
|
||||
constexpr uint8_t CSX = ClusterType::cluster_size_x;
|
||||
constexpr uint8_t CSY = ClusterType::cluster_size_y;
|
||||
constexpr int col_radius = CSX / 2;
|
||||
constexpr int row_radius = CSY / 2;
|
||||
|
||||
// Squared threshold constants; avoids sqrt at runtime
|
||||
// c2^2 is for the 2x2 quadrant test
|
||||
// c3^2 is for the full-cluster total test
|
||||
constexpr int pow2_c2 = ((CSY + 1) / 2) * ((CSX + 1) / 2);
|
||||
constexpr int pow2_c3 = CSX * CSY;
|
||||
|
||||
// Thread/pixel mapping
|
||||
auto col_global =
|
||||
static_cast<ssize_t>(threadIdx.x + blockDim.x * blockIdx.x);
|
||||
auto row_global =
|
||||
static_cast<ssize_t>(threadIdx.y + blockDim.y * blockIdx.y);
|
||||
auto global_tid = static_cast<ssize_t>(col_global + ncols * row_global);
|
||||
auto local_tid = threadIdx.x + blockDim.x * threadIdx.y;
|
||||
|
||||
// ====================
|
||||
// Shared memory layout
|
||||
// ====================
|
||||
// The tile is laid out contiguously in a 1D configuration:
|
||||
// [0 ... tile_size-1] = pedestal-subtracted frame values
|
||||
//
|
||||
// Each tile has (blockDim.x + 2*col_radius) x (blockDim.y + 2*row_radius)
|
||||
// elements, with a halo of col_radius/row_radius pixels on each side.
|
||||
|
||||
// CUDA prefers raw bytes + aligned cast
|
||||
// Compile error happens when using: `extern __shared__ T sh[];`
|
||||
extern __shared__ __align__(sizeof(COMPUTE_TYPE)) unsigned char smem[];
|
||||
COMPUTE_TYPE *shmem = reinterpret_cast<COMPUTE_TYPE *>(smem);
|
||||
|
||||
// Stride includes halo on both sides
|
||||
auto shmem_stride = static_cast<int>(blockDim.x) + 2 * col_radius;
|
||||
auto tile_size =
|
||||
shmem_stride * (static_cast<int>(blockDim.y) + 2 * row_radius);
|
||||
|
||||
// Offset so that thread (0,0) maps to shared-memory position
|
||||
// (row_radius, col_radius) i.e. past the top-left halo.
|
||||
auto shmem_tid =
|
||||
(static_cast<int>(threadIdx.y) + row_radius) * shmem_stride +
|
||||
(static_cast<int>(threadIdx.x) + col_radius);
|
||||
|
||||
// Cooperative zero-fill
|
||||
for (int idx = static_cast<int>(local_tid); idx < tile_size;
|
||||
idx += static_cast<int>(blockDim.x * blockDim.y)) {
|
||||
shmem[idx] = COMPUTE_TYPE{0};
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
// OOB flag
|
||||
bool valid_pixel = col_global < static_cast<ssize_t>(ncols) &&
|
||||
row_global < static_cast<ssize_t>(nrows);
|
||||
|
||||
// ======================================================
|
||||
// Load pedestal-subtracted frame data into shared memory (MIXED PRECISION)
|
||||
// ======================================================
|
||||
|
||||
// Helper: read (frame - pedestal_mean) from global memory, or 0 if OOB.
|
||||
// gr, gc are the global row/col of the pixel to load.
|
||||
// Returns the pedestal-subtracted value.
|
||||
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]);
|
||||
};
|
||||
|
||||
// A. Interior: every valid thread loads its own pixel
|
||||
if (valid_pixel) {
|
||||
shmem[shmem_tid] = load_pixel(row_global, col_global);
|
||||
}
|
||||
|
||||
// B. Halo regions (Boundaries)
|
||||
// B.1 Top rows of the halo
|
||||
if (threadIdx.y == 0 && valid_pixel) {
|
||||
for (int i = 1; i <= row_radius; ++i) {
|
||||
if (row_global - i >= 0)
|
||||
shmem[shmem_tid - i * shmem_stride] =
|
||||
load_pixel(row_global - i, col_global);
|
||||
}
|
||||
// Top-left corner rectangle
|
||||
if (threadIdx.x == 0) {
|
||||
for (int i = 1; i <= row_radius; ++i)
|
||||
for (int j = 1; j <= col_radius; ++j)
|
||||
if (row_global - i >= 0 && col_global - j >= 0)
|
||||
shmem[shmem_tid - i * shmem_stride - j] =
|
||||
load_pixel(row_global - i, col_global - j);
|
||||
}
|
||||
// Top-right corner rectangle
|
||||
if (threadIdx.x == blockDim.x - 1) {
|
||||
for (int i = 1; i <= row_radius; ++i)
|
||||
for (int j = 1; j <= col_radius; ++j)
|
||||
if (row_global - i >= 0 &&
|
||||
col_global + j < static_cast<ssize_t>(ncols))
|
||||
shmem[shmem_tid - i * shmem_stride + j] =
|
||||
load_pixel(row_global - i, col_global + j);
|
||||
}
|
||||
}
|
||||
|
||||
// B.2 Left column of the halo
|
||||
if (threadIdx.x == 0 && valid_pixel) {
|
||||
for (int j = 1; j <= col_radius; ++j)
|
||||
if (col_global - j >= 0)
|
||||
shmem[shmem_tid - j] = load_pixel(row_global, col_global - j);
|
||||
}
|
||||
|
||||
// B.3 Right column of the halo
|
||||
if (threadIdx.x == blockDim.x - 1 && valid_pixel) {
|
||||
for (int j = 1; j <= col_radius; ++j)
|
||||
if (col_global + j < static_cast<ssize_t>(ncols))
|
||||
shmem[shmem_tid + j] = load_pixel(row_global, col_global + j);
|
||||
}
|
||||
|
||||
// B.4 Bottom rows of the halo
|
||||
if (threadIdx.y == blockDim.y - 1 && valid_pixel) {
|
||||
for (int i = 1; i <= row_radius; ++i) {
|
||||
if (row_global + i < static_cast<ssize_t>(nrows))
|
||||
shmem[shmem_tid + i * shmem_stride] =
|
||||
load_pixel(row_global + i, col_global);
|
||||
}
|
||||
// Bottom-left corner rectangle
|
||||
if (threadIdx.x == 0) {
|
||||
for (int i = 1; i <= row_radius; ++i)
|
||||
for (int j = 1; j <= col_radius; ++j)
|
||||
if (row_global + i < static_cast<ssize_t>(nrows) &&
|
||||
col_global - j >= 0)
|
||||
shmem[shmem_tid + i * shmem_stride - j] =
|
||||
load_pixel(row_global + i, col_global - j);
|
||||
}
|
||||
// Bottom-right corner rectangle
|
||||
if (threadIdx.x == blockDim.x - 1) {
|
||||
for (int i = 1; i <= row_radius; ++i)
|
||||
for (int j = 1; j <= col_radius; ++j)
|
||||
if (row_global + i < static_cast<ssize_t>(nrows) &&
|
||||
col_global + j < static_cast<ssize_t>(ncols))
|
||||
shmem[shmem_tid + i * shmem_stride + j] =
|
||||
load_pixel(row_global + i, col_global + j);
|
||||
}
|
||||
}
|
||||
|
||||
__syncthreads();
|
||||
|
||||
// =====================
|
||||
// Cluster-finding logic
|
||||
// =====================
|
||||
if (!valid_pixel)
|
||||
return;
|
||||
|
||||
// Per-pixel variance from global pedestal arrays
|
||||
// Variance = rms^2 = E[X^2] - E[X]^2
|
||||
// NOTE: Keep thresholds squared to avoid one sqrtf() per pixel.
|
||||
PEDESTAL_TYPE mean_px = d_pd_mean[global_tid];
|
||||
PEDESTAL_TYPE var_px =
|
||||
d_pd_sum2[global_tid] / n_pd_samples - mean_px * mean_px;
|
||||
PEDESTAL_TYPE rms_sq = max(var_px, PEDESTAL_TYPE{0}); // variance = rms^2
|
||||
PEDESTAL_TYPE nSig_sq_rms_sq = static_cast<PEDESTAL_TYPE>(m_nSigma) *
|
||||
static_cast<PEDESTAL_TYPE>(m_nSigma) *
|
||||
rms_sq;
|
||||
|
||||
// Pedestal-subtracted value of the center pixel (already in shmem)
|
||||
COMPUTE_TYPE val_pixel = shmem[shmem_tid];
|
||||
|
||||
// Negative pedestal early exit:
|
||||
// val_pixel < -nSigma * rms
|
||||
// is equivalent to:
|
||||
// val_pixel < 0 && val_pixel^2 > nSigma^2 * rms^2
|
||||
if (val_pixel < COMPUTE_TYPE{0} &&
|
||||
static_cast<PEDESTAL_TYPE>(val_pixel) *
|
||||
static_cast<PEDESTAL_TYPE>(val_pixel) >
|
||||
nSig_sq_rms_sq)
|
||||
return; // NOTE: pedestal update for this pixel is skipped (same as
|
||||
// sequential)
|
||||
|
||||
// Stencil reduction: total, max, quadrant sums
|
||||
COMPUTE_TYPE total = 0.0f;
|
||||
COMPUTE_TYPE max_val = -HUGE_VALF;
|
||||
|
||||
// // Quandrants
|
||||
// PEDESTAL_TYPE tl = PEDESTAL_TYPE{0}; // top-left quadrant (ir<=0,
|
||||
// ic<=0) PEDESTAL_TYPE tr = PEDESTAL_TYPE{0}; // top-right quadrant
|
||||
// (ir<=0, ic>=0) PEDESTAL_TYPE bl = PEDESTAL_TYPE{0}; // bottom-left
|
||||
// (ir>=0, ic<=0) PEDESTAL_TYPE br = PEDESTAL_TYPE{0}; // bottom-right
|
||||
// (ir>=0, ic>=0)
|
||||
|
||||
#pragma unroll
|
||||
for (int ir = -row_radius; ir <= row_radius; ++ir) {
|
||||
#pragma unroll
|
||||
for (int ic = -col_radius; ic <= col_radius; ++ic) {
|
||||
COMPUTE_TYPE val = shmem[shmem_tid + ir * shmem_stride + ic];
|
||||
|
||||
total += val;
|
||||
max_val = fmaxf(max_val, val);
|
||||
|
||||
// // Quadrant accumulation (pixels on the axes contribute to two
|
||||
// quadrants) if (ir <= 0 && ic <= 0) tl += val; if (ir <= 0 && ic
|
||||
// >= 0) tr += val; if (ir >= 0 && ic <= 0) bl += val; if (ir >= 0
|
||||
// && ic >= 0) br += val;
|
||||
}
|
||||
}
|
||||
|
||||
// Three-way classification (mirrors ClusterFinder's logic)
|
||||
//
|
||||
// 1. Single-pixel significance: max_val > nSigma * rms
|
||||
// -> only the pixel that IS the max gets recorded (local-max
|
||||
// suppression)
|
||||
//
|
||||
// 2. Quadrant significance: max(tl,tr,bl,br) > c2 * nSigma * rms
|
||||
// -> charge-sharing events where a 2x2 sub-region is significant
|
||||
// NOTE: This test is absent in the serial ClusterFinder!
|
||||
//
|
||||
// 3. Total significance: total > c3 * nSigma * rms
|
||||
// -> distributed events where the full cluster sum is significant
|
||||
|
||||
bool is_photon = false;
|
||||
|
||||
// Test 1: single-pixel significance
|
||||
// max_val > nSigma * rms
|
||||
// is equivalent to:
|
||||
// max_val > 0 && max_val^2 > nSigma^2 * rms^2
|
||||
if (max_val > COMPUTE_TYPE{0} &&
|
||||
static_cast<PEDESTAL_TYPE>(max_val) *
|
||||
static_cast<PEDESTAL_TYPE>(max_val) >
|
||||
nSig_sq_rms_sq) {
|
||||
// Local-max suppression: only the center-pixel thread records the
|
||||
// cluster
|
||||
if (val_pixel < max_val)
|
||||
return; // some other pixel in the neighborhood is brighter
|
||||
is_photon = true;
|
||||
}
|
||||
|
||||
/* // Test 2: quadrant significance (only if test 1 didn't fire)
|
||||
if (!is_photon) {
|
||||
PEDESTAL_TYPE max_quad = max(max(tl, tr), max(bl, br));
|
||||
if (max_quad > 0 && max_quad * max_quad > pow2_c2 * nSig_sq_rms_sq) {
|
||||
is_photon = true;
|
||||
}
|
||||
} */
|
||||
|
||||
// Test 3: total significance (only if tests 1 & 2 didn't fire)
|
||||
if (!is_photon) {
|
||||
if (total > 0 && static_cast<PEDESTAL_TYPE>(total) *
|
||||
static_cast<PEDESTAL_TYPE>(total) >
|
||||
pow2_c3 * nSig_sq_rms_sq) {
|
||||
// Local-max suppression: only the center-pixel thread records the
|
||||
// cluster, so an extended charge-shared event yields one cluster,
|
||||
// not one per pixel. (Backported from the current kernel to make
|
||||
// the opt2 counts match; the historical snapshot lacked this gate.)
|
||||
if (val_pixel < max_val)
|
||||
return; // some other pixel in the neighborhood is brighter
|
||||
is_photon = true;
|
||||
}
|
||||
}
|
||||
|
||||
// Pedestal update (if not a photon)
|
||||
// In the sequential code, non-photon pixels feed back into the running
|
||||
// pedestal via push_fast(). In this kernel, the GPU updates all pixels in a
|
||||
// frame simultaneously. So the updated pedestal will only be used starting
|
||||
// from the next frame. -> This avoids a/serialization and b/global mem I/O.
|
||||
if (!is_photon && valid_pixel) {
|
||||
PEDESTAL_TYPE raw_val = static_cast<PEDESTAL_TYPE>(d_frame[global_tid]);
|
||||
PEDESTAL_TYPE sum = d_pd_sum[global_tid];
|
||||
PEDESTAL_TYPE sum2 = d_pd_sum2[global_tid];
|
||||
|
||||
sum += raw_val - sum / n_pd_samples;
|
||||
sum2 += raw_val * raw_val - sum2 / n_pd_samples;
|
||||
|
||||
d_pd_sum[global_tid] = sum;
|
||||
d_pd_sum2[global_tid] = sum2;
|
||||
d_pd_mean[global_tid] = sum / n_pd_samples;
|
||||
return;
|
||||
}
|
||||
|
||||
/*
|
||||
if (!is_photon) return; // Debugging
|
||||
*/
|
||||
|
||||
// Delay building clusterData until we know this thread will write a photon.
|
||||
// This avoids CSX*CSY conversions/rounds for the overwhelmingly common
|
||||
// background pixels.
|
||||
CT clusterData[CSX * CSY];
|
||||
int idx = 0;
|
||||
|
||||
#pragma unroll
|
||||
for (int ir = -row_radius; ir <= row_radius; ++ir) {
|
||||
#pragma unroll
|
||||
for (int ic = -col_radius; ic <= col_radius; ++ic) {
|
||||
COMPUTE_TYPE val = shmem[shmem_tid + ir * shmem_stride + ic];
|
||||
if constexpr (std::is_integral_v<CT>)
|
||||
clusterData[idx] = static_cast<CT>(lroundf(val));
|
||||
else
|
||||
clusterData[idx] = static_cast<CT>(val);
|
||||
idx++;
|
||||
}
|
||||
}
|
||||
|
||||
// Write cluster to global output buffer using atomic index
|
||||
// for coordination across all blocks
|
||||
uint32_t write_idx = atomicAdd(d_cluster_count, 1u);
|
||||
|
||||
// Guard against overflowing the pre-allocated cluster buffer
|
||||
if (write_idx >= max_clusters)
|
||||
return;
|
||||
|
||||
ClusterType cluster{};
|
||||
cluster.x = static_cast<decltype(cluster.x)>(col_global);
|
||||
cluster.y = static_cast<decltype(cluster.y)>(row_global);
|
||||
|
||||
memcpy(reinterpret_cast<CT *>(&cluster.data), clusterData,
|
||||
sizeof(CT) * CSX * CSY);
|
||||
|
||||
d_clusters[write_idx] = cluster;
|
||||
}
|
||||
|
||||
} // namespace aare::device_opt2
|
||||
@@ -166,7 +166,34 @@ def ClusterFinderCUDAGraph(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.i
|
||||
n_streams=n_streams)
|
||||
|
||||
|
||||
def ClusterCollector(clusterfindermt, dtype=np.int32):
|
||||
def ClusterFinderCUDAOpt2(image_size, cluster_size=(3, 3), n_sigma=5, dtype=np.int32,
|
||||
max_clusters_per_frame=3000, n_streams=4):
|
||||
"""
|
||||
Factory for the OPT2 snapshot finder — the pre-refactor pipeline (per-frame
|
||||
pinned staging, round-robin streams with sync barriers, variable-length D2H),
|
||||
kept only for benchmarking the optimization arc against the current finder.
|
||||
|
||||
It uses its own kernel snapshot (clusterfinder_kernel_opt2.cuh): f32 stencil
|
||||
/ f64 pedestal. The Test3 local-max gate has been backported so its cluster
|
||||
counts match the CPU and current finders (correctness held constant across
|
||||
the opt arc; only the pipeline differs).
|
||||
|
||||
Only the 3x3 cluster size is registered.
|
||||
"""
|
||||
if not _cuda_available():
|
||||
raise RuntimeError(
|
||||
"ClusterFinderCUDAOpt2 is not available in this build of aare. "
|
||||
"Rebuild with -DAARE_CUDA=ON (and -DAARE_PYTHON_BINDINGS=ON)."
|
||||
)
|
||||
|
||||
cls = _get_class("ClusterFinderCUDAOpt2", cluster_size, dtype)
|
||||
return cls(image_size,
|
||||
n_sigma=n_sigma,
|
||||
max_clusters_per_frame=max_clusters_per_frame,
|
||||
n_streams=n_streams)
|
||||
|
||||
|
||||
def ClusterCollector(clusterfindermt, dtype=np.int32):
|
||||
"""
|
||||
Factory function to create a ClusterCollector object. Provides a cleaner syntax for
|
||||
the templated ClusterCollector in C++.
|
||||
|
||||
@@ -32,7 +32,7 @@ from ._aare import corner
|
||||
|
||||
from ._version import __version__
|
||||
from .ClusterFinder import ClusterFinder, ClusterFinderFrozen, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
|
||||
from .ClusterFinder import ClusterFinderCUDA, ClusterFinderCUDAGraph, _cuda_available
|
||||
from .ClusterFinder import ClusterFinderCUDA, ClusterFinderCUDAGraph, ClusterFinderCUDAOpt2, _cuda_available
|
||||
from .ClusterVector import ClusterVector
|
||||
from .Cluster import Cluster
|
||||
|
||||
|
||||
@@ -0,0 +1,102 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#pragma once
|
||||
#include "aare/ClusterFinderCUDAOpt2.hpp"
|
||||
#include "aare/ClusterVector.hpp"
|
||||
#include "aare/NDView.hpp"
|
||||
#include "aare/Pedestal.hpp"
|
||||
#include "np_helper.hpp"
|
||||
|
||||
#include <cstdint>
|
||||
#include <pybind11/pybind11.h>
|
||||
#include <pybind11/stl.h>
|
||||
|
||||
namespace py = pybind11;
|
||||
using pd_type = double;
|
||||
|
||||
using namespace aare;
|
||||
|
||||
#pragma GCC diagnostic push
|
||||
#pragma GCC diagnostic ignored "-Wunused-parameter"
|
||||
|
||||
namespace aare {
|
||||
|
||||
// Binding for the OPT2 snapshot finder (pre-refactor pipeline: per-frame pinned
|
||||
// staging, round-robin streams with sync barriers, variable-length D2H). Kept
|
||||
// only for benchmarking the optimization arc; not part of the shipped API.
|
||||
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
|
||||
typename CoordType = uint16_t>
|
||||
void define_ClusterFinderCUDAOpt2(py::module &m, const std::string &typestr) {
|
||||
auto class_name = fmt::format("ClusterFinderCUDAOpt2_{}", typestr);
|
||||
|
||||
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
|
||||
using CF = ClusterFinderCUDAOpt2<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())
|
||||
// ctor: (image_size, n_sigma, capacity, n_streams) — capacity is the
|
||||
// per-stream device cluster buffer (upper bound on clusters/frame).
|
||||
.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") = 3000, 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.)")
|
||||
|
||||
.def("push_pedestal_frame",
|
||||
[](CF &self, ContigArr frame) {
|
||||
auto view = make_view_2d(frame);
|
||||
self.push_pedestal_frame(view);
|
||||
})
|
||||
|
||||
.def("clear_pedestal", &CF::clear_pedestal)
|
||||
|
||||
.def_property_readonly("pedestal",
|
||||
[](CF &self) {
|
||||
auto pd = new NDArray<pd_type, 2>{};
|
||||
*pd = self.pedestal();
|
||||
return return_image_data(pd);
|
||||
})
|
||||
|
||||
.def_property_readonly("noise",
|
||||
[](CF &self) {
|
||||
auto arr = new NDArray<pd_type, 2>{};
|
||||
*arr = self.noise();
|
||||
return return_image_data(arr);
|
||||
})
|
||||
|
||||
.def(
|
||||
"steal_clusters",
|
||||
[](CF &self, bool realloc_same_capacity) {
|
||||
return self.steal_clusters(realloc_same_capacity);
|
||||
},
|
||||
py::arg("realloc_same_capacity") = true)
|
||||
|
||||
.def(
|
||||
"find_clusters",
|
||||
[](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::call_guard<py::gil_scoped_release>())
|
||||
|
||||
.def(
|
||||
"find_clusters_batched",
|
||||
[](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,
|
||||
py::call_guard<py::gil_scoped_release>(),
|
||||
R"(Process a 3D array (n_frames, nrows, ncols) round-robin across
|
||||
n_streams. Returns a list of ClusterVector, one per input frame.)")
|
||||
|
||||
.def("avg_kernel_time_ms", &CF::avg_kernel_time_ms)
|
||||
.def("reset_timers", &CF::reset_timers);
|
||||
}
|
||||
|
||||
} // namespace aare
|
||||
|
||||
#pragma GCC diagnostic pop
|
||||
@@ -8,6 +8,7 @@
|
||||
#include "bind_Cluster.hpp"
|
||||
#include "bind_ClusterFinderCUDA.hpp"
|
||||
#include "bind_ClusterFinderCUDAGraph.hpp"
|
||||
#include "bind_ClusterFinderCUDAOpt2.hpp"
|
||||
#include "bind_ClusterVector.hpp"
|
||||
|
||||
#include <pybind11/pybind11.h>
|
||||
@@ -30,6 +31,10 @@ namespace py = pybind11;
|
||||
aare::define_ClusterFinderCUDAGraph<T, N, M, U>(m, "Cluster" #N \
|
||||
"x" #M #TYPE_CODE);
|
||||
|
||||
#define DEFINE_BINDINGS_CLUSTERFINDER_CUDA_OPT2(T, N, M, U, TYPE_CODE) \
|
||||
aare::define_ClusterFinderCUDAOpt2<T, N, M, U>(m, "Cluster" #N \
|
||||
"x" #M #TYPE_CODE);
|
||||
|
||||
PYBIND11_MODULE(_aare_cuda, m) {
|
||||
|
||||
// Types first — finders reference them in their signatures.
|
||||
@@ -83,8 +88,14 @@ PYBIND11_MODULE(_aare_cuda, m) {
|
||||
DEFINE_BINDINGS_CLUSTERFINDER_CUDA_GRAPH(int, 9, 9, uint16_t, i);
|
||||
DEFINE_BINDINGS_CLUSTERFINDER_CUDA_GRAPH(double, 9, 9, uint16_t, d);
|
||||
DEFINE_BINDINGS_CLUSTERFINDER_CUDA_GRAPH(float, 9, 9, uint16_t, f);
|
||||
|
||||
// OPT2 snapshot finder (benchmark only) — 3x3 is what the deck uses.
|
||||
DEFINE_BINDINGS_CLUSTERFINDER_CUDA_OPT2(int, 3, 3, uint16_t, i);
|
||||
DEFINE_BINDINGS_CLUSTERFINDER_CUDA_OPT2(double, 3, 3, uint16_t, d);
|
||||
DEFINE_BINDINGS_CLUSTERFINDER_CUDA_OPT2(float, 3, 3, uint16_t, f);
|
||||
}
|
||||
|
||||
#undef DEFINE_CUDA_CLUSTER_TYPES
|
||||
#undef DEFINE_BINDINGS_CLUSTERFINDER_CUDA
|
||||
#undef DEFINE_BINDINGS_CLUSTERFINDER_CUDA_GRAPH
|
||||
#undef DEFINE_BINDINGS_CLUSTERFINDER_CUDA_GRAPH
|
||||
#undef DEFINE_BINDINGS_CLUSTERFINDER_CUDA_OPT2
|
||||
File diff suppressed because one or more lines are too long
@@ -0,0 +1,32 @@
|
||||
# Minimal probe for nsys: train pedestal, run one batched pass, print summary.
|
||||
# Usage: nsys_kernel_probe.py [n_streams] [n_frames]
|
||||
import sys
|
||||
sys.path.append('/home/ferjao_k/aare/build')
|
||||
|
||||
from pathlib import Path
|
||||
import time
|
||||
from aare import File, ClusterFinderCUDA
|
||||
|
||||
n_streams = int(sys.argv[1]) if len(sys.argv) > 1 else 8
|
||||
N = int(sys.argv[2]) if len(sys.argv) > 2 else 2000
|
||||
|
||||
base = Path('/mnt/sls_det_storage/moench_data/2603_MaxIVBeamtime/2026032408/process/xrf/')
|
||||
f = File(base / 'Cu_factor_10_data_master_0.json')
|
||||
pd = File(base / 'Cu_factor_10_pedestal_master_0.json')
|
||||
|
||||
cf = ClusterFinderCUDA((f.rows, f.cols), (3, 3), n_sigma=5,
|
||||
max_clusters_per_frame=3000, n_streams=n_streams)
|
||||
for _ in range(1000):
|
||||
cf.push_pedestal_frame(pd.read_frame().copy())
|
||||
|
||||
data = f.read_n(N)
|
||||
cf.register_input_buffer(data)
|
||||
|
||||
t0 = time.perf_counter()
|
||||
res = cf.find_clusters_batched(data, first_frame=0)
|
||||
t = time.perf_counter() - t0
|
||||
|
||||
cf.unregister_input_buffer()
|
||||
n = sum(cv.size for cv in res)
|
||||
print(f'n_streams={n_streams} N={N} wall={t:.3f}s ({N/t:.0f} FPS) '
|
||||
f'clusters/frame={n/N:.2f} event kernel_ms={cf.avg_kernel_time_ms():.3f}')
|
||||
Reference in New Issue
Block a user