feat(sampler): choose the sampling device with a device argument

.claude/skills/benchmarking-tsim/SKILL.md

@@ -100,10 +100,25 @@ uv run tsim-dev bench --side .:device=gpu --side .:device=cpu,threads=10   # GPU
 ```
 
 - The `device` column must say `metal` on the GPU side. It says `cpu` where Metal falls
-  back: a component wider than 11, noise beyond the GPU tables, no GPU memory.
+  back: a component wider than 11, a rotation branch below float32's resolution (an angle
+  below about 6e-5 π), noise beyond the GPU tables, no GPU memory. Under the default
+  `device="auto"` these fallbacks are silent; in a script of your own, compile with
+  `device="metal"` to make them raise instead.
+- `TSIM_DEVICE=cpu` (the bench's default) keeps `device="auto"` on the CPU; an explicit
+  `device=` argument overrides it.
+- The GPU's rate depends on call size much more than the CPU's: each call pays a fixed cost
+  (setting up the launches, waiting for the GPU to finish), and each batch of up to 2^20 shots is one
+  launch. The bench calibrates calls to about 0.3 s, which is enough. A GPU rate measured
+  with calls of 10^5 shots or fewer, or with a small `batch_size`, measures overhead, not
+  sampling; never report one. Leave `batch_size` unset on the GPU; it can only lower the
+  batch, and it changes a seed's GPU samples.
 - A reference without the Metal backend samples on the CPU even with `--device gpu`; its
   rows end in `DEVICE: cpu -> metal`, and the ratio compares devices, not code.
 - Compare the GPU with all CPU cores (`--threads N`, from the `NOTE` line), not one.
+- Say whether rates are bit-packed (the default) or `--bools`, above all on a GPU: booleans
+  cost it 2.3-3.2x on circuits with hundreds of detectors (surface_d7_r7: 757M packed,
+  237M booleans on an M4 Pro), against 1.0-1.8x on the CPU, so a boolean comparison
+  understates the GPU's lead. Compare Tsim with other simulators on the same output format.
 - The GPU's own time per kernel: set `tsim.statevector.metal.runtime.PROFILE = {}` before
   a `sample` call, which then runs each dispatch alone and adds its GPU seconds under the
   kernel's name (`noise`, `tape`, `gather_rows`, ...).

.claude/skills/profiling-tsim/SKILL.md

@@ -31,19 +31,24 @@ Tsim is slow:
   sampler returns its one observable at 60.9M shots/s per core, the measurement
   sampler its 10 measurements at 11.6M (Apple M4 Pro).
 - **Output.** The harness times bit-packed rows, Tsim's own output. `--bools`
-  times `sample()`'s default boolean output, which adds a conversion: 1.6x on
-  surface_d7_r7 on the CPU. A slow boolean rate next to a fast bit-packed one
-  points at the conversion, not the sampler.
+  times `sample()`'s default boolean output, which adds a conversion: 1.6-1.8x on
+  surface_d7_r7 on the CPU, and 2.3-3.2x on a GPU for circuits with hundreds of
+  detectors (Metal on an M4 Pro: surface_d7_r7 3.2x, end2end_d3_dsurface15 2.8x,
+  hook_inject_T_d15_pd5 2.3x; under 1.3x for circuits with few detectors). A slow
+  boolean rate next to a fast bit-packed one points at the conversion, not the
+  sampler.
 - **Threads.** One pinned thread gives the rate per core. Tsim samples on all
   performance cores by default (`numba.get_num_threads()`); the `NOTE` line says
   how many. Quote all-core rates with `--threads N`.
 - **Device.** The harness sets `TSIM_DEVICE=cpu`, so every side samples on the
   CPU unless you pass `--device gpu` (or `--side REF:device=gpu`). On an
   Apple-silicon Mac, Tsim samples on the GPU by default for k_max <= 11, 4-19x
   faster than all CPU cores; there `--device gpu` measures what users get. The
-  `device` column says where each row sampled: a GPU falls back to the CPU for
-  components wider than it runs, and a row that ends in `DEVICE: cpu -> metal`
-  compares a GPU with the CPU.
+  `device` column says where each row sampled: under `device="auto"` a GPU falls
+  back to the CPU silently (components wider than 11, tiny rotations, noise beyond
+  its tables), and a row that ends in `DEVICE: cpu -> metal` compares a GPU with the
+  CPU. A slow GPU row with small calls or a small `batch_size` is launch overhead;
+  see "Apple-silicon Macs" in the `benchmarking-tsim` skill.
 - **Compile.** `compile` includes the first use of a backend in the process (about
   0.1 s for Metal's PyObjC) but not the first call, whose kernel compilation goes
   into `warmup_s` in the JSON.

.claude/skills/testing-tsim-correctness/SKILL.md

@@ -27,7 +27,9 @@ handful of T gates; its memory grows exponentially.
 
 ## Invariants to test when you touch them
 
-- Same seed gives the same samples for any thread count and any `batch_size`.
+- Same seed gives the same samples for any thread count and any `batch_size` on the
+  CPU. On a GPU the samples depend on `batch_size`; same seed and same `batch_size`
+  give the same samples.
 - Bit-packed and unpacked outputs agree. Every `sample()` call draws new shots, so
   compare two fresh samplers with the same seed, not two calls on one sampler.
 - Postselection returns only the accepted rows; compare accepted-row distributions.
@@ -44,6 +46,11 @@ handful of T gates; its memory grows exponentially.
   test is a broken test.
 - Stim merges consecutive identical instructions into one with several targets, so
   write circuits that exercise both single- and multi-target paths.
+- GPU tests: `test/conftest.py` keeps every test outside
+  `test/integration/test_metal_sampler.py` on the CPU. In the Metal tests compile the
+  GPU side with `device="metal"`, so a silent fallback cannot make the test pass on the
+  CPU, and the reference with `device="cpu"`; compare them with `assert_counts_match`
+  (GPU samples never equal CPU samples bit for bit).
 - Keep circuits small (a few qubits, tens of thousands of shots) so the suite stays
   near 90 s. Long randomized campaigns belong in a job, not in the suite:
   `uv run tsim-dev job start fuzz --mem 4G -- <command>`.

CHANGELOG.md

@@ -11,7 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `strategy="active_statevector"` samples on the GPU of an Apple-silicon Mac, in
   float32, whenever every component has at most 11 active qubits (Tsim now depends
   on `pyobjc-framework-Metal` there). A seed's GPU samples differ from the CPU's and
-  depend on `batch_size`; `TSIM_DEVICE=cpu` keeps the samplers on the CPU, and
+  depend on `batch_size`. The `device` argument of `compile_sampler`,
+  `compile_detector_sampler` and the sampler classes chooses where they sample:
+  `"auto"` (the default) uses the GPU when it can run the circuit and otherwise
+  the CPU, without a warning; `"cpu"` keeps them on the CPU; `"metal"` raises when
+  the GPU cannot run the circuit. `TSIM_DEVICE=cpu` keeps `"auto"` on the CPU, and
   `SamplerStats.device` and the repr say where they sample.
 - The docs pages have a menu to copy a page or open it in ChatGPT or Claude, and
   the README, the install page and the agent skill point to the docs MCP server

docs/install.mdx

@@ -38,30 +38,47 @@ pip install "bloqade-tsim[cuda13]"
 ## Apple GPUs
 
 On an Apple-silicon Mac, the default strategy (`"active_statevector"`) samples on the GPU
-through Metal; the plain install includes PyObjC's Metal bindings for it. The samplers use the
-GPU whenever every component of the circuit has at most 11 active qubits, and otherwise sample
-on the CPU, with a warning. Postselection runs on the GPU, so heavily postselected circuits
-gain the most. A sampler's repr says where it samples (`samples on metal` or `samples on
-cpu`), and so does `sampler.stats.device`.
+through Metal; the plain install includes PyObjC's Metal bindings for it. Where it samples is
+the `device` argument of `compile_sampler` and `compile_detector_sampler`:
+
+| `device` | Samples |
+| --- | --- |
+| `"auto"` (default) | On the GPU when it can run the circuit, else on the CPU, without a warning |
+| `"cpu"` | On the CPU, in double precision |
+| `"metal"` | On the GPU, or raises an error that says why it cannot |
+
+```python
+sampler = circuit.compile_detector_sampler(device="metal")
+```
+
+The GPU runs circuits whose components have at most 11 active qubits. That includes Clifford
+circuits, whose detectors depend on the noise alone: the GPU samples the noise as well as the
+state vectors. Postselection runs on
+the GPU, so heavily postselected circuits gain the most. A sampler's repr says where it
+samples (`samples on metal` or `samples on cpu`), and so does `sampler.stats.device`.
 
 The GPU computes outcome probabilities in float32, so a sample follows the CPU's distribution
 to within float32 rounding. Float32 cannot tell a measurement outcome of probability below
 about 1e-11 from rounding error, so the GPU never draws one; a circuit with a rotation small
-enough to make such an outcome (an `R_Z` or `R_X` by less than about 1e-4 π) samples on the
-CPU, with a warning. A seed's samples on the GPU differ from the CPU's and depend on
-`batch_size`.
+enough to make such an outcome (an `R_Z` or `R_X` by less than about 6e-5 π) samples on the
+CPU (or raises, with `device="metal"`).
+
+A seed's samples on the GPU differ from the CPU's and depend on `batch_size`, so with
+`device="auto"` a seed reproduces its samples only on the same kind of machine. Use
+`device="cpu"` to reproduce samples across machines.
 
 ### Sampling on the CPU instead
 
-Set `TSIM_DEVICE=cpu` to keep the samplers on the CPU, in double precision, for example to
-reproduce samples from another machine:
+Set `TSIM_DEVICE=cpu` to keep the samplers of `device="auto"` on the CPU without changing
+code, for example to reproduce samples from another machine:
 
 ```bash
 TSIM_DEVICE=cpu python my_script.py
 ```
 
 Tsim reads it whenever it compiles a sampler, so `os.environ["TSIM_DEVICE"] = "cpu"` before
-`compile_sampler` or `compile_detector_sampler` works too.
+`compile_sampler` or `compile_detector_sampler` works too. An explicit `device` argument
+overrides it.
 
 ## Agent skill
 

skills/simulating-circuits-with-tsim/SKILL.md

@@ -104,8 +104,9 @@ If `sampler.active_width` raises `AttributeError`, the installed Tsim is older t
 upgrade `bloqade-tsim` before going further.
 
 Before a long run, time it: one small warm-up call (the first call compiles kernels), then
-one call of about 10^6 shots, and extrapolate. Small calls understate the rate. Report the
-width and the measured shots per second, with the thread count, along with any plan.
+one call of about 10^6 shots, and extrapolate. Small calls understate the rate, on a GPU even
+more. Report the width, the measured shots per second and where it ran (`sampler.stats.device`,
+and the thread count on the CPU), along with any plan.
 
 If the width is too large: postselect inside Tsim (discarded shots stop early, see below),
 simulate a smaller instance (lower distance or fewer rounds) and report how the width grows,
@@ -176,27 +177,49 @@ yourself.
 - **Shots per call:** each `sample()` call has a fixed overhead. For fast circuits, ask for
   about 10^6 or more shots per call, limited by memory, instead of looping over small calls.
   Calls of 10^5 shots have been measured 2 to 4 times slower than calls of 10^7, and 100-shot
-  calls about 80 times slower.
+  calls about 80 times slower. On a GPU the fixed cost per call is larger (setting up the
+  launches and waiting for the GPU to finish), so ask for 10^6 to 10^7 shots per call there.
 - **`batch_size`:** leave it at the default, which is sized from free memory and is close to
   optimal. Small values (a few thousand or less) are several times slower. Lower it only to
-  cap memory; it does not change the samples.
+  cap memory; on the CPU it does not change the samples. On a GPU the default batch is up to 2^20
+  shots (fewer for very wide rows), one GPU launch, all launched back to back; `batch_size` can only lower it, a small
+  one makes launch overhead dominate, and it changes which samples a seed gives.
 - **`bit_packed=True`** returns little-endian packed `uint8` rows, as in Stim, using 8 times
   less memory. On a distance-11 surface code memory (1,320 detectors) it was twice as fast as
-  booleans, at distance 5 about 1.3 times. For small circuits it gains little and is slower
+  booleans, at distance 5 about 1.3 times. On an Apple GPU the gap is larger, since booleans
+  move 8 times more bytes off the GPU: 2.3 to 3.2 times on circuits with hundreds of
+  detectors (surface_d7_r7: 757M against 237M shots/s on an M4 Pro). Use packed rows for
+  large studies and decode them packed. For small circuits it gains little and is slower
   with `separate_observables` (the observables are split off and re-packed); use
   `append_observables=True` with packed output.
 - **Threads:** the default strategy runs on `numba.get_num_threads()` threads, all cores by
   default. Use `numba.set_num_threads(k)` before `sample` (or `NUMBA_NUM_THREADS=k` before
   importing Tsim) and leave cores free on a shared machine. Do not run several multi-threaded
   Tsim processes side by side; they oversubscribe the CPU.
 - **Apple GPUs:** on an Apple-silicon Mac the default strategy samples on the GPU when every
-  component has at most 11 active qubits, often 5 to 10 times faster than all CPU cores. It
+  component has at most 11 active qubits, often 5 to 10 times faster than all CPU cores. That
+  includes Clifford circuits, whose detectors come from the noise alone: the GPU samples the
+  noise too. It
   computes outcome probabilities in float32 and never draws an outcome below about 1e-11
-  (circuits with rotations by less than about 1e-4 π stay on the CPU, with a warning); set
-  `TSIM_DEVICE=cpu` before compiling to stay on the CPU, in double precision. The sampler's repr says where it samples (`samples on metal`).
+  (circuits with rotations by less than about 6e-5 π stay on the CPU). Thread settings do not
+  affect GPU sampling.
+- **`device`** (in `compile_sampler` and `compile_detector_sampler`) chooses where it samples:
+  `"auto"` (default) uses the GPU when it can and otherwise the CPU, silently; `"cpu"` stays
+  on the CPU in double precision; `"metal"` raises with the reason if the GPU cannot run the
+  circuit (no Apple GPU, width above 11, a tiny rotation, out of GPU memory). Check
+  `sampler.stats.device` (or the repr, `samples on metal`) after compiling, since `"auto"`
+  does not warn. Pass `device="metal"` for long runs that must not drop to the CPU unnoticed,
+  and `device="cpu"` to reproduce samples across machines or to compare against the GPU.
+  `TSIM_DEVICE=cpu` in the environment keeps `"auto"` on the CPU; it is read when a sampler
+  is compiled, so an unpickled sampler keeps the device choice it was compiled with.
+- **Processes:** sample from one process per GPU. Threads sharing a GPU sampler run its calls
+  one at a time, and several processes compete for the same GPU. A process forked after its parent used the GPU cannot use
+  it (it samples on the CPU, or raises with `device="metal"`); start workers with the
+  `"spawn"` or `"forkserver"` method (macOS's default is `"spawn"`).
 - **Seeds:** a `seed` reproduces samples for the same Tsim and dependency versions and
-  backend. On the CPU, thread count and `batch_size` do not affect the samples; on a GPU,
-  `batch_size` does.
+  device. On the CPU, thread count and `batch_size` do not affect the samples; on a GPU,
+  `batch_size` does, and GPU samples differ from CPU samples, so with `device="auto"` a seed
+  reproduces only on the same kind of machine.
 
 ## Postselection
 
@@ -286,7 +309,7 @@ tick it off:
 
 ```
 Tsim results:
-- [ ] Active width read; rate timed with a warm-up and a large call; thread count noted
+- [ ] Active width read; rate timed with a warm-up and a large call; device and threads noted
 - [ ] Noiseless run: detectors used for postselection or decoding never fire
 - [ ] Reference flags used only if detectors and observables are deterministic
 - [ ] At least one check below passed (stand-in vs Stim, to_tensor, independent simulator)

skills/simulating-circuits-with-tsim/references/qec-experiments.md

@@ -56,7 +56,7 @@ sinter can call Tsim, but it costs speed and can cost correctness:
 - sinter's decoders expect Stim's flips relative to a noiseless reference; Tsim returns raw
   values, so a detector that is 1 without noise looks like an error.
 - Each sinter worker process compiles its own sampler and starts its own numba thread pool, so
-  several workers oversubscribe the CPU.
+  several workers oversubscribe the CPU, or queue on the one GPU of an Apple-silicon Mac.
 
 Write the short loop below instead. sinter's helpers are still useful for the statistics
 (`sinter.fit_binomial`) and for plotting.
@@ -124,7 +124,8 @@ Notes:
   `max_seconds` as caps, rather than fixing the shot count in advance: logical error rates are
   often orders of magnitude away from the first guess.
 - Calls are sized from a memory budget (`max_bytes_per_call`, 1 GB by default), not a fixed
-  shot count: 10^6 shots of a circuit with 50,000 detectors is a 50 GB array.
+  shot count: 10^6 shots of a circuit with 50,000 detectors is a 50 GB array. Keep calls
+  large on a GPU as well: below about 10^6 shots per call its fixed cost per call shows.
 - `run_point` works on boolean rows, which the decoder recipes expect. For circuits with
   hundreds of detectors or more, packed rows use 8 times less memory and can be twice as fast;
   to use them, change the loop to sample with `bit_packed=True` and `append_observables=True`,