cuda: add a fused candidate-scoring kernel for CLM and Kev

docs/benchmarks/scoring-fusion/README.md

@@ -0,0 +1,116 @@
+# What fusing the candidate scoring is worth
+
+`src/backends/cuda/scoring/` argues that a cosine needs each candidate's norm and
+its dot with the query, that both read the same elements, and that accumulating
+them in one pass over `candidates` is therefore cheaper than composing a
+normalisation pass with a scoring pass. This measures it, and the measurement
+changed the kernel.
+
+## Protocol and controls
+
+The hypothesis was that the fused kernel beats an unfused equivalent that reads
+the candidate matrix twice. Stop conditions were any failed correctness check, a
+missing library, or no GPU; a correctness failure aborts the run rather than
+producing a number, because a kernel that computes the wrong thing can be
+arbitrarily fast.
+
+- One shared RTX 4090, 1 of 8 on the host, otherwise idle (0 % utilisation,
+  30 °C) for the recorded run. The host is shared and its other GPUs carry
+  unrelated work, so the **ratios** are the result here, not the absolute
+  microseconds. The fused and unfused measurements are taken in the same
+  process, on the same buffers, and were repeated with the order swapped to rule
+  out a clock-ramp bias: the ratios moved by at most 0.01.
+- Both paths are checked against `reference.py` first, per shape, at a `max_abs`
+  of 1e-4. Both passed before anything was timed.
+- Per shape: 30 warmup calls, then 300 calls timed individually with CUDA events;
+  the median is reported.
+- Every timing in a row is taken on **one created stream**, including the graph
+  replay. Measuring on the legacy default stream instead adds its synchronisation
+  semantics to both and makes the columns incomparable — it cost the fused kernel
+  about 5 µs at K = 255, which is larger than the effect being measured.
+- `bench.cu` is the baseline. It is `candidate_scoring.cu`'s own block with the
+  row norm loaded from a staged array instead of accumulated beside the dot, plus
+  a separate norm pass. The norm pass runs only when the similarity is a cosine,
+  because a plain dot has no norm to compute and running it anyway would slow the
+  baseline for a reason the fused kernel never removed.
+- It is a **conservative** baseline: a real unfused pipeline built from a GEMM
+  would also round-trip the similarity matrix through global memory, which this
+  does not. The differences below are a floor.
+
+```sh
+src/backends/cuda/scoring/build.sh ./out 89
+python3 src/backends/cuda/scoring/bench.py --library ./out/libscoring.so --iterations 300 \
+    --out docs/benchmarks/scoring-fusion/results.json
+```
+
+## The finding: the kernel was occupancy-bound, not bandwidth-bound
+
+At the block size the kernel was written with, 256 threads, fusion **lost** at
+every shape that mattered:
+
+| shape | similarity | K | D | questions | fused | unfused | ratio |
+| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: |
+| `clm-cosine-5x512` | cosine | 5 | 512 | 1 | 9.2 µs | 11.3 µs | 1.23× |
+| `clm-cosine-64x512` | cosine | 64 | 512 | 1 | 19.4 µs | 19.3 µs | 1.00× |
+| `clm-cosine-255x512` | cosine | 255 | 512 | 1 | 52.9 µs | 43.2 µs | **0.82×** |
+| `clm-cosine-batch8-255x512` | cosine | 255 | 512 | 8 | 53.1 µs | 44.0 µs | **0.83×** |
+| `kev-dot-255x2560` | dot | 255 | 2560 | 1 | 167.8 µs | 96.1 µs | **0.57×** |
+
+Two things said why. The achieved bandwidth was 1 to 190 GB/s against a device
+peak near 1000, so nothing was bandwidth-bound and there was no read to save. And
+a single question and a batch of eight took **the same 53 µs** — eight times the
+work in the same wall clock, which is only possible if each question is running
+on its own SM. One thread block per question means a one-question request, the
+common serving case, occupies **one of the 128 SMs**, and the norm accumulation
+sat on that block's critical path instead of being spread across the device.
+
+Raising the block to 1024 threads — the maximum on sm_89, and one line — moves
+the work onto 32 warps instead of 8 and reverses every result:
+
+| shape | similarity | K | D | questions | fused | unfused | ratio |
+| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: |
+| `clm-cosine-5x512` | cosine | 5 | 512 | 1 | 9.2 µs | 11.3 µs | 1.22× |
+| `clm-cosine-64x512` | cosine | 64 | 512 | 1 | 12.1 µs | 18.3 µs | 1.51× |
+| `clm-cosine-255x512` | cosine | 255 | 512 | 1 | 23.6 µs | 42.2 µs | **1.79×** |
+| `clm-cosine-batch8-255x512` | cosine | 255 | 512 | 8 | 23.8 µs | 43.8 µs | **1.84×** |
+| `kev-dot-255x2560` | dot | 255 | 2560 | 1 | 55.3 µs | 100.3 µs | **1.81×** |
+
+Kev is the control. It is a plain dot, so there is no norm for the fusion to
+remove from the traffic, and its gain is entirely the second kernel and the
+global similarity round-trip it does not need. That it lands at the same ~1.8× as
+the cosine shapes is the useful part of the result.
+
+Correctness is unchanged by the block size: the parity harness reports the same
+`1.835e-07` worst case, and `compute-sanitizer --tool memcheck` reports
+`ERROR SUMMARY: 0 errors`.
+
+## CUDA Graph replay
+
+The ABI is written for capture — it queues on the caller's stream and never
+synchronizes — so the same launch can be replayed from a graph. It is timed here
+because the smallest shapes are launch-bound, and that is the part a graph
+removes:
+
+| shape | direct | graph replay | gain |
+| --- | ---: | ---: | ---: |
+| `clm-cosine-5x512` | 9.2 µs | 6.1 µs | **1.49x** |
+| `clm-cosine-64x512` | 12.1 µs | 9.2 µs | 1.31x |
+| `clm-cosine-255x512` | 24.2 µs | 20.7 µs | 1.17x |
+| `clm-cosine-batch8-255x512` | 24.4 µs | 20.5 µs | 1.19x |
+| `kev-dot-255x2560` | 52.4 µs | 49.3 µs | 1.06x |
+
+The saving is a flat ~3 µs, which is what a launch costs; it is 1.49x of a 9 µs
+call and 1.06x of a 52 µs one. That is the expected shape and it explains the
+floor at K = 5, but it also means the graph gain **shrinks as the shapes grow** —
+so it does not change which configuration wins at any shape measured here. A
+graph and the fused kernel are complementary rather than alternatives, and the
+two effects multiply.
+
+## Limits
+
+- One GPU, one architecture: sm_89 only. The library builds for whatever
+  `build.sh` is given, but nothing else has been run.
+- The tiny shapes are launch-bound, not compute-bound — `clm-cosine-5x512` is 9 µs
+  either way. The interesting range starts around K = 64.
+- No end-to-end number. This is the kernel, with the projections left to the
+  caller as the design intends.

docs/benchmarks/scoring-fusion/results.json

@@ -0,0 +1,111 @@
+{
+  "measured_at": "2026-10-03T13:15:08Z",
+  "provenance": {
+    "gpu": "NVIDIA GeForce RTX 4090",
+    "driver": "580.76.05",
+    "nvcc": "Cuda compilation tools, release 13.0, V13.0.88",
+    "library": "./out/libscoring.so",
+    "baseline_build": "nvcc -O3 -std=c++17 -gencode arch=compute_89,code=sm_89 -shared -Xcompiler -fPIC -Xcompiler -Wall,-Wextra bench.cu -lcudart -o libscoring_bench.so",
+    "iterations": 300,
+    "arch": "89",
+    "kernel_definition": "THREADS = 1024 in candidate_scoring.cu"
+  },
+  "results": [
+    {
+      "shape": "clm-cosine-5x512",
+      "questions": 1,
+      "K": 5,
+      "D": 512,
+      "cosine": true,
+      "fused_ms": 0.009151999838650227,
+      "unfused_ms": 0.011455999687314034,
+      "speedup": 1.2517482396507134,
+      "fused_gbs": 1.3448426810522358,
+      "unfused_gbs": 1.9682263107050233,
+      "fused_bytes": 12308,
+      "unfused_bytes": 22548,
+      "fused_max_abs": 6.932993346087102e-09,
+      "unfused_max_abs": 6.932993346087102e-09,
+      "direct_ms": 0.009184000082314014,
+      "graph_replay_ms": 0.006144000217318535,
+      "graph_gain": 1.4947916271920716
+    },
+    {
+      "shape": "clm-cosine-64x512",
+      "questions": 1,
+      "K": 64,
+      "D": 512,
+      "cosine": true,
+      "fused_ms": 0.012160000391304493,
+      "unfused_ms": 0.018432000651955605,
+      "speedup": 1.515789478521412,
+      "fused_gbs": 10.968420699671686,
+      "unfused_gbs": 14.34722171474872,
+      "fused_bytes": 133376,
+      "unfused_bytes": 264448,
+      "fused_max_abs": 1.2430758852821633e-09,
+      "unfused_max_abs": 1.2430758852821633e-09,
+      "direct_ms": 0.01206399966031313,
+      "graph_replay_ms": 0.009184000082314014,
+      "graph_gain": 1.3135888014140202
+    },
+    {
+      "shape": "clm-cosine-255x512",
+      "questions": 1,
+      "K": 255,
+      "D": 512,
+      "cosine": true,
+      "fused_ms": 0.02425600029528141,
+      "unfused_ms": 0.043007999658584595,
+      "speedup": 1.7730870355798547,
+      "fused_gbs": 21.656826913140733,
+      "unfused_gbs": 24.357050044546877,
+      "fused_bytes": 525308,
+      "unfused_bytes": 1047548,
+      "fused_max_abs": 1.2364590072297399e-09,
+      "unfused_max_abs": 1.2364590072297399e-09,
+      "direct_ms": 0.024240000173449516,
+      "graph_replay_ms": 0.020735999569296837,
+      "graph_gain": 1.168981514126811
+    },
+    {
+      "shape": "clm-cosine-batch8-255x512",
+      "questions": 8,
+      "K": 255,
+      "D": 512,
+      "cosine": true,
+      "fused_ms": 0.024383999407291412,
+      "unfused_ms": 0.04398399963974953,
+      "speedup": 1.803805803349767,
+      "fused_gbs": 172.34514854619624,
+      "unfused_gbs": 190.53255885411613,
+      "fused_bytes": 4202464,
+      "unfused_bytes": 8380384,
+      "fused_max_abs": 2.1920026549090976e-09,
+      "unfused_max_abs": 2.1920026549090976e-09,
+      "direct_ms": 0.024383999407291412,
+      "graph_replay_ms": 0.020479999482631683,
+      "graph_gain": 1.1906250011368684
+    },
+    {
+      "shape": "kev-dot-255x2560",
+      "questions": 1,
+      "K": 255,
+      "D": 2560,
+      "cosine": false,
+      "fused_ms": 0.05427199974656105,
+      "unfused_ms": 0.09728000313043594,
+      "speedup": 1.7924528962395585,
+      "fused_gbs": 48.320681239798475,
+      "unfused_gbs": 26.9578527509269,
+      "fused_bytes": 2622460,
+      "unfused_bytes": 2622460,
+      "fused_max_abs": 7.07020300361183e-08,
+      "unfused_max_abs": 7.07020300361183e-08,
+      "direct_ms": 0.05238400027155876,
+      "graph_replay_ms": 0.04934399947524071,
+      "graph_gain": 1.0616083177011915
+    }
+  ],
+  "note": "Ratios, not absolute microseconds: the host is shared (1 of 8 GPUs) and its other GPUs carry unrelated work. Every timing in a row is taken on one created stream. The same harness at THREADS = 256 gave 0.82x, 0.83x and 0.57x on clm-cosine-255x512, clm-cosine-batch8-255x512 and kev-dot-255x2560; see README.md."
+}

src/backends/cuda/scoring/README.md

@@ -0,0 +1,159 @@
+# Fused candidate scoring
+
+One CUDA kernel that turns a query vector and a matrix of candidate vectors into
+a probability distribution, without materialising the similarities in between.
+
+This is the readout both models that need it share the shape of:
+
+| Model | Similarity | Notes |
+| --- | --- | --- |
+| Kev ([#27](https://github.com/ThinkFlowLab/system1-omni/issues/27)) | `scale * <k_proj(c_i), q_proj(s)>` | candidates are hidden states, one question per prefill row |
+| CLM ([#9](https://github.com/ThinkFlowLab/system1-omni/issues/9)) | `exp(logit_scale) * cos(state_head(s), action_head(c))` | cosine, so `normalize = 1` |
+
+They differ in the similarity and in whether a projection happens first; what
+they share is the primitive. `#9` describes CLM's answer as "the softmax over the
+question's candidates" of `exp(logit_scale) * cos(state_head(s), action_head(c))`,
+and `#27` describes Kev's as `scale * (k(h_opts) @ q(h_decide))` with
+`scale = 1/sqrt(256)`, softmaxed over that question's candidates. Neither is in
+`#19`'s `gemm.cu`, which is a bf16 cuBLASLt GEMM.
+
+## What is fused, precisely
+
+Not merely "fewer launches". When the similarity is a cosine — CLM's case, not
+Kev's, which is a plain dot — each candidate's norm and its dot with the query
+need **the same elements**, so both accumulate in a single read of `C`. The
+unfused version composes a normalisation pass with a scoring pass and reads `C`
+twice; at CLM's widest, 255 candidates of 512 floats, that is 522 KB per question
+read twice and once here. Beyond that:
+
+- One thread block per question, 1024 threads; the stable softmax runs over
+  shared memory, so there is no second kernel, no atomics, and no global round
+  trip for the similarities. The block size is not cosmetic — see
+  [the benchmark](../../../docs/benchmarks/scoring-fusion/README.md), where 256
+  threads made the fusion a **loss** at every shape that mattered and 1024
+  threads made it a 1.8x win.
+- `K ≤ 255` (the serving API's limit), which is what makes the in-block softmax
+  possible at all — the whole similarity vector fits in shared memory.
+- A batch is a grid over questions, not a loop of launches.
+
+The query norm is computed once per block rather than once per candidate.
+
+## Files
+
+| File | What it is | Verified |
+| --- | --- | --- |
+| `reference.py` | The oracle: float64, deliberately slow and obvious | **yes**, 17 tests |
+| `kernel_simulation.py` | The kernel's algorithm in numpy: float32, warp tree order, same norm flooring | **yes** |
+| `candidate_scoring.cu` | The kernel and its C ABI | **yes**, compiled and checked against the reference on sm_89 |
+| `bench.cu` | The unfused baseline the kernel is measured against | **yes**, timed on sm_89 |
+| `bench.py` | Drives both and writes the report | **yes** |
+
+## Verification state
+
+**Verified on real hardware**, on two driver stacks. RTX 4090 (sm_89), CUDA 13.0
+(V13.0.88), compiled with `./build.sh ./out 89`:
+
+| driver | worst absolute difference (tolerance 1e-4) |
+| --- | --- |
+| 595.58.03 | 1.835e-07 |
+| 580.76.05 | 1.835e-07 |
+
+`gpu_parity.py` checks three things, and all three run:
+
+- **13 fixed cases**, including the degenerate ones (single candidate, identical
+  candidates, a zero candidate) and the one whose raw logits overflow a naive
+  softmax.
+- **The batch path**, as a grid over questions: 1, 8, 4 and 3 questions with `K`
+  up to 255 and `D` up to 256. Each row is compared with the reference and checked
+  to sum to 1. The single-question call goes through `cs_score_candidates` and the
+  rest through `cs_score_candidates_batch`, so the two entry points have to agree.
+- **Rejected arguments**: `K` above 255, `K` of zero, `D` of zero, a zero scale, a
+  negative temperature and a null query must each return an error rather than a
+  scored answer. `K` above 255 is the one that matters: the shared array is sized
+  `MAX_K`, so scoring the first 255 of a larger set would be a confident wrong
+  answer instead of a failure.
+
+### Two defects the first GPU run found
+
+Both were in how the kernel was called or reduced; neither was visible without a
+GPU, and neither was caught by the simulation below.
+
+1. **The parity harness passed host pointers.** A numpy buffer is host memory, so
+   the kernel dereferenced it as device memory — `cudaErrorIllegalAddress`, not a
+   wrong answer. `compute-sanitizer` located it as an invalid global read with
+   consecutive lanes on consecutive addresses, which is the coalesced load
+   pattern. Fixed by allocating with `cudaMalloc`, copying in, synchronising and
+   copying out.
+2. **The query norm was counted `warps` times.** The inner loop
+   `d = lane; d < D; d += WARP` already covers every element of `D` across one
+   warp's 32 lanes, but every warp computed the same partial and the partials
+   were then summed across warps — so the norm came out `sqrt(8)` too large on a
+   256-thread block. Every unnormalized case sat at float32 rounding while the
+   normalize cases were off by ~1e-2, which is what pointed at it. Fixed by
+   having warp 0 compute it alone.
+
+That second one is worth recording precisely: **the simulation did not catch it.**
+`kernel_simulation.py` computed the query norm correctly, because it was written
+from the intent rather than transcribed from the `.cu`. A simulation is only as
+good as its fidelity, and the hardware is what settles it.
+
+## What the reference and simulation established before hardware
+
+1. **The reference is right.** All cases sum to 1; a single candidate gets 1.0;
+   identical candidates get exactly `1/K`; a zero candidate does not produce NaN;
+   logits large enough to overflow a naive softmax give a finite distribution.
+2. **The algorithm agrees with the reference.** `kernel_simulation.py` reproduces
+   the intended arithmetic — float32 accumulation, the max-subtracted softmax,
+   the norm floored as `max(sqrt(sum), eps)` and **not** as
+   `sqrt(max(sum, eps))` — and matched the float64 oracle to 1.8e-07 over the
+   fixed cases. The measured hardware worst case is the same order, 1.835e-07.
+3. **The reduction order is a tree**, matching `__shfl_down_sync`, tested
+   structurally and shown to differ from a sequential sum on a concrete input.
+
+## Declared tolerance
+
+`max_abs: 1e-4` on probabilities. The measured worst case on the RTX 4090 is
+**1.835e-07**, so the declared bound has 545x of headroom and a failure is a real
+failure rather than tolerance noise. Declared in `scoring.backend.json`.
+
+## What the fusion is worth
+
+Measured against an unfused baseline that reads the candidate matrix twice, at
+1024 threads: **1.22x** at CLM's narrowest shape, **1.51x** at 64 candidates and
+**1.79x** to **1.84x** at 255. Kev's plain dot, which has no norm to fuse, lands
+at the same **1.81x** from the second kernel and the global similarity round trip
+it does not need.
+
+The measurement also found a defect the correctness checks could not: at the 256
+threads this kernel was first written with, a one-question request occupies one of
+the 128 SMs, and the norm accumulation sat on that block's critical path — so the
+fusion *lost*, down to 0.57x on Kev's shape. Protocol, both result sets and the
+limits are in [the benchmark](../../../docs/benchmarks/scoring-fusion/README.md).
+
+## Reproducing
+
+On a machine with `nvcc` and an NVIDIA GPU:
+
+```sh
+src/backends/cuda/scoring/build.sh ./out 89   # or CUDA_ARCH_LIST="80 89 90"
+python3 src/backends/cuda/scoring/gpu_parity.py --library ./out/libscoring.so
+python3 src/backends/cuda/scoring/test_scoring.py   # the reference, on CPU
+python3 src/backends/cuda/scoring/bench.py --library ./out/libscoring.so   # compiled kernel vs baseline
+```
+
+The inputs are not committed. They are deterministic from a fixed seed, so the
+harness regenerates them on every run; `reference.py --json vectors.json` dumps
+them if a run needs to be archived. Committing a megabyte of generated floats to
+check a kernel that regenerates them is weight without evidence.
+
+`gpu_parity.py` skips rather than fails when there is no library or no device, so
+a machine without a GPU reports that it could not check, not that it passed.
+
+## Why not a generic kernel
+
+The two models disagree on details that a generic tensor abstraction would have
+to hide and then re-expose: whether there is a projection before the similarity,
+whether the norm floor applies to the query, the row, or both, and what the
+per-question candidate count is. This kernel takes the *shape* they share and
+leaves the projections to the caller, which is where `#6` puts kernel selection
+anyway — with the model engine.