[rafaelha]

compile_sampler and compile_detector_sampler take device="cuda" (with strategy="active_statevector"). With it, the noise sampler and every statevector component run as CUDA kernels. CuPy compiles them with NVRTC, in float32, for components of up to 11 active qubits.

What's in it

• Noise kernel. Each block of shots keeps its shot words in shared memory. Its threads walk each fire-probability group's pooled trials with geometric skips and record the faults. A second pass applies them, one 32-byte record per outcome. The block screens its shots against the postselected direct detectors and queues the survivors.
• Tape kernels. Each component's tape is generated as straight-line code, with the amplitudes in registers. A shot uses up to 32 threads, and warp shuffles fetch partner amplitudes. A tape longer than a span budget is emitted as loops instead, so wide, long tapes compile in seconds. Only the screened survivors run a tape. Outcomes are drawn per (seed, shot, action), as the CPU executor draws them. An outcome whose weight is at most 1e-11 of the total counts as impossible, so float32 rounding cannot create branches; the others are drawn against all 64 bits of a draw.
• CompiledDetectorSampler.sample_postselected(shots, *, postselection_mask, observable_expectations, ...). This is the only new public method. It returns only the shots that pass postselection, bit-packed, so the GPU copies only those rows, once, into pinned memory. shots counts every shot, so the acceptance rate is len(rows) / shots. It also works on the CPU. With observable_expectations=True (CUDA only), it also returns each observable's expectation of (−1)^obs on the shot's state before the measurement that settles it, as Clifft's EXP_VAL does. That expectation ignores postselection that comes later, so when a postselected detector of the same component is tested at or after that measurement, the shot's sampled ±1 is returned instead. Either way, the mean over the returned shots is an unbiased estimate of the observable's expectation.
• sample and sample(..., postselection_mask=...). Both keep their contract on the GPU and return every row.
• Install. A new gpu extra (cupy-cuda13x[ctk]), which needs CUDA 13 (driver 580 or newer). runtime.py preloads the pip NVRTC 13, because CuPy otherwise loads the first NVRTC on the library path.

Not supported on the GPU:

• trajectory sharing;
• components wider than 11 active qubits (they raise ValueError).

The stabilizer-rank strategies and the JAX cuda12/cuda13 extras are unchanged.

Performance (RTX 5090)

Measured with near-clifford-benchmarks: 2^22 shots per call, excluding compilation and warmup, with the device-to-host copy of the returned rows included. Rates are in shots/s.

circuit	k	all shots	postselected
cultivation_d3	4	1.3e9	1.9e9
cultivation_d5	6	6.9e8	3.2e9
cultivation_d5_end2end	6	8.1e7	3.8e8
cultivation_d7	8	1.3e8	1.4e9
distillation_d5	3	1.1e9	7.2e9
coherent_d5_r5	11	8.4e5	9.7e6

• Postselected runs are faster because discarded shots stop at screening and are never copied.
• End-to-end d=5: the postselected rate (3.8e8) is about 450× the 1-thread CPU executor with trajectory sharing (8.4e5) on a Ryzen 9 9900X. That CPU figure predates the numba-cache fix that landed with Speed up the statevector amplitude kernels #213, so the ratio needs re-measuring.
• Compile time: a few seconds per program (NVRTC), cached per postselection mask.

Testing

• test/integration/test_cuda_sampler.py has 45 tests, which skip without CuPy and a GPU. They compare detector, measurement and postselected samples against the stabilizer-rank strategy and the CPU executor. They also cover reference samples, programs without statevector components, expectations (including their fallback under later postselection), widths 8, 9 and 11 (8, 16 and 32 threads per shot), the width limit, and small outcome probabilities drawn at their rate (p = 1e-10 at 2^28 shots).
• The full suite passes locally except test_kernels.py::test_the_calling_thread_keeps_parallel_kernels_after_run_tasks, which fails the same way without this PR under numba's TBB threading layer (fixed separately on fix/thread-spread-test).
• Statistical cross-checks at scale:
  • d=5 end-to-end cultivation reproduces Clifft's published curve at every compared gap threshold (1.1e11 shots).
  • tsim GPU, tsim CPU and stim agree on Hetényi's d=3 gauging circuit at 1e10 shots each.

CI has no GPU, so the CUDA tests only run locally.

🤖 Generated with Claude Code

[github-actions]

PR Preview Action v1.8.2
🚀 View preview at https://queracomputing.github.io/tsim/pr-preview/pr-219/
Built to branch gh-pages at 2026-10-28 04:46 UTC. Preview will be ready when the GitHub Pages deployment is complete.

[github-actions]

☂️ Code Coverage

current status: ✅

Overall Coverage

Statements	Covered	Coverage	Threshold	Status
4554	3915	86%	0%	🟢

New Files

File	Coverage	Status
src/tsim/statevector/cuda/_init_.py	100%	🟢
src/tsim/statevector/cuda/runtime.py	69%	🟢
TOTAL	84%	🟢

Modified Files

File	Coverage	Status
src/tsim/circuit.py	92%	🟢
src/tsim/core/types.py	100%	🟢
src/tsim/sampler.py	84%	🟢
TOTAL	92%	🟢

updated for commit: 840be38 by action🐍

[rafaelha]

Review of the CUDA backend. The first three comments change sampled statistics relative to device="cpu"; the next two are memory-scaling failures; the last is test coverage.

[rafaelha]

Geometric skip is biased for small fire probabilities. unit(r) is a 24-bit uniform and __logf is the fast approximation, so the per-segment fire probability has a floor of 2^-24 (u = 1.0 always gives skip 0). The CPU next_skip is exact down to 2^-64.

Measured on an RTX 5090 over all 2^24 unit() values, 256-trial segment, one channel per block:

p	P(fire) / exact
1e-8	1.02x
1e-10	4.7x
1e-12	233x

A circuit with many negligible-but-nonzero channels (e.g. 1000 idle channels at p=1e-12) gets ~230x the CPU fault rate, with no warning. __logf is never positive for any unit() output, so there is no negative-skip / out-of-bounds issue, only the bias. Suggest a 53-bit uniform with log/log1p in double, or at least a 32-bit uniform on (0,1).

[rafaelha]

Samples depend on batch_size. Component seeds and the noise seed are drawn from self._rng once per device batch, and the noise stream is keyed by blockIdx within the batch. So compile_sampler(seed=1, device="cuda").sample(2**23) and .sample(2**23, batch_size=2**22) return different samples.

This contradicts the active_statevector docstring in sampler.py ("the sampled outcomes of a seed do not depend on the batch size") and ShotSampler.sample's batch_size doc. Either key seeds/noise by the global shot index, or document that device="cuda" is exempt.

[rafaelha]

RNG stream collides with the reference sample. shots.py uses default_rng([seed, 2]) for _reference_seeds. Here the same stream is used, and both the host reference() and the first _run_batch start with integers(0, 2**64, size=num_components), so they get identical component seeds. draw(seed, 0, ...) then matches the CPU _draw for shot 0.

With use_detector_reference_sample=True, shot 0 of the first batch reproduces the noiseless reference whenever noise misses the component, so its detectors are deterministically 0 after the XOR. Use a distinct stream, e.g. [seed, 3].

[rafaelha]

Whole-call host allocation. _sample_words allocates (shots, num_words) uint64 up front, including key words, and copies every word back per batch. ShotSampler.sample on CPU splits into _part_shots() parts bounded by memory_budget() and keeps only output columns.

E.g. compile_detector_sampler(device="cuda").sample(10**9) on a 12-word layout tries to allocate 96 GB and raises MemoryError before any kernel runs; the same call on CPU succeeds. The gather kernel (compiled but never launched) looks like it was meant for this: gather only output_words per batch and size the host array accordingly.

[rafaelha]

Batch size ignores device memory and word count. A layout with 150 shot words (the noise kernel admits up to ~176) needs 2^22 × 150 × 8 B ≈ 5 GB for words alone, plus queues and the float32 values buffer. On an 8–12 GB GPU this hits OutOfMemoryError unless the user guesses a smaller batch_size. Suggest deriving the default from cp.cuda.Device().mem_info and layout.num_words, analogous to memory_budget() on CPU.

[rafaelha]

Module-wide skip hides device-independent tests in CI. test_sample_postselected_returns_only_passing_shots[cpu], test_expectations_need_the_cuda_device, and the device="tpu" / strategy="cat5" ValueError checks in test_cuda_needs_the_statevector_strategy don't need CuPy, but CI (no GPU) skips them. A regression in the CPU branch of sample_postselected would go unnoticed. Suggest moving these to test/unit/test_sampler.py or applying the skip only to the CUDA cases.