perf: add nsys probes for PixelHistogramCUDA

perf/hist/ mirrors perf/cf/: a driver that runs each config under nsys,
the traced process, settings and a README. The fill (H2D + fill_bin_major
per frame) and the readout (bin_to_pixel_major + one D2H per values())
are reported separately. Only one chunk of frames is held in host memory
and re-queued until the requested frame count.

gpu_span.analyze() takes an optional engines subset, so the one-off
readout D2H doesn't count towards the fill window or roofline. Default
output is unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

python/tests/perf/README.md

@@ -7,13 +7,14 @@ with its own driver, traced process, README and `results/`:
 | folder | class |
 |---|---|
 | [`cf/`](cf/README.md) | `ClusterFinderCUDA` |
+| [`hist/`](hist/README.md) | `PixelHistogramCUDA` |
 
 Shared by all of them:
 
 | file | role |
 |---|---|
 | `common.py` | paths (`AARE_REPO`, `AARE_BUILD`, `AARE_PERF_DATA`, `NSYS`), guards (fresh build, idle GPU), `env.json` stamp, results folder |
-| `gpu_span.py` | reads an nsys SQLite export: busy time, duty cycle and roofline per engine, optionally for one kernel name only |
+| `gpu_span.py` | reads an nsys SQLite export: busy time, duty cycle and roofline per engine, optionally for one kernel name and a subset of engines only |
 
 Requirements (a CUDA build of the Python extension, Nsight Systems, the
 dataset, an idle GPU) are listed in each folder's README.

python/tests/perf/gpu_span.py

@@ -10,6 +10,9 @@ pedestal-upload copies that precede any kernel are excluded.
 
 `kernel_name` (optional) keeps only kernels whose short name contains it, for
 a probe whose per-frame kernel shares the trace with one-off kernels.
+`engines` (optional, Python only) limits the engines that are reported and
+that count for the window and the roofline: a probe whose only D2H is a
+one-off readout at the end passes ("kernel", "H2D").
 
 Roofline = tallest engine's busy time per frame. PCIe is full-duplex and the
 engines overlap, so the floor is the max, never the sum.
@@ -36,16 +39,19 @@ def union(rows):
     return busy + ce - cs
 
 
-def analyze(sqlite_path, n_frames, kernel_name=None) -> dict:
+ENGINES = ("kernel", "H2D", "D2H")
+
+
+def analyze(sqlite_path, n_frames, kernel_name=None, engines=ENGINES) -> dict:
     db = sqlite3.connect(str(sqlite_path))
     kernels = [(s, e) for s, e, name in db.execute(KERNELS)
                if kernel_name is None or kernel_name in name]
     if not kernels:
         raise RuntimeError(f"no kernel matching {kernel_name!r} in {sqlite_path}")
     ops = {"kernel": kernels}
-    ops.update({k: db.execute(q).fetchall() for k, q in Q.items()})
+    ops.update({k: db.execute(q).fetchall() for k, q in Q.items() if k in engines})
     lo = ops["kernel"][0][0]
-    hi = max([ops["kernel"][-1][1]] + [r[1] for r in ops["D2H"][-1:]])
+    hi = max([ops["kernel"][-1][1]] + [r[1] for r in ops.get("D2H", [])[-1:]])
     span = hi - lo
 
     out = {"n_frames": n_frames, "window_ms": span / 1e6,
@@ -59,7 +65,7 @@ def analyze(sqlite_path, n_frames, kernel_name=None) -> dict:
         out[f"{k}_overlap"] = tot / busy if busy else 0.0
         out[f"{k}_duty_pct"] = 100 * busy / span
         out[f"{k}_us_per_frame"] = busy / 1e3 / n_frames
-    tallest = max(("kernel", "H2D", "D2H"), key=lambda k: out[f"{k}_us_per_frame"])
+    tallest = max(ops, key=lambda k: out[f"{k}_us_per_frame"])
     out["bottleneck"] = tallest
     out["roofline_us_per_frame"] = out[f"{tallest}_us_per_frame"]
     out["roofline_fps"] = 1e6 / out["roofline_us_per_frame"]

python/tests/perf/hist/README.md

@@ -0,0 +1,89 @@
+# PixelHistogramCUDA nsys probes
+
+GPU time of `PixelHistogramCUDA`, traced with Nsight Systems on real MOENCH
+frames. Two costs, reported separately:
+
+- **Fill, per frame**: the H2D copy of the frame and the `fill_bin_major`
+  kernel. They overlap (two frame slots, a copy stream and a compute stream),
+  so throughput is bound by the taller of the two, not their sum.
+- **Readout, once per `values()` call**: the `bin_to_pixel_major` kernel, then
+  one D2H of `rows × cols × n_bins` counts. It doesn't depend on the number of
+  frames.
+
+These are *not* throughput numbers. nsys inflates wall clock by tracing every
+CUDA call. Use them for "what does the kernel cost" and "which engine binds",
+and measure sustained FPS separately without a profiler.
+
+## Files
+
+| file | role |
+|---|---|
+| `run_probes.py` | driver: runs each config under nsys, exports SQLite, writes `probes.csv` |
+| `nsys_hist_probe.py` | the traced process: fill one histogram with N frames, read it out, exit |
+| `hist_common.py` | histogram settings: headers for the rebuild guard, kernel names, results folder |
+| `../gpu_span.py`, `../common.py` | shared with the other probes, see [`../README.md`](../README.md) |
+
+## Requirements
+
+As for the cluster-finder probes ([`../cf/README.md`](../cf/README.md)): a
+fresh CUDA build of the extension, a Python of the same minor version with
+numpy and matplotlib, Nsight Systems, the MOENCH dataset, an idle GPU. Also:
+
+- **Host memory**: only `--chunk` frames are loaded (2000 by default, 640 MB as
+  uint16, twice that as float32), plus the output histogram. The chunk is
+  pinned once and queued repeatedly until `--frames` frames have gone through:
+  the histogram doesn't depend on frame order, and every pass is a real H2D
+  copy.
+- **GPU memory**: about `2 × rows × cols × n_bins × 4` bytes, plus two frame
+  batches. The accumulator stays allocated and `values()` allocates a
+  temporary of the same size for the transpose. 1000 bins of 400 × 400 pixels
+  needs about 1.4 GB.
+
+## Run
+
+```bash
+cd python/tests/perf/hist
+python run_probes.py                # full sweep: 4 configs, 20k frames each
+python run_probes.py --frames 2000  # smoke test only; clocks do not ramp
+python run_probes.py --only u16_b1000_bs64 --tag mychange
+```
+
+| label | frames | bins | batch | what it shows |
+|---|---|---|---|---|
+| `u16_b1000_bs64` | uint16 | 1000 | 64 | reference: raw frames, default batch size |
+| `f32_b1000_bs64` | float32 | 1000 | 64 | twice the H2D bytes per frame |
+| `u16_b100_bs64` | uint16 | 100 | 64 | 64 MB histogram, about the size of L2: increments mostly hit cache |
+| `u16_b1000_bs512` | uint16 | 1000 | 512 | fewer, longer copies and kernels |
+
+The range is taken from the data (0.1 to 99.9 percentile of the first frames),
+so nearly every value is histogrammed; the probe prints the fraction in range.
+Counts are uint32.
+
+Output goes to `results/<date>[_tag]/`: `env.json`, one `.nsys-rep` +
+`.sqlite` per config, and `probes.csv`. Re-running one config replaces its own
+row and keeps the others. `probes.csv` and `env.json` can be committed with a
+result; the traces are gitignored.
+
+To inspect one trace by hand:
+
+```bash
+nsys stats --report cuda_gpu_sum results/<dir>/probe_u16_b1000_bs64.nsys-rep
+python ../gpu_span.py results/<dir>/probe_u16_b1000_bs64.sqlite 20000 fill_bin_major
+```
+
+Without the kernel name, `gpu_span.py` also counts the readout kernel and the
+readout D2H as if they were per-frame work.
+
+## Reading the numbers
+
+- The fill window starts at the first fill kernel, so the copy of the very
+  first batch is left out: one batch in `frames / batch`, under 0.5 % at the
+  defaults.
+- A kernel below the H2D bar is invisible end to end; kernel work only pays
+  once the kernel is the taller bar.
+- Readout cost grows with `rows × cols × n_bins` and is paid per `values()`
+  call: amortise it by reading out rarely, and pass a registered `out` array
+  so the D2H runs at full pinned bandwidth.
+
+No reference results yet: the first sweep on an idle GPU will be the
+baseline.

python/tests/perf/hist/hist_common.py

@@ -0,0 +1,27 @@
+"""PixelHistogramCUDA-specific settings for the probes in this folder.
+
+Importing this module also makes ../common.py and ../gpu_span.py importable.
+"""
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+HERE = Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE.parent))
+import common  # noqa: E402
+
+# Rebuild-guard: the extension must be newer than all of these.
+HEADERS = ["include/aare/hist/PixelHistogramCUDA.hpp",
+           "include/aare/hist/pixelhistogram_kernel.cuh",
+           "include/aare/RAII_device.hpp"]
+
+# Kernel short names in the trace. The fill kernel runs once per batch; the
+# readout kernel once per values() call.
+FILL_KERNEL = "fill_bin_major"
+READOUT_KERNEL = "bin_to_pixel_major"
+
+
+def results_dir(tag: str = "") -> Path:
+    """hist/results/<date>[_tag]/"""
+    return common.results_dir(HERE, "", tag)

python/tests/perf/hist/nsys_hist_probe.py

@@ -0,0 +1,70 @@
+"""The process nsys traces: fill one histogram with N frames, read it out, exit.
+
+    nsys_hist_probe.py <dtype> <n_frames> <n_bins> <batch> <chunk>
+
+dtype is u16, f32 or f64 (the frame dtype copied to the GPU). Run by
+run_probes.py; run directly only under `nsys profile --trace=cuda`. The wall
+times printed are profiler-inflated and are not throughput numbers.
+
+Only `chunk` frames are read and held in host memory. They are registered once
+and passed to fill_histogram() repeatedly until n_frames have been queued: the
+histogram doesn't depend on frame order, and every pass is a real H2D copy.
+The chunk is rounded down to a multiple of the batch size so that no pass ends
+in a partial batch.
+"""
+import math
+import sys
+import time
+from pathlib import Path
+
+import numpy as np
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+import hist_common  # noqa: E402,F401  (puts ../ and the build dir on sys.path)
+from hist_common import common  # noqa: E402
+
+from aare import File, PixelHistogramCUDA  # noqa: E402
+
+DTYPES = {"u16": np.uint16, "f32": np.float32, "f64": np.float64}
+
+dtype = DTYPES[sys.argv[1]]
+N, n_bins, batch, chunk = (int(a) for a in sys.argv[2:6])
+chunk = max(batch, min(chunk, N) // batch * batch)
+
+f = File(common.DATA_FILE)
+data = f.read_n(chunk)
+# Range from the data itself, so that nearly every value is histogrammed.
+lo, hi = np.percentile(data[:min(chunk, 200)], [0.1, 99.9])
+xmin, xmax = float(math.floor(lo)), float(math.ceil(hi) + 1)
+data = np.ascontiguousarray(data.astype(dtype, copy=False))
+rows, cols = data.shape[1:]
+
+hist = PixelHistogramCUDA(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin,
+                          xmax=xmax, dtype=data.dtype, storage_dtype=np.uint32,
+                          batch_size=batch)
+out = np.empty((rows, cols, n_bins), dtype=np.uint32)
+hist.register_input_buffer(data)
+hist.register_output_buffer(out)
+
+t0 = time.perf_counter()
+queued = 0
+while queued < N:
+    n = min(chunk, N - queued)
+    hist.fill_histogram(data[:n])
+    queued += n
+hist.flush()
+t_fill = time.perf_counter() - t0
+
+t0 = time.perf_counter()
+hist.values(out=out)
+t_read = time.perf_counter() - t0
+hist.unregister_input_buffer()
+hist.unregister_output_buffer()
+
+in_range = out.sum(dtype=np.uint64) / (N * rows * cols)
+print(f"dtype={data.dtype}  N={N}  n_bins={n_bins}  batch={batch}  chunk={chunk}  "
+      f"range=[{xmin:g}, {xmax:g})")
+print(f"  H2D/frame={rows * cols * data.itemsize:,} B   "
+      f"readout={out.nbytes:,} B   in range={100 * in_range:.2f}%")
+print(f"  wall fill={t_fill:.3f}s ({N / t_fill:.0f} FPS)  readout={t_read * 1e3:.1f} ms  "
+      f"(profiler-inflated)")

python/tests/perf/hist/run_probes.py

@@ -0,0 +1,150 @@
+#!/usr/bin/env python3
+"""nsys probe sweep for PixelHistogramCUDA: fill cost per frame, readout cost.
+
+    python3.11 run_probes.py                       # full sweep, 20k frames
+    python3.11 run_probes.py --frames 2000         # quick check (clocks may not ramp)
+    python3.11 run_probes.py --only u16_b1000_bs64 --tag mytag
+
+Per config: an .nsys-rep + .sqlite in results/<date>[_tag]/ and one row in
+probes.csv (re-running a config replaces its own row only).
+
+Two separate costs:
+  fill     per frame: H2D of the frame and fill_bin_major, overlapping. The
+           roofline is the taller of the two.
+  readout  once per values() call: bin_to_pixel_major, then one D2H of
+           rows * cols * n_bins counts. Independent of the number of frames.
+"""
+from __future__ import annotations
+
+import argparse
+import csv
+import sqlite3
+import subprocess
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+import hist_common  # noqa: E402  (puts ../ and the build dir on sys.path)
+from hist_common import common  # noqa: E402
+import gpu_span  # noqa: E402
+
+HERE = Path(__file__).resolve().parent
+
+# (dtype, n_bins, batch, label)
+#   u16_b1000_bs64   reference: raw frames, default batch size
+#   f32_b1000_bs64   twice the H2D bytes per frame
+#   u16_b100_bs64    histogram of 160k pixels x 100 bins x 4 B = 64 MB, about
+#                    the size of L2: the fill kernel's scattered increments
+#                    mostly hit cache
+#   u16_b1000_bs512  larger batches: fewer, longer copies and kernels
+CONFIGS = [
+    ("u16", 1000, 64, "u16_b1000_bs64"),
+    ("f32", 1000, 64, "f32_b1000_bs64"),
+    ("u16", 100, 64, "u16_b100_bs64"),
+    ("u16", 1000, 512, "u16_b1000_bs512"),
+]
+
+
+def readout(sqlite_path) -> dict:
+    """Total time of the readout kernel and of all D2H copies."""
+    db = sqlite3.connect(str(sqlite_path))
+    k = [e - s for s, e, name in db.execute(gpu_span.KERNELS)
+         if hist_common.READOUT_KERNEL in name]
+    d2h_ns, d2h_bytes = db.execute(
+        "SELECT COALESCE(SUM(end-start), 0), COALESCE(SUM(bytes), 0) "
+        "FROM CUPTI_ACTIVITY_KIND_MEMCPY WHERE copyKind=2").fetchone()
+    return {"readout_kernel_n": len(k), "readout_kernel_ms": sum(k) / 1e6,
+            "readout_D2H_ms": d2h_ns / 1e6, "readout_D2H_MB": d2h_bytes / 1e6,
+            "readout_D2H_GBps": d2h_bytes / d2h_ns if d2h_ns else 0.0}
+
+
+def run_one(dtype, n_bins, batch, label, n_frames, chunk, outdir) -> dict | None:
+    rep = outdir / f"probe_{label}"
+    print(f"\n--- {label}: {dtype} frames, {n_bins} bins, batch={batch}, N={n_frames} ---")
+
+    prof = [common.NSYS, "profile", "--trace=cuda", "--sample=none",
+            "--cpuctxsw=none", "--force-overwrite=true", "-o", str(rep),
+            sys.executable, str(HERE / "nsys_hist_probe.py"),
+            dtype, str(n_frames), str(n_bins), str(batch), str(chunk)]
+    p = subprocess.run(prof, capture_output=True, text=True)
+    for line in p.stdout.splitlines():
+        if line.strip().startswith(("dtype=", "H2D/frame", "wall")):
+            print("   ", line.strip())
+    # nsys writes a report even when the traced process crashes, and returns
+    # the child's exit code: check it, or a dead probe surfaces later as a
+    # missing CUPTI table in the export.
+    if p.returncode != 0 or not rep.with_suffix(".nsys-rep").exists():
+        print(f"    FAILED (exit {p.returncode}): {p.stderr.strip()[-400:]}")
+        return None
+
+    subprocess.run([common.NSYS, "stats", "--force-export=true", "--report",
+                    "cuda_gpu_sum", str(rep.with_suffix(".nsys-rep"))],
+                   capture_output=True, text=True)
+    sq = rep.with_suffix(".sqlite")
+    if not sq.exists():
+        print("    FAILED: no sqlite export")
+        return None
+
+    r = gpu_span.analyze(sq, n_frames, hist_common.FILL_KERNEL,
+                         engines=("kernel", "H2D"))
+    r.update(readout(sq))
+    r.update(label=label, dtype=dtype, n_bins=n_bins, batch=batch)
+    print(f"    fill: kernel {r['kernel_us_per_frame']:5.2f} us (duty {r['kernel_duty_pct']:4.1f}%)  "
+          f"H2D {r['H2D_us_per_frame']:5.2f} ({r['H2D_duty_pct']:4.1f}%)  "
+          f"-> {r['bottleneck']}-bound at {r['roofline_us_per_frame']:.2f} us/frame "
+          f"= {r['roofline_fps']:,.0f} FPS")
+    print(f"    readout: kernel {r['readout_kernel_ms']:.2f} ms  D2H {r['readout_D2H_ms']:.2f} ms "
+          f"({r['readout_D2H_MB']:,.0f} MB at {r['readout_D2H_GBps']:.1f} GB/s)")
+    return r
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description=__doc__,
+                                 formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--frames", type=int, default=20_000)
+    ap.add_argument("--chunk", type=int, default=2000,
+                    help="frames held in host memory and re-queued until --frames")
+    ap.add_argument("--tag", default="")
+    ap.add_argument("--only", nargs="+", default=None, help="subset of labels")
+    args = ap.parse_args()
+
+    common.assert_build_fresh(hist_common.HEADERS)
+    common.assert_idle_gpu()
+    env = common.capture_env()
+    outdir = hist_common.results_dir(args.tag)
+    common.write_env(outdir / "env.json", env)
+    print(f"build: git={env['git_rev']}{' (dirty)' if env['git_dirty'] else ''}")
+    print(f"out:   {outdir}")
+
+    rows = []
+    for dtype, n_bins, batch, label in CONFIGS:
+        if args.only and label not in args.only:
+            continue
+        r = run_one(dtype, n_bins, batch, label, args.frames, args.chunk, outdir)
+        if r:
+            rows.append(r)
+
+    if not rows:
+        return 1
+    out = outdir / "probes.csv"
+    prior = list(csv.DictReader(out.open())) if out.exists() else []
+    fresh = {r["label"] for r in rows}
+    merged = [r for r in prior if r["label"] not in fresh] + rows
+    with out.open("w", newline="") as fh:
+        w = csv.DictWriter(fh, fieldnames=list(rows[0]))
+        w.writeheader()
+        w.writerows(merged)
+
+    print("\n=== fill (us/frame) and readout (ms, once per values()) ===")
+    print(f"{'config':<18} {'kernel':>8} {'H2D':>8}  {'bound':<6} {'FPS':>10}  "
+          f"{'rd kernel':>9} {'rd D2H':>8}")
+    for r in rows:
+        print(f"{r['label']:<18} {r['kernel_us_per_frame']:8.2f} {r['H2D_us_per_frame']:8.2f}  "
+              f"{r['bottleneck']:<6} {r['roofline_fps']:10,.0f}  "
+              f"{r['readout_kernel_ms']:9.2f} {r['readout_D2H_ms']:8.2f}")
+    print(f"\n-> {out}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())