feat(statevector): sample on Apple GPUs through Metal

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

@@ -82,6 +82,32 @@ uv run tsim-dev bench --side .:threads=4 --side origin/main:threads=4
 - Compile time is compared too; when the sides differ only in threads, read the rate
   column.
 
+## Apple-silicon Macs (Metal)
+
+On a Mac with an M-series chip, the default strategy samples on the GPU through Metal
+whenever every component has at most 11 active qubits (k_max <= 11). A plain `uv sync`
+installs it. That is what users of a MacBook get, and it is much faster than the CPU: on an
+M4 Pro, 26-96x one CPU thread and 4-19x all ten performance cores, bit-packed. End-to-end
+cultivation d=5 ran at 125M shots/s, 19x the ten cores. Circuits without noise whose CPU
+trajectories share well, such as qv_q10, run at the CPU's rate.
+
+The bench sets `TSIM_DEVICE=cpu` by default, so on a Mac a plain run measures the CPU only.
+Say which you measured, and on a Mac measure both:
+
+```bash
+uv run tsim-dev bench --device gpu --ref REF   # a GPU change; REF must have the Metal backend
+uv run tsim-dev bench --side .:device=gpu --side .:device=cpu,threads=10   # GPU vs all cores
+```
+
+- 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.
+- 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.
+- 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`, ...).
+
 ## What the command controls for you
 
 - One thread, pinned to a core (`--threads N` to change; say so when you report).

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

@@ -38,8 +38,9 @@ Tsim is slow:
   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`). With a GPU
-  backend installed, `--device gpu` measures what users get by default. 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.
@@ -61,7 +62,8 @@ uv run tsim-dev profile cultivation_d5 --seconds 5
 The run is a fresh process, one pinned thread, warm numba cache, under a memory cap,
 with bit-packed output on the CPU (`--threads`, `--device` and `--sampler` as in
 step 2). py-spy sees the host only: on a GPU the sampling phase shows the wait for
-the device, and the time inside its kernels needs the backend's own timers.
+the device. For Metal's time per kernel, see "Apple-silicon Macs" in the
+`benchmarking-tsim` skill (`runtime.PROFILE`).
 The output splits samples by phase (`compile`, `sample`) and prints three tables:
 
 - **innermost Tsim frame (kernel callers):** read this first. Jitted numba code has

AGENTS.md

@@ -66,9 +66,10 @@ more than a handful of T gates: their compile memory grows exponentially.
   may change in any release; every other module is internal.
 - **Raw outputs:** detectors and observables are raw parities, unlike Stim's flips relative to
   a reference. The reference-sample flags change only that.
-- **Seeds:** with `active_statevector`, a seed's samples do not depend on the thread count or
-  `batch_size`. Do not break this; noise blocks and per-shot draws are seeded by call and shot
-  index (`shots.py`).
+- **Seeds:** with `active_statevector` on the CPU, a seed's samples do not depend on the
+  thread count or `batch_size`. Do not break this; noise blocks and per-shot draws are seeded
+  by call and shot index (`shots.py`). On a GPU (`statevector/metal/`) they depend on
+  `batch_size`, as documented, and differ from the CPU's.
 - **Direct outputs stay direct:** outputs that depend on the noise alone never go through a
   state vector. Most detectors of QEC circuits are direct, which is why the detector sampler is
   fast.

CHANGELOG.md

@@ -8,6 +8,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- `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
+  `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
   at `https://bloqade.quera.com/tsim/mcp`. The README and the install page also

README.md

@@ -23,8 +23,8 @@ Or with [uv](https://docs.astral.sh/uv/):
 uv add bloqade-tsim
 ```
 
-The default sampling strategy runs on the CPU. To run the stabilizer-rank strategies on an NVIDIA GPU,
-install the JAX CUDA extra:
+The default sampling strategy runs on the CPU, and on the GPU of an Apple-silicon Mac (no extra
+needed). To run the stabilizer-rank strategies on an NVIDIA GPU, install the JAX CUDA extra:
 ```bash
 pip install "bloqade-tsim[cuda13]"   # or bloqade-tsim[cuda12]
 ```
@@ -82,7 +82,7 @@ Choose a strategy with the `strategy` argument of `compile_sampler` and `compile
 
 | Strategy | Method | Runs on | Cost grows with | Use it for |
 |---|---|---|---|---|
-| `"active_statevector"` (default) | Dense state vector over the active qubits of each component, as in [Clifft](https://arxiv.org/abs/2604.27058) and [SymFT](https://arxiv.org/abs/2607.28600) | CPU (numba, multi-threaded) | Number of active qubits (at most 30) | Most circuits, especially with postselection |
+| `"active_statevector"` (default) | Dense state vector over the active qubits of each component, as in [Clifft](https://arxiv.org/abs/2604.27058) and [SymFT](https://arxiv.org/abs/2607.28600) | CPU (numba, multi-threaded), or an Apple GPU (Metal) | Number of active qubits (at most 30) | Most circuits, especially with postselection |
 | `"cat5"`, `"bss"`, `"cutting"` | Stabilizer rank decomposition | CPU or GPU (JAX) | Number of non-Clifford gates | Few non-Clifford gates acting on many qubits; GPU throughput |
 
 ```python

codecov.yml

@@ -0,0 +1,4 @@
+# NOTE: the Metal executor needs an Apple GPU, which CI does not have; its tests
+# run on Apple silicon only.
+ignore:
+  - "src/tsim/statevector/metal"

docs/index.mdx

@@ -39,8 +39,8 @@ samples = detector_sampler.sample(shots=100)
 uv add bloqade-tsim
 ```
 
-The default sampling strategy runs on the CPU. To run the stabilizer-rank strategies on an NVIDIA GPU,
-install the JAX CUDA extra:
+The default sampling strategy runs on the CPU, and on the GPU of an Apple-silicon Mac (no extra
+needed). To run the stabilizer-rank strategies on an NVIDIA GPU, install the JAX CUDA extra:
 
 ```bash
 uv add "bloqade-tsim[cuda13]"
@@ -55,7 +55,7 @@ Choose a strategy with the `strategy` argument of `compile_sampler` and `compile
 
 | Strategy | Method | Runs on | Cost grows with | Use it for |
 |---|---|---|---|---|
-| `"active_statevector"` (default) | Dense state vector over the active qubits of each component, as in [Clifft](https://arxiv.org/abs/2604.27058) and [SymFT](https://arxiv.org/abs/2607.28600) | CPU (numba, multi-threaded) | Number of active qubits (at most 30) | Most circuits, especially with postselection |
+| `"active_statevector"` (default) | Dense state vector over the active qubits of each component, as in [Clifft](https://arxiv.org/abs/2604.27058) and [SymFT](https://arxiv.org/abs/2607.28600) | CPU (numba, multi-threaded), or an Apple GPU (Metal) | Number of active qubits (at most 30) | Most circuits, especially with postselection |
 | `"cat5"`, `"bss"`, `"cutting"` | Stabilizer rank decomposition | CPU or GPU (JAX) | Number of non-Clifford gates | Few non-Clifford gates acting on many qubits; GPU throughput |
 
 ```python

docs/install.mdx

@@ -1,6 +1,6 @@
 ---
 title: "Installation"
-description: "Install Tsim with uv or pip, with optional CUDA support."
+description: "Install Tsim with uv or pip, with optional GPU support."
 ---
 
 Tsim requires Python 3.11 or later.
@@ -35,6 +35,34 @@ For GPU acceleration with CUDA:
 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`.
+
+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`.
+
+### 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:
+
+```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.
+
 ## Agent skill
 
 Tsim ships an agent skill that teaches coding agents (Claude Code, Codex, Cursor and others)

pyproject.toml

@@ -43,6 +43,8 @@ dependencies = [
     "stim>=1.15.0",
     "stim-timeline-text>=0.1.0",
     "zxparam>=0.1.0",
+    # NOTE: the statevector executor on Apple GPUs (about 7 MB).
+    "pyobjc-framework-Metal>=10.0; sys_platform == 'darwin' and platform_machine == 'arm64'",
 ]
 
 [project.urls]
@@ -167,6 +169,9 @@ extend_exclude = [
     "^test/",
 ]
 
+[tool.deptry.package_module_name_map]
+pyobjc-framework-Metal = "Metal"
+
 [tool.deptry.per_rule_ignores]
 DEP002 = ["ipywidgets"]
 

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

@@ -189,8 +189,14 @@ yourself.
   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
+  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`).
 - **Seeds:** a `seed` reproduces samples for the same Tsim and dependency versions and
-  backend. Thread count and `batch_size` do not affect the samples.
+  backend. On the CPU, thread count and `batch_size` do not affect the samples; on a GPU,
+  `batch_size` does.
 
 ## Postselection
 

src/tsim/circuit.py

@@ -954,7 +954,9 @@ def compile_sampler(
             strategy: "active_statevector" (the default) to simulate each
                 non-direct component on a dense state vector over its active
                 qubits, or a stabilizer rank decomposition strategy ("cat5",
-                "bss", "cutting"). "active_statevector" samples on the CPU, on
+                "bss", "cutting"). "active_statevector" samples in float32 on an
+                Apple GPU when every component has at most 11 active qubits
+                (``TSIM_DEVICE=cpu`` turns this off), else on the CPU, on
                 ``numba.get_num_threads()`` threads. The stabilizer rank
                 strategies are provisional.
             seed: Random seed for the sampler. If None, a random seed is drawn. A seed
@@ -964,7 +966,8 @@ def compile_sampler(
                 non-direct components, the samples of a seed depend on neither the batch
                 size nor the thread count. With a stabilizer-rank strategy they also
                 need the same explicit ``batch_size``, since the default batch size
-                depends on the free memory.
+                depends on the free memory. On a GPU, they do depend on the batch
+                size.
 
         Returns:
             A CompiledMeasurementSampler that can be used to sample measurements.
@@ -991,7 +994,9 @@ def compile_detector_sampler(
             strategy: "active_statevector" (the default) to simulate each
                 non-direct component on a dense state vector over its active
                 qubits, or a stabilizer rank decomposition strategy ("cat5",
-                "bss", "cutting"). "active_statevector" samples on the CPU, on
+                "bss", "cutting"). "active_statevector" samples in float32 on an
+                Apple GPU when every component has at most 11 active qubits
+                (``TSIM_DEVICE=cpu`` turns this off), else on the CPU, on
                 ``numba.get_num_threads()`` threads. The stabilizer rank
                 strategies are provisional.
             seed: Random seed for the sampler. If None, a random seed is drawn. A seed
@@ -1002,6 +1007,7 @@ def compile_detector_sampler(
                 size nor the thread count. With a stabilizer-rank strategy they also
                 need the same explicit ``batch_size`` and the same reference sample
                 settings, since the default batch size depends on the free memory.
+                On a GPU, they do depend on the batch size.
 
         Returns:
             A CompiledDetectorSampler that can be used to sample detectors and observables.