diff --git a/.github/dependabot.yml b/.github/dependabot.yml
index 89a7a70b8..d795bc1a0 100644
--- a/.github/dependabot.yml
+++ b/.github/dependabot.yml
@@ -61,3 +61,12 @@ updates:
       day: "monday"
     labels:
       - "dependencies"
+    ignore:
+      # The image's Python must be one pyproject.toml supports (requires-python);
+      # a minor bump such as 3.13 -> 3.14 broke that once (#810). Move it by hand
+      # together with requires-python and the CI matrix.
+      # test/repro/test_packaging_consistency.py checks the three agree.
+      - dependency-name: "python"
+        update-types:
+          - "version-update:semver-major"
+          - "version-update:semver-minor"
diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml
new file mode 100644
index 000000000..b4657477a
--- /dev/null
+++ b/.github/workflows/docker.yml
@@ -0,0 +1,37 @@
+# Builds the Dockerfile and runs the reproducibility tests inside the image, so
+# the container is checked, not only described. Runs when the image's inputs change.
+name: Docker image
+
+on:
+  push:
+    branches: [ master ]
+    paths: [ "Dockerfile", "pyproject.toml", "poetry.lock", ".github/workflows/docker.yml" ]
+  pull_request:
+    branches: [ master ]
+    paths: [ "Dockerfile", "pyproject.toml", "poetry.lock", ".github/workflows/docker.yml" ]
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  docker:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v7
+      # The image is about 7 GB (torch bundles its CUDA libraries); make room.
+      - name: Free runner disk space
+        run: sudo rm -rf /usr/share/dotnet /usr/local/lib/android /opt/ghc /opt/hostedtoolcache/CodeQL
+      - name: Build image
+        run: docker build -t bindsnet .
+      # Runners have no GPU, and the image's CUDA torch build can segfault there while
+      # preloading CUDA (see .github/scripts/use_cpu_torch.sh), so the test container
+      # swaps in the CPU build of the same version; everything else is the image as built.
+      - name: Run reproducibility tests in the image
+        run: >-
+          docker run --rm bindsnet bash -c
+          ".github/scripts/use_cpu_torch.sh python && python -m pytest -q test/repro"
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 6b248e8d6..2f1c58446 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -24,7 +24,7 @@ jobs:
     - name: Check release tag matches pyproject version
       if: github.event_name == 'release'
       run: |
-        version=$(python -c "import tomllib; print(tomllib.load(open('pyproject.toml','rb'))['tool']['poetry']['version'])")
+        version=$(python -c "import tomllib; print(tomllib.load(open('pyproject.toml','rb'))['project']['version'])")
         tag="${GITHUB_REF_NAME#v}"
         if [ "$version" != "$tag" ]; then
           echo "Release tag $GITHUB_REF_NAME does not match pyproject version $version"
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8a6708bf5..100b1a678 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,41 @@ see the [GitHub releases / tags](https://github.com/BindsNET/bindsnet/releases).
 ## [Unreleased]
 
 ### Fixed
+- `--seed` now makes GPU runs repeatable in `eth_mnist.py`, `batch_eth_mnist.py`,
+  `supervised_mnist.py`, `conv_mnist.py`, `conv1d_MNIST.py`, `conv3d_MNIST.py`,
+  `SOM_LM-SNNs.py` and `breakout/random_network_baseline.py`: they call
+  `torch.manual_seed(seed)` unconditionally, which seeds the CPU generator (weight
+  initialisation, Poisson encoding, data shuffling) and every GPU. Before, the GPU path
+  called only `torch.cuda.manual_seed_all`; two runs of
+  `batch_eth_mnist.py --seed 0 --n_train 640 --n_test 320` gave 0.37 and 0.47.
+  `dot_tracing.py` now uses its `--seed` for torch and the dot simulator (default
+  seed 1, since the simulator reads 0 as "seed from the clock").
+- `supervised_mnist.py` plots input-to-excitatory weights through
+  `feature_index["weight"]` (`MulticompartmentConnection` has no `.w`), so it runs to
+  the end.
+- `reservoir.py` scores its readout on the MNIST test set (it reused the training
+  loader before) and reports the number of images it actually scored.
+- `conv_mnist.py` now tests the network: it labels output neurons from training spikes
+  (`assign_labels`) and reports all-activity and proportion-weighting test accuracy,
+  like `eth_mnist.py`.
+- `--plot` works as a switch in `eth_mnist.py`, `supervised_mnist.py`, `conv_mnist.py`,
+  `reservoir.py`, `conv1d_MNIST.py`, `conv3d_MNIST.py` and `SOM_LM-SNNs.py`: plotting is
+  off unless `--plot` is given (it was `store_true` with a default of `True`).
+  `conv_mnist.py` no longer pauses 1 s per sample.
+- Training and test loops in the MNIST examples run exactly `--n_train` / `--n_test`
+  samples (`>=` instead of `>`, which ran one more).
+- The Breakout examples run with the locked `ale-py` 0.12.1: `GymEnvironment` maps
+  `BreakoutDeterministic-v4` to its ale-py 0.8.1 definition, `Breakout-v4` with
+  `frameskip=4` (newer ale-py no longer register that id).
+- `play_breakout_from_ANN.py` loads `trained_shallow_ANN.pt` under PyTorch 2.6+: the
+  file is now a `state_dict` with identical tensors, loaded with `weights_only=True`
+  (it was a pickled `Net` object).
+- `publish.yml` reads the version from `[project]`, where it has lived since 0.3.4, so
+  a published GitHub release passes its tag check (it read `[tool.poetry]`; the last
+  PyPI upload was started by hand, which skips the check).
+- Dockerfile on `python:3.13-slim`, inside `requires-python = ">=3.11,<3.14"` (#810 had
+  moved it to 3.14); Dependabot no longer proposes minor or major Python bumps for the
+  image.
 - API reference on Read the Docs was empty: the build never installed `bindsnet`, so
   every `automodule` failed to import (`No module named 'matplotlib'`), and
   `docs/pyproject.toml` downgraded Sphinx to 7.2.6. `.readthedocs.yaml` now installs a
@@ -22,6 +57,26 @@ see the [GitHub releases / tags](https://github.com/BindsNET/bindsnet/releases).
   `environment.cue_reward`, `environment.dot_simulator` and
   `analysis.dotTrace_plotter`, which were missing.
 
+### Added
+- `test/repro/test_packaging_consistency.py`: checks that the Dockerfile's Python is in
+  `requires-python`, that CI tests every supported Python, that `CITATION.cff` and
+  `pyproject.toml` carry the same version, that `publish.yml` reads that version, and
+  that every example with `--seed` seeds torch unconditionally.
+- `.github/workflows/docker.yml`: builds the image and runs `test/repro` inside it when
+  the Dockerfile, `pyproject.toml` or `poetry.lock` change.
+- `examples/benchmark/bindsnet_scaling.py`: the README scaling benchmark for BindsNET
+  alone (no other simulator needed).
+- `examples/benchmark/benchmark.py` rewritten: every simulator (BindsNET CPU/GPU,
+  BRIAN2, brian2genn, PyNEST, ANNarchy CPU/GPU) gets the network the README describes,
+  each in its own process, and a simulator that is not installed is skipped. It
+  replaces the 2018 script, which imported an `experiments` helper package that was not
+  released and used its own input settings per simulator: BindsNET 0-1 Hz with U(0, 1)
+  weights, BRIAN2 15 Hz with U(0, 0.01) weights and no `on_pre` (so its synapses did not
+  drive the membrane), PyNEST 60 Hz connected one-to-one, ANNarchy 50 Hz with U(0, 1)
+  weights. `annarchy.py` and `gpu_annarchy.py` merged into it; `plot_benchmark.py` reads
+  its CSV.
+- `test/environment/test_gym_environment.py`.
+
 ### Changed
 - CI: one test workflow (`python-app.yml`: job `build` on Python 3.13 plus a 3.11/3.12
   matrix, Poetry 2.4.3 with dependency caching, superseded runs cancelled);
diff --git a/DATA.md b/DATA.md
index 2a0e5b307..47e33b3f1 100644
--- a/DATA.md
+++ b/DATA.md
@@ -31,7 +31,9 @@ the library provides.
 
 ### Atari — Breakout (and Space Invaders)
 - **Loader:** `bindsnet.environment.GymEnvironment("BreakoutDeterministic-v4")`
-  (see `examples/breakout/*.py`).
+  (see `examples/breakout/*.py`). `GymEnvironment` maps this id to its definition,
+  `Breakout-v4` with `frameskip=4`, so it works with ale-py 0.9 and later, which
+  register Breakout under `Breakout-v4`.
 - **Upstream source:** Arcade Learning Environment via `gymnasium[atari]` + `ale-py`
   (declared in `pyproject.toml`). ROMs are provided through the ALE/AutoROM tooling.
 - **Obtained by:** the Gymnasium/ALE runtime; not stored in this repo.
@@ -49,14 +51,18 @@ the library provides.
 ## 2. Synthetic stimuli (no external dataset)
 
 ### Scaling-benchmark Poisson drive
-Used by `examples/benchmark/benchmark.py` and reported in the README "Benchmarking"
-section and Hazan et al. 2018:
+Described in the README "Benchmarking" section and Hazan et al. 2018:
 - Population of **n** Poisson input neurons, firing rates drawn from **U(0, 100) Hz**.
 - Connected all-to-all to an equally sized population of LIF neurons; connection
   weights sampled from **N(0, 1)**.
 - **n** varied 250 → 10,000 in steps of 250; each run simulated **1,000 ms** at
   **dt = 1.0 ms**.
 
+`examples/benchmark/benchmark.py` (every simulator) and
+`examples/benchmark/bindsnet_scaling.py` (BindsNET only) implement exactly this, since
+2026-10; the 2018 script used its own input settings for each simulator (see
+[REPRODUCING.md](REPRODUCING.md#scaling-benchmark)).
+
 This stimulus is generated programmatically; there is no dataset to download.
 
 ---
diff --git a/Dockerfile b/Dockerfile
index f3be6f205..b4446b9e3 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -4,7 +4,7 @@
 #   docker run --rm -it --gpus all bindsnet bash    # --gpus needs the NVIDIA container toolkit
 # The torch wheels from the cu130 index bundle their CUDA libraries, so no CUDA base
 # image is needed; the host only needs an NVIDIA driver.
-FROM python:3.14-slim
+FROM python:3.13-slim
 
 ARG DEPS=dev
 ARG POETRY_VERSION=2.4.3
diff --git a/README.md b/README.md
index 1d8ed5c12..2e4cf289d 100644
--- a/README.md
+++ b/README.md
@@ -114,10 +114,16 @@ preprocessing — is declared in [DATA.md](https://github.com/BindsNET/bindsnet/
 ## Reproducing results
 
 [REPRODUCING.md](https://github.com/BindsNET/bindsnet/blob/master/REPRODUCING.md) maps each shipped model and published claim to its
-model class, example script, exact command, seed, and expected output (e.g. the
+model class, example script, exact command, seed, and measured output (e.g. the
 Diehl & Cook 2015 MNIST replication via `examples/mnist/eth_mnist.py`, and the
 Hazan et al. 2018 scaling benchmark).
 
+Each BindsNET release is an archived snapshot. BindsNET is a research package that
+keeps gaining features and fixes, and results can also depend on the hardware, so
+results may differ between releases. A published result is reproduced with the
+release it was produced on; please name that release in the publication, ideally by
+its version DOI (see [Citing the software](#citing-the-software)).
+
 ## Background
 
 The simulation of biologically plausible spiking neuron dynamics can be challenging. It is typically done by solving ordinary differential equations (ODEs) which describe said dynamics. PyTorch does not explicitly support the solution of differential equations (as opposed to [`brian2`](https://github.com/brian-team/brian2), for example), but we can convert the ODEs defining the dynamics into difference equations and solve them at regular, short intervals (a `dt` on the order of 1 millisecond) as an approximation. Of course, under the hood, packages like `brian2` are doing the same thing. Doing this in [`PyTorch`](https://pytorch.org/) is exciting for a few reasons:
@@ -137,6 +143,8 @@ We simulated a network with a population of n Poisson input neurons with firing
 
 Several packages, including BRIAN and PyNEST, allow the setting of certain global preferences; e.g., the number of CPU threads, the number of OpenMP processes, etc. We chose these settings for our benchmark study in an attempt to maximize each library's speed, but note that BindsNET requires no setting of such options. Our approach, inheriting the computational model of PyTorch, appears to make the best use of the available hardware, and therefore makes it simple for practitioners to get the best performance from their system with the least effort.
 
+The figure above is the 2018 result. `examples/benchmark/benchmark.py` runs exactly the setup described above for BindsNET, BRIAN2, PyNEST and ANNarchy in one command (BindsNET alone: `examples/benchmark/bindsnet_scaling.py`), and [REPRODUCING.md](REPRODUCING.md#scaling-benchmark) lists the times measured with it. It replaces the script written in 2018, which used its own input settings for each simulator.
+
 <p align="middle">
 <img src="https://raw.githubusercontent.com/BindsNET/bindsnet/master/docs/BindsNET%20benchmark.png" alt="BindsNET%20Benchmark"  width="503" height="403">
 </p>
diff --git a/REPRODUCING.md b/REPRODUCING.md
index 88c71ee4d..7e2b707ab 100644
--- a/REPRODUCING.md
+++ b/REPRODUCING.md
@@ -1,48 +1,136 @@
 # Reproducing results with BindsNET
 
 This table traces each model BindsNET describes or ships back to executable code: the
-model class, the example script, an exact command, the seed, the expected output, and
+model class, the example script, an exact command, the seed, the measured output, and
 the data it needs (declared in [DATA.md](DATA.md)).
 
-> **Honesty note.** Commands, defaults, seeds, and model classes below are verified
-> against the source. The **Expected output** cells describe *what the script reports*
-> and the qualitative trend; cells marked *(not measured here)* have **not** been run to
-> a final metric in producing this table — run the command to obtain the number for your
-> hardware. No accuracy/timing figure is asserted that was not measured.
+> **Measured, not quoted.** Every number in the **Measured output** column was
+> produced by the command next to it, at commit `f392a65`, on the machine described
+> under [Measurement setup](#measurement-setup). The runs are short (one pass over
+> 20,000 training images, 100 neurons), so they serve as reference points for checking
+> a reproduction; Diehl & Cook (2015) report up to 95% with 6,400 neurons.
+
+> **Which release.** Each BindsNET release is an archived snapshot. BindsNET is a
+> research package that keeps gaining features and fixes, and results can also depend
+> on the hardware, so results may differ between releases. A published result is
+> reproduced with the release it was produced on; please name that release in the
+> publication, ideally by its version DOI (see "Citing the software" in the README).
 
 ## Model → code → command map
 
-| Claim / source | Model class | Script | Command (defaults shown) | Seed | Expected output | Data |
-|----------------|-------------|--------|--------------------------|------|-----------------|------|
-| Diehl & Cook 2015 MNIST replication (DOI `10.3389/fncom.2015.00099`) | `DiehlAndCook2015` | `examples/mnist/eth_mnist.py` | `python examples/mnist/eth_mnist.py --n_neurons 100 --n_epochs 1 --time 250 --seed 0` | `--seed 0` (`torch.manual_seed`) | Prints test accuracy at end. **Measured:** all-activity **0.81**, proportion-weighting **0.82** at seed 0 with `--n_train 20000 --n_test 10000` (GPU, torch 2.6, ~7.8 h on an RTX 2070). Accuracy rises with `--n_neurons` and with the full 60000-sample train set (Diehl & Cook report up to ~95% at 6400 neurons). | MNIST |
-| Batched ETH MNIST | `DiehlAndCook2015` | `examples/mnist/batch_eth_mnist.py` | `python examples/mnist/batch_eth_mnist.py --n_neurons 100 --batch_size 32 --time 100 --seed 0` | `--seed 0` | Prints test accuracy; faster per-epoch via batching. *(not measured here)* | MNIST |
-| Supervised MNIST (label-clamped) | `DiehlAndCook2015` | `examples/mnist/supervised_mnist.py` | `python examples/mnist/supervised_mnist.py --n_neurons 100 --time 250 --intensity 32 --seed 0` | `--seed 0` | Prints test accuracy. *(not measured here)* | MNIST |
-| Convolutional SNN on MNIST | (in-script conv network) | `examples/mnist/conv_mnist.py` | `python examples/mnist/conv_mnist.py --time 50 --batch_size 1 --seed 0` | `--seed 0` | Prints accuracy. *(not measured here)* | MNIST |
-| Reservoir / liquid-state MNIST | (in-script reservoir) | `examples/mnist/reservoir.py` | `python examples/mnist/reservoir.py --n_neurons 500 --n_epochs 100 --time 250 --seed 0` | `--seed 0` | Prints accuracy after readout training. *(not measured here)* | MNIST |
-| Scaling benchmark (Hazan et al. 2018, DOI `10.3389/fninf.2018.00089`) | `Input` + `LIFNodes` via `Connection` | `examples/benchmark/benchmark.py` | **Not single-command** — see note below | n/a (timing study) | Runtime-vs-`n` curve; published figure is `docs/BindsNET benchmark.png` | synthetic Poisson drive (DATA.md) |
-| Atari Breakout (ANN→SNN demo) | trained ANN + SNN pipeline | `examples/breakout/play_breakout_from_ANN.py` | `python examples/breakout/play_breakout_from_ANN.py` | set in script | Plays Breakout from the shipped `trained_shallow_ANN.pt` | Atari Breakout (DATA.md) |
+Run the MNIST commands from `examples/mnist/` (the scripts read and download MNIST at
+`../../data/MNIST`). `CUDA_VISIBLE_DEVICES=` runs the script on the CPU; leave it out
+to use a GPU. Runs with the same seed on the same device repeat exactly, and CPU and
+GPU runs follow different learning trajectories (see [Determinism](#determinism)).
+
+| Claim / source | Model class | Script | Command | Seed | Measured output | Data |
+|----------------|-------------|--------|---------|------|-----------------|------|
+| Diehl & Cook 2015 MNIST replication (DOI `10.3389/fncom.2015.00099`) | `DiehlAndCook2015` | `examples/mnist/eth_mnist.py` | `CUDA_VISIBLE_DEVICES= python eth_mnist.py --n_neurons 100 --n_epochs 1 --time 250 --seed 0 --n_train 20000 --n_test 10000` | `--seed 0` | CPU. Test accuracy: all-activity **0.75**, proportion-weighting **0.76**. 2.06 h. | MNIST |
+| Batched ETH MNIST | `DiehlAndCook2015` | `examples/mnist/batch_eth_mnist.py` | `python batch_eth_mnist.py --n_neurons 100 --batch_size 32 --time 100 --seed 0 --n_train 20000 --n_test 10000` | `--seed 0` | GPU (RTX 2070). Test accuracy: all-activity **0.81**, proportion-weighting **0.83**. 7.6 min. | MNIST |
+| Supervised MNIST (label-clamped) | `DiehlAndCook2015` | `examples/mnist/supervised_mnist.py` | `CUDA_VISIBLE_DEVICES= python supervised_mnist.py --n_neurons 100 --time 250 --intensity 32 --seed 0 --n_train 20000 --n_test 10000` | `--seed 0` | CPU. Test accuracy: all-activity **0.64**, proportion-weighting **0.64**. 2.06 h. | MNIST |
+| Convolutional SNN on MNIST | (in-script conv network) | `examples/mnist/conv_mnist.py` | `CUDA_VISIBLE_DEVICES= python conv_mnist.py --time 50 --batch_size 1 --seed 0 --n_train 20000 --n_test 10000` | `--seed 0` | CPU. Test accuracy: all-activity **0.12**, proportion-weighting **0.14**. It learns early (0.27 / 0.28 after 1,000 images) and then declines; see [conv_mnist over training](#conv_mnist-over-training). 1.58 h. | MNIST |
+| Reservoir / liquid-state MNIST | (in-script reservoir) | `examples/mnist/reservoir.py` | `python reservoir.py --n_neurons 500 --n_epochs 100 --time 250 --seed 0` | `--seed 0` | CPU. Readout accuracy on 500 held-out test images: **53.00 %** (trained on 500 images, the script's default `--examples`). 5.8 min. | MNIST |
+| Scaling benchmark (Hazan et al. 2018, DOI `10.3389/fninf.2018.00089`) | `Input` + `LIFNodes` via `Connection` | `examples/benchmark/bindsnet_scaling.py`, `examples/benchmark/benchmark.py` | `python examples/benchmark/bindsnet_scaling.py --csv bindsnet_scaling.csv` | `--seed 0` | Simulation time for 1,000 ms; see [Scaling benchmark](#scaling-benchmark). | synthetic Poisson drive (DATA.md) |
+| Atari Breakout (ANN→SNN demo) | trained ANN + SNN pipeline | `examples/breakout/play_breakout_from_ANN.py` | `cd examples/breakout && python play_breakout_from_ANN.py` | `@@seed 42` (default) | GPU (RTX 2070). Reward per episode over the default 10 episodes: 5, 5, 9, 5, 5, 4, 5, 6, 4, 6 (mean **5.4**), identical in two runs. 73 min (the two runs ran at the same time). Writes `play-breakout_results.csv` in the working directory. | Atari Breakout (DATA.md) |
+
+Times are wall-clock times for training plus testing. The four CPU runs and the GPU run
+ran at the same time, each pinned to its own 12 CPU cores; a run with the machine to
+itself can be faster.
+
+## Measurement setup
+
+- Commit `f392a65` (`git worktree` of the branch, so no uncommitted change was used).
+- AMD Ryzen Threadripper 2990WX (32 cores), NVIDIA GeForce RTX 2070 (driver 615.71.09).
+- Python 3.11.9, PyTorch 2.14.0 (CUDA 13.0 build), the dependency versions in
+  `poetry.lock`.
+- MNIST as downloaded by torchvision.
 
 ## Notes
 
 ### Determinism
-- Each MNIST example accepts `--seed` (default `0`) and calls `torch.manual_seed(seed)`
-  and `torch.cuda.manual_seed_all(seed)`. Pass the same `--seed` to repeat a run.
-- Residual nondeterminism can come from CUDA atomic operations and first-run dataset
-  download ordering. For stricter determinism run on CPU and, where feasible, set
-  `torch.use_deterministic_algorithms(True)`.
-- An automated, seeded smoke-reproduction test
-  (`test/repro/test_smoke_repro.py`) runs a tiny network end-to-end on CPU and asserts
-  an exact pre-measured output, so determinism is checked continuously in CI.
-
-### Scaling benchmark is a multi-simulator study
-`examples/benchmark/benchmark.py` compares BindsNET against **BRIAN2, PyNEST, ANNarchy,
-BRIAN2genn, and Nengo**, and imports those packages plus an `experiments` helper module.
-It is therefore **not** a single-command reproduction: it requires those external
-simulators installed and the benchmark harness. The published BindsNET result is the
-figure `docs/BindsNET benchmark.png` and the parameters in the README "Benchmarking"
-section (Poisson inputs U(0,100) Hz, weights N(0,1), dt = 1.0 ms, 1000 ms/run, n from
-250 to 10,000). A BindsNET-only timing reproduction (no external simulators) can be
-built from `Input` + `LIFNodes` + `Connection`.
+- Each MNIST example accepts `--seed` (default `0`) and calls `torch.manual_seed(seed)`,
+  which seeds the CPU generator and every GPU (scripts that also draw from numpy seed
+  it too), so weight initialisation, Poisson encoding and data shuffling all repeat.
+  This holds since 2026-10; earlier releases seeded only the GPU generator on GPU
+  runs. The `eth_mnist.py` figure listed in earlier versions of this file (0.81 /
+  0.82) came from such a run, and the table above replaces it with a repeatable one.
+- Checked with the code of commit `8dd308c`: each of the five MNIST examples, run
+  twice with the same command on the GPU, printed identical accuracies (300 training
+  and 100 test images; `batch_eth_mnist.py` 640 and 320; `reservoir.py` 100 and 100).
+  `eth_mnist.py` and `supervised_mnist.py` also repeat exactly on the CPU, and
+  `play_breakout_from_ANN.py` gave the same reward in all 10 episodes in two GPU runs.
+  The input-to-excitatory weights of `DiehlAndCook2015` after 20 samples on the GPU
+  were bit-for-bit equal across two runs (`torch.equal`).
+- CPU and GPU runs with the same seed give different numbers (for example
+  `eth_mnist.py` with 300 training images: 0.30 on the CPU, 0.19 on the GPU):
+  floating-point sums round differently on the two devices, and learning carries the
+  difference forward. Report the device with a result.
+- `test/repro/test_smoke_repro.py` runs a tiny network end-to-end on CPU and asserts an
+  exact pre-measured output in CI. `test/repro/test_packaging_consistency.py` checks
+  that every example with `--seed` seeds torch unconditionally.
+
+### Scaling benchmark
+The README and the 2018 paper describe: n Poisson inputs with rates from U(0, 100) Hz,
+all-to-all to n LIF neurons, weights from N(0, 1), 1,000 ms at dt = 1.0 ms, n from 250
+to 10,000. `examples/benchmark/benchmark.py` runs exactly this setup in every simulator,
+with the LIF parameters of BindsNET's `LIFNodes` (rest and reset -65 mV, threshold
+-52 mV, membrane time constant 100 ms, refractory 5 ms; input spikes add the weight in
+mV to the membrane potential), each simulator in its own process.
+`examples/benchmark/bindsnet_scaling.py` runs BindsNET alone. They replace the script
+written for the 2018 paper (commit `6a1168d`, first in release 0.2.2), which used its
+own input settings for each simulator and a helper package that was not released, so
+the published figure (`docs/BindsNET benchmark.png`) stands as the 2018 result;
+details in `CHANGELOG.md`.
+
+Measured simulation time (`network.run` or the simulator's equivalent, excluding
+network construction) for 1,000 ms, same machine as above. BindsNET from
+`bindsnet_scaling.py` in the project environment; the other simulators from
+`benchmark.py` in a separate environment (BRIAN2 2.9.0, NEST 3.10.0, ANNarchy 5.0,
+PyTorch 2.14.1 CPU build, numpy 2.3.5), each run with the machine to itself.
+
+| Simulator (threads) | n = 1,000 | n = 5,000 | n = 10,000 |
+|---------------------|-----------|-----------|------------|
+| BindsNET GPU, RTX 2070 (`bindsnet_scaling.py`) | 0.43 s | 0.51 s | 1.27 s |
+| BindsNET CPU, PyTorch default 32 threads (`bindsnet_scaling.py`) | 0.32 s | 6.10 s | 17.24 s |
+| BindsNET CPU, PyTorch CPU build (`benchmark.py`) | 0.29 s | 5.02 s | 18.54 s |
+| BRIAN2, runtime mode, 1 thread | 0.41 s | 6.69 s | 33.51 s |
+| PyNEST, 64 threads | 3.71 s | 12.54 s | 22.82 s |
+| ANNarchy CPU, 1 thread (its default) | 1.33 s | 36.60 s | 163.54 s |
+| ANNarchy CPU, 64 threads (a test with `benchmark.py` changed to 64 threads, not the committed code) | 5.33 s | 149.30 s | 341.38 s |
+
+Each simulator ran with its default thread setting, except NEST with every core; none
+was tuned, so the table shows these settings on this machine rather than a ranking of
+the simulators. ANNarchy was also tried with all 64 threads, which was slower here and
+uneven between sizes (5,000: 149 s, 6,000: 58 s).
+
+brian2genn (needs the GeNN code generator) and ANNarchy on GPU (needs the CUDA
+compiler) are included in `benchmark.py` and were not run here. This run confirms that
+each simulator executes the network and reports its time; comparing spike statistics
+across simulators was outside its scope. Full results:
+`examples/benchmark/results_2026-10/`.
+
+### conv_mnist over training
+Same command with different training lengths (CPU, seed 0, 1,000 test images unless
+noted):
+
+| Training images | All-activity | Proportion-weighting |
+|-----------------|--------------|----------------------|
+| 1,000 | 0.27 | 0.28 |
+| 3,000 | 0.25 | 0.28 |
+| 10,000 | 0.15 | 0.17 |
+| 20,000 (10,000 test images) | 0.12 | 0.14 |
+
+Logging the convolutional layer during the first 2,750 training images located the
+cause: the mean adaptive threshold `theta` kept rising, from 19.4 mV to 68.8 mV, while
+recorded spikes per image went from 36.3 to 20.4 and the neurons that fired at all
+from 400 to 353. In `DiehlAndCookNodes.forward` (`bindsnet/network/nodes.py`) every
+neuron that crosses threshold is reset and raises `theta` before `one_spike` keeps a
+single spike per time step, so neurons whose spike is not kept also raise their
+threshold. Charging only the neuron that fires was measured with seeds 0, 1 and 2
+(20,000 training images, 10,000 test images, all-activity accuracy). `conv_mnist`
+averaged 0.17 against 0.15 with the current order, and `eth_mnist` averaged 0.74
+against 0.77, lower on every seed. This release keeps the current order; the
+`conv_mnist` threshold growth remains open.
 
 ### Data
 All datasets and synthetic stimuli these scripts use are declared in
diff --git a/bindsnet/environment/environment.py b/bindsnet/environment/environment.py
index d10c4ed1c..56beafeb6 100644
--- a/bindsnet/environment/environment.py
+++ b/bindsnet/environment/environment.py
@@ -9,6 +9,10 @@
 from bindsnet.datasets.preprocess import binary_image, crop, gray_scale, subsample
 from bindsnet.encoding import Encoder, NullEncoder
 
+# Gym ids that ale-py 0.9 removed, mapped to the id and settings they stood for
+# (ale-py 0.8.1, ale_py/gym.py: "Deterministic" = the base id with frameskip 4).
+REMOVED_GYM_IDS = {"BreakoutDeterministic-v4": ("Breakout-v4", {"frameskip": 4})}
+
 
 class Environment(ABC):
     # language=rst
@@ -89,7 +93,8 @@ def __init__(
             2D inputs.
         """
         self.name = name
-        self.env = gym.make(id=name, render_mode=render_mode)
+        gym_id, gym_kwargs = REMOVED_GYM_IDS.get(name, (name, {}))
+        self.env = gym.make(id=gym_id, render_mode=render_mode, **gym_kwargs)
         self.action_space = self.env.action_space
 
         self.encoder = encoder
diff --git a/examples/benchmark/annarchy.py b/examples/benchmark/annarchy.py
deleted file mode 100755
index c58763da6..000000000
--- a/examples/benchmark/annarchy.py
+++ /dev/null
@@ -1,150 +0,0 @@
-from __future__ import print_function
-
-import argparse
-import os
-from time import time as t
-
-import ANNarchy
-import numpy as np
-import pandas as pd
-
-plots_path = os.path.join("..", "..", "figures")
-benchmark_path = os.path.join("..", "..", "benchmark")
-if not os.path.isdir(benchmark_path):
-    os.makedirs(benchmark_path)
-
-
-def ANNarchy_cpu(n_neurons, time):
-    ANNarchy.setup(paradigm="openmp", dt=1.0)
-    ANNarchy.clear()
-
-    t1 = t()
-
-    IF = ANNarchy.Neuron(
-        parameters="""
-            tau_m = 10.0
-            tau_e = 5.0
-            vt = -54.0
-            vr = -60.0
-            El = -74.0
-            Ee = 0.0
-        """,
-        equations="""
-            tau_m * dv/dt = El - v + g_exc *  (Ee - vr) : init = -60.0
-            tau_e * dg_exc/dt = - g_exc
-        """,
-        spike="""
-            v > vt
-        """,
-        reset="""
-            v = vr
-        """,
-    )
-
-    Input = ANNarchy.PoissonPopulation(name="Input", geometry=n_neurons, rates=50.0)
-    Output = ANNarchy.Population(name="Output", geometry=n_neurons, neuron=IF)
-    proj = ANNarchy.Projection(pre=Input, post=Output, target="exc", synapse=None)
-    proj.connect_all_to_all(weights=ANNarchy.Uniform(0.0, 1.0))
-
-    ANNarchy.compile()
-    ANNarchy.simulate(duration=time)
-
-    return t() - t1
-
-
-def ANNarchy_gpu(n_neurons, time):
-    ANNarchy.setup(paradigm="cuda", dt=1.0)
-    ANNarchy.clear()
-
-    t1 = t()
-
-    IF = ANNarchy.Neuron(
-        parameters="""
-            tau_m = 10.0
-            tau_e = 5.0
-            vt = -54.0
-            vr = -60.0
-            El = -74.0
-            Ee = 0.0
-        """,
-        equations="""
-            tau_m * dv/dt = El - v + g_exc *  (Ee - vr) : init = -60.0
-            tau_e * dg_exc/dt = - g_exc
-        """,
-        spike="""
-            v > vt
-        """,
-        reset="""
-            v = vr
-        """,
-    )
-
-    Input = ANNarchy.PoissonPopulation(name="Input", geometry=n_neurons, rates=50.0)
-    Output = ANNarchy.Population(name="Output", geometry=n_neurons, neuron=IF)
-    proj = ANNarchy.Projection(pre=Input, post=Output, target="exc", synapse=None)
-    proj.connect_all_to_all(weights=ANNarchy.Uniform(0.0, 1.0))
-
-    ANNarchy.compile()
-    ANNarchy.simulate(duration=time)
-
-    return t() - t1
-
-
-def main(start=100, stop=1000, step=100, time=1000, interval=100, plot=False):
-    times = {"ANNarchy_cpu": []}
-
-    f = os.path.join(
-        benchmark_path, "benchmark_{start}_{stop}_{step}_{time}.csv".format(**locals())
-    )
-    if not os.path.isfile(f):
-        raise Exception("{0} not found.".format(f))
-
-    for n_neurons in range(start, stop + step, step):
-        print("\nRunning benchmark with {0} neurons.".format(n_neurons))
-        for framework in times.keys():
-            if framework == "ANNarchy_cpu" and n_neurons > 5000:
-                times[framework].append(np.nan)
-                continue
-
-            print("- {0}:".format(framework), end=" ")
-
-            fn = globals()[framework]
-            elapsed = fn(n_neurons=n_neurons, time=time)
-            times[framework].append(elapsed)
-
-            print("(elapsed: {0:.4f})".format(elapsed))
-
-    df = pd.read_csv(f, index_col=0)
-
-    for framework in times.keys():
-        print(pd.Series(times[framework]))
-        df[framework] = times[framework]
-
-    print()
-    print(df)
-    print()
-
-    df.to_csv(f)
-
-
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--start", type=int, default=100)
-    parser.add_argument("--stop", type=int, default=1000)
-    parser.add_argument("--step", type=int, default=100)
-    parser.add_argument("--time", type=int, default=1000)
-    parser.add_argument("--interval", type=int, default=100)
-    parser.add_argument("--plot", dest="plot", action="store_true")
-    parser.set_defaults(plot=False)
-    args = parser.parse_args()
-
-    print(args)
-
-    main(
-        start=args.start,
-        stop=args.stop,
-        step=args.step,
-        time=args.time,
-        interval=args.interval,
-        plot=args.plot,
-    )
diff --git a/examples/benchmark/benchmark.py b/examples/benchmark/benchmark.py
index 75f9bf92f..48894ce8a 100755
--- a/examples/benchmark/benchmark.py
+++ b/examples/benchmark/benchmark.py
@@ -1,251 +1,236 @@
+"""Scaling benchmark of Hazan et al. (2018), README "Benchmarking" section.
+
+Every simulator gets the same network: n Poisson input neurons with rates drawn
+from U(0, 100) Hz, connected all-to-all to n leaky integrate-and-fire neurons
+with weights drawn from N(0, 1) (in mV, added to the membrane potential on each
+input spike), simulated for 1000 ms at dt = 1.0 ms, for n = 250 to 10,000 in
+steps of 250. The LIF parameters are BindsNET's ``LIFNodes`` defaults: rest and
+reset -65 mV, threshold -52 mV, membrane time constant 100 ms, refractory 5 ms.
+
+    python examples/benchmark/benchmark.py --simulators bindsnet_cpu bindsnet_gpu brian2 nest
+    python examples/benchmark/plot_benchmark.py benchmark.csv
+
+Threads (not tuned): NEST uses every core; BindsNET uses PyTorch's default (one
+thread per physical core); ANNarchy uses its default of one thread (on a 64-thread
+Threadripper 2990WX, all 64 made it slower); BRIAN2 runs in its default runtime mode,
+which is single-threaded (its C++ standalone mode with OpenMP is not used here).
+
+Each simulator is imported only when requested; one that is not installed is
+skipped with a message. BindsNET needs nothing extra
+(``bindsnet_scaling.py`` holds its timing code). BindsNET, BRIAN2, NEST and ANNarchy
+(CPU) were run with this file in 2026-10 (``results_2026-10/``); the brian2genn and
+ANNarchy GPU paths are included and were not run. This file replaces the 2018
+script, which used its own input settings for each simulator (see REPRODUCING.md).
+"""
+
 import argparse
 import os
-from time import time as t
+import subprocess
+import sys
+import tempfile
+from time import perf_counter
 
-import brian2genn
-import matplotlib.pyplot as plt
-import nengo
 import numpy as np
-import pandas as pd
-import torch
-from brian2 import *
-from experiments import ROOT_DIR
-from experiments.benchmark import plot_benchmark
-from nest import *
-
-from bindsnet.encoding import poisson
-from bindsnet.network import Network
-from bindsnet.network.nodes import Input, LIFNodes
-from bindsnet.network.topology import Connection
-
-plots_path = os.path.join(ROOT_DIR, "figures")
-benchmark_path = os.path.join(ROOT_DIR, "benchmark")
-if not os.path.isdir(benchmark_path):
-    os.makedirs(benchmark_path)
-
-# "Warm up" the GPU.
-torch.set_default_tensor_type("torch.cuda.FloatTensor")
-x = torch.rand(1000)
-del x
-
-# BRIAN2 clock
-defaultclock = 1.0 * ms
-
-
-def BindsNET_cpu(n_neurons, time):
-    t0 = t()
-
-    torch.set_default_tensor_type("torch.FloatTensor")
-
-    t1 = t()
-
-    network = Network()
-    network.add_layer(Input(n=n_neurons), name="X")
-    network.add_layer(LIFNodes(n=n_neurons), name="Y")
-    network.add_connection(
-        Connection(source=network.layers["X"], target=network.layers["Y"]),
-        source="X",
-        target="Y",
-    )
-
-    data = {"X": poisson(datum=torch.rand(n_neurons), time=time)}
-    network.run(inputs=data, time=time)
-
-    return t() - t0, t() - t1
-
-
-def BindsNET_gpu(n_neurons, time):
-    if torch.cuda.is_available():
-        t0 = t()
-
-        torch.set_default_tensor_type("torch.cuda.FloatTensor")
+from bindsnet_scaling import run_once as bindsnet_run_once
 
-        t1 = t()
+REST, RESET, THRESH, TAU_M, REFRAC = -65.0, -65.0, -52.0, 100.0, 5.0  # mV, ms
 
-        network = Network()
-        network.add_layer(Input(n=n_neurons), name="X")
-        network.add_layer(LIFNodes(n=n_neurons), name="Y")
-        network.add_connection(
-            Connection(source=network.layers["X"], target=network.layers["Y"]),
-            source="X",
-            target="Y",
-        )
 
-        data = {"X": poisson(datum=torch.rand(n_neurons), time=time)}
-        network.run(inputs=data, time=time)
+def bindsnet_cpu(n, time, dt, rng):
+    return bindsnet_run_once(n, time, dt, "cpu")
 
-        return t() - t0, t() - t1
 
+def bindsnet_gpu(n, time, dt, rng):
+    import torch
 
-def BRIAN2(n_neurons, time):
-    t0 = t()
+    if not torch.cuda.is_available():
+        raise ImportError("CUDA is not available")
+    return bindsnet_run_once(n, time, dt, "cuda")
 
-    set_device("runtime", build_on_run=False)
-    defaultclock = 1.0 * ms
-    device.build()
 
-    t1 = t()
+def _brian2(n, time, dt, rng, genn):
+    import brian2 as b2
 
-    eqs_neurons = """
-        dv/dt = (ge * (-60 * mV) + (-74 * mV) - v) / (10 * ms) : volt
-        dge/dt = -ge / (5 * ms) : 1
-    """
+    t0 = perf_counter()
+    if genn:
+        import brian2genn  # noqa: F401  (registers the "genn" device)
 
-    input = PoissonGroup(n_neurons, rates=15 * Hz)
-    neurons = NeuronGroup(
-        n_neurons,
-        eqs_neurons,
-        threshold="v > (-54 * mV)",
-        reset="v = -60 * mV",
+        build_dir = tempfile.TemporaryDirectory(prefix="brian2genn_")
+        b2.set_device("genn", directory=build_dir.name)
+    b2.defaultclock.dt = dt * b2.ms
+    inputs = b2.PoissonGroup(n, rates=rng.uniform(0.0, 100.0, n) * b2.Hz)
+    lif = b2.NeuronGroup(
+        n,
+        f"dv/dt = ({REST}*mV - v) / ({TAU_M}*ms) : volt (unless refractory)",
+        threshold=f"v >= {THRESH}*mV",
+        reset=f"v = {RESET}*mV",
+        refractory=REFRAC * b2.ms,
         method="exact",
     )
-    S = Synapses(input, neurons, """w: 1""")
-    S.connect()
-    S.w = "rand() * 0.01"
-
-    run(time * ms)
-
-    return t() - t0, t() - t1
-
-
-def BRIAN2GENN(n_neurons, time):
-    t0 = t()
-
-    set_device("genn", build_on_run=False)
-    defaultclock = 1.0 * ms
-    device.build()
-
-    t1 = t()
-
-    eqs_neurons = """
-        dv/dt = (ge * (-60 * mV) + (-74 * mV) - v) / (10 * ms) : volt
-        dge/dt = -ge / (5 * ms) : 1
-    """
-
-    input = PoissonGroup(n_neurons, rates=15 * Hz)
-    neurons = NeuronGroup(
-        n_neurons,
-        eqs_neurons,
-        threshold="v > (-54 * mV)",
-        reset="v = -60 * mV",
-        method="exact",
+    lif.v = REST * b2.mV
+    syn = b2.Synapses(inputs, lif, "w : 1", on_pre="v_post += w*mV")
+    syn.connect()
+    syn.w = rng.standard_normal(n * n)
+    net = b2.Network(inputs, lif, syn)
+    t1 = perf_counter()
+    net.run(time * b2.ms, namespace={})
+    t2 = perf_counter()
+    if genn:
+        b2.device.reinit()
+        b2.device.activate()
+        b2.set_device("runtime")
+        build_dir.cleanup()
+    return t2 - t0, t2 - t1
+
+
+def brian2(n, time, dt, rng):
+    return _brian2(n, time, dt, rng, genn=False)
+
+
+def brian2genn(n, time, dt, rng):
+    return _brian2(n, time, dt, rng, genn=True)
+
+
+def nest(n, time, dt, rng):
+    import nest as nest_
+
+    t0 = perf_counter()
+    nest_.ResetKernel()
+    nest_.set_verbosity("M_ERROR")
+    nest_.resolution = dt
+    nest_.local_num_threads = os.cpu_count()
+    generators = nest_.Create(
+        "poisson_generator", n, params={"rate": rng.uniform(0.0, 100.0, n).tolist()}
     )
-    S = Synapses(input, neurons, """w: 1""")
-    S.connect()
-    S.w = "rand() * 0.01"
-
-    run(time * ms)
-
-    device.reinit()
-    device.activate()
-
-    return t() - t0, t() - t1
-
-
-def PyNEST(n_neurons, time):
-    t0 = t()
-
-    ResetKernel()
-    SetKernelStatus({"local_num_threads": 8, "resolution": 1.0})
-
-    t1 = t()
-
-    r_ex = 60.0  # [Hz] rate of exc. neurons
-
-    neuron = Create("iaf_psc_alpha", n_neurons)
-    noise = Create("poisson_generator", n_neurons)
-
-    SetStatus(noise, [{"rate": r_ex}])
-    Connect(noise, neuron)
-
-    Simulate(time)
-
-    return t() - t0, t() - t1
-
-
-def Nengo(n_neurons, time):
-    t0 = t()
-    t1 = t()
-
-    model = nengo.Network()
-    with model:
-        X = nengo.Ensemble(n_neurons, dimensions=1, neuron_type=nengo.LIF())
-        Y = nengo.Ensemble(n_neurons, dimensions=2, neuron_type=nengo.LIF())
-        nengo.Connection(X, Y, transform=np.random.rand(n_neurons, n_neurons))
-
-    with nengo.Simulator(model) as sim:
-        sim.run(time / 1000)
-
-    return t() - t0, t() - t1
-
-
-def main(start=100, stop=1000, step=100, time=1000, interval=100, plot=False):
-    f = os.path.join(benchmark_path, f"benchmark_{start}_{stop}_{step}_{time}.csv")
-    if os.path.isfile(f):
-        os.remove(f)
-
-    times = {
-        "BindsNET_cpu": [],
-        "BindsNET_gpu": [],
-        "BRIAN2": [],
-        "BRIAN2GENN": [],
-        "BRIAN2GENN comp.": [],
-        "PyNEST": [],  # , 'Nengo': []
-    }
-
-    for n_neurons in range(start, stop + step, step):
-        print(f"\nRunning benchmark with {n_neurons} neurons.")
-        for framework in times.keys():
-            if n_neurons > 2500 and framework == "PyNEST":
-                times[framework].append(np.nan)
-                continue
-
-            if framework == "BRIAN2GENN comp.":
-                continue
-
-            print(f"- {framework}:", end=" ")
-
-            fn = globals()[framework]
-            total, sim = fn(n_neurons=n_neurons, time=time)
-            times[framework].append(sim)
-
-            if framework == "BRIAN2GENN":
-                times["BRIAN2GENN comp."].append(total - sim)
+    # Parrot neurons repeat each generator's train, so every LIF neuron sees the
+    # same n input trains (a poisson_generator sends independent trains per target).
+    inputs = nest_.Create("parrot_neuron", n)
+    nest_.Connect(generators, inputs, "one_to_one")
+    lif = nest_.Create(
+        "iaf_psc_delta",
+        n,
+        params={
+            "E_L": REST,
+            "V_reset": RESET,
+            "V_th": THRESH,
+            "V_m": REST,
+            "tau_m": TAU_M,
+            "t_ref": REFRAC,
+        },
+    )
+    nest_.Connect(
+        inputs,
+        lif,
+        "all_to_all",
+        syn_spec={"weight": nest_.random.normal(0.0, 1.0), "delay": dt},
+    )
+    t1 = perf_counter()
+    nest_.Simulate(float(time))
+    t2 = perf_counter()
+    return t2 - t0, t2 - t1
+
+
+def _annarchy(n, time, dt, rng, paradigm):
+    import ANNarchy as ann
+
+    t0 = perf_counter()
+    ann.setup(paradigm=paradigm)
+    net = ann.Network(dt=dt)
+    neuron = ann.Neuron(
+        parameters=f"tau_m = {TAU_M}\nEl = {REST}\nvt = {THRESH}\nvr = {RESET}",
+        equations=f"tau_m * dv/dt = El - v : init = {REST}",
+        spike="v >= vt",
+        reset="v = vr",
+        refractory=REFRAC,
+    )
+    inputs = net.create(ann.PoissonPopulation(n, rates=0.0))
+    inputs.rates = rng.uniform(0.0, 100.0, n)
+    lif = net.create(geometry=n, neuron=neuron)
+    proj = net.connect(
+        inputs, lif, target="exc", synapse=ann.Synapse(pre_spike="v_post += w")
+    )
+    proj.all_to_all(weights=ann.Normal(0.0, 1.0))
+    with tempfile.TemporaryDirectory(prefix="annarchy_") as build_dir:
+        net.compile(directory=build_dir, silent=True)
+        t1 = perf_counter()
+        net.simulate(float(time))
+        t2 = perf_counter()
+    return t2 - t0, t2 - t1
 
-            print(f"(total: {total:.4f}; sim: {sim:.4f})")
 
-    print(times)
+def annarchy_cpu(n, time, dt, rng):
+    return _annarchy(n, time, dt, rng, "openmp")
 
-    df = pd.DataFrame.from_dict(times)
-    df.index = list(range(start, stop + step, step))
 
-    print()
-    print(df)
-    print()
+def annarchy_gpu(n, time, dt, rng):
+    return _annarchy(n, time, dt, rng, "cuda")
 
-    df.to_csv(f)
 
-    plot_benchmark.main(
-        start=start, stop=stop, step=step, time=time, interval=interval, plot=plot
+SIMULATORS = {
+    f.__name__: f
+    for f in (
+        bindsnet_cpu,
+        bindsnet_gpu,
+        brian2,
+        brian2genn,
+        nest,
+        annarchy_cpu,
+        annarchy_gpu,
     )
+}
 
 
-if __name__ == "__main__":
+def main():
     parser = argparse.ArgumentParser()
-    parser.add_argument("--start", type=int, default=100)
-    parser.add_argument("--stop", type=int, default=1000)
-    parser.add_argument("--step", type=int, default=100)
-    parser.add_argument("--time", type=int, default=1000)
-    parser.add_argument("--interval", type=int, default=100)
-    parser.add_argument("--plot", dest="plot", action="store_true")
-    parser.set_defaults(plot=False)
+    parser.add_argument("--start", type=int, default=250)
+    parser.add_argument("--stop", type=int, default=10000)
+    parser.add_argument("--step", type=int, default=250)
+    parser.add_argument("--time", type=int, default=1000, help="ms")
+    parser.add_argument("--dt", type=float, default=1.0, help="ms")
+    parser.add_argument(
+        "--simulators", nargs="+", default=list(SIMULATORS), choices=list(SIMULATORS)
+    )
+    parser.add_argument("--seed", type=int, default=0)
+    parser.add_argument("--csv", type=str, default="benchmark.csv")
+    parser.add_argument("--append", action="store_true", help=argparse.SUPPRESS)
     args = parser.parse_args()
 
-    main(
-        start=args.start,
-        stop=args.stop,
-        step=args.step,
-        time=args.time,
-        interval=args.interval,
-        plot=args.plot,
-    )
+    if not args.append:
+        with open(args.csv, "w") as f:
+            f.write("n,simulator,sim_s,total_s\n")
+        # One process per simulator: NEST 3.10 and ANNarchy 5.0 corrupt each
+        # other's memory when loaded into the same process.
+        common = [
+            f"--{k}={getattr(args, k)}"
+            for k in ("start", "stop", "step", "time", "dt", "seed", "csv")
+        ]
+        for name in args.simulators:
+            subprocess.run(
+                [sys.executable, __file__, *common, "--simulators", name, "--append"]
+            )
+        return
+
+    import torch
+
+    torch.manual_seed(args.seed)
+    rng = np.random.default_rng(args.seed)
+
+    available = []
+    for name in args.simulators:
+        try:
+            SIMULATORS[name](args.start, 10, args.dt, rng)  # import check + warm-up
+            available.append(name)
+        except Exception as e:  # not installed, no CUDA, or failed to build
+            print(f"skipping {name}: {e}")
+
+    with open(args.csv, "a") as f:
+        for n in range(args.start, args.stop + 1, args.step):
+            for name in available:
+                total, sim = SIMULATORS[name](n, args.time, args.dt, rng)
+                f.write(f"{n},{name},{sim:.6f},{total:.6f}\n")
+                f.flush()
+                print(f"{n:>6} {name:>13} sim {sim:9.4f} s  total {total:9.4f} s")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/benchmark/bindsnet_scaling.py b/examples/benchmark/bindsnet_scaling.py
new file mode 100644
index 000000000..27d55e1fb
--- /dev/null
+++ b/examples/benchmark/bindsnet_scaling.py
@@ -0,0 +1,85 @@
+"""BindsNET-only version of the scaling benchmark described in the README.
+
+n Poisson input neurons with rates drawn from U(0, 100) Hz, connected all-to-all
+to n LIF neurons with weights drawn from N(0, 1), simulated for 1000 ms at
+dt = 1.0 ms, for n = 250 to 10,000 in steps of 250. Needs no other simulator.
+
+    python examples/benchmark/bindsnet_scaling.py --csv bindsnet_scaling.csv
+
+Prints, per n and device, the wall time of ``network.run`` alone ("sim") and of
+building the network plus encoding plus running ("total").
+"""
+
+import argparse
+from time import perf_counter
+
+import torch
+
+from bindsnet.encoding import poisson
+from bindsnet.network import Network
+from bindsnet.network.nodes import Input, LIFNodes
+from bindsnet.network.topology import Connection
+
+
+def run_once(n, time, dt, device):
+    def sync():
+        if device == "cuda":
+            torch.cuda.synchronize()
+
+    sync()
+    t0 = perf_counter()
+    network = Network(dt=dt)
+    network.add_layer(Input(n=n), name="X")
+    network.add_layer(LIFNodes(n=n), name="Y")
+    network.add_connection(
+        Connection(network.layers["X"], network.layers["Y"], w=torch.randn(n, n)),
+        source="X",
+        target="Y",
+    )
+    network.to(device)
+    rates = torch.rand(n) * 100.0  # Hz
+    inputs = {"X": poisson(datum=rates, time=time, dt=dt, device=device)}
+    sync()
+    t1 = perf_counter()
+    network.run(inputs=inputs, time=time)
+    sync()
+    t2 = perf_counter()
+    return t2 - t0, t2 - t1
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--start", type=int, default=250)
+    parser.add_argument("--stop", type=int, default=10000)
+    parser.add_argument("--step", type=int, default=250)
+    parser.add_argument("--time", type=int, default=1000, help="ms")
+    parser.add_argument("--dt", type=float, default=1.0, help="ms")
+    parser.add_argument("--devices", nargs="+", default=["cpu", "cuda"])
+    parser.add_argument("--seed", type=int, default=0)
+    parser.add_argument("--csv", type=str, default=None)
+    args = parser.parse_args()
+
+    torch.manual_seed(args.seed)
+    devices = [d for d in args.devices if d != "cuda" or torch.cuda.is_available()]
+
+    # Warm up each device so the first size does not pay one-off start-up costs.
+    for device in devices:
+        run_once(args.start, 10, args.dt, device)
+
+    rows = []
+    print(f"{'n':>6} {'device':>6} {'sim_s':>9} {'total_s':>9}")
+    for n in range(args.start, args.stop + 1, args.step):
+        for device in devices:
+            total, sim = run_once(n, args.time, args.dt, device)
+            rows.append((n, device, sim, total))
+            print(f"{n:>6} {device:>6} {sim:>9.4f} {total:>9.4f}", flush=True)
+
+    if args.csv:
+        with open(args.csv, "w") as f:
+            f.write("n,device,sim_s,total_s\n")
+            for row in rows:
+                f.write("%d,%s,%.6f,%.6f\n" % row)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/benchmark/gpu_annarchy.py b/examples/benchmark/gpu_annarchy.py
deleted file mode 100755
index b65988a68..000000000
--- a/examples/benchmark/gpu_annarchy.py
+++ /dev/null
@@ -1,115 +0,0 @@
-from __future__ import print_function
-
-import argparse
-import os
-from time import time as t
-
-import ANNarchy
-import numpy as np
-import pandas as pd
-
-plots_path = os.path.join("..", "..", "figures")
-benchmark_path = os.path.join("..", "..", "benchmark")
-if not os.path.isdir(benchmark_path):
-    os.makedirs(benchmark_path)
-
-
-def ANNarchy_gpu(n_neurons, time):
-    t0 = t()
-
-    ANNarchy.setup(paradigm="cuda", dt=1.0)
-    ANNarchy.clear()
-
-    IF = ANNarchy.Neuron(
-        parameters="""
-            tau_m = 10.0
-            tau_e = 5.0
-            vt = -54.0
-            vr = -60.0
-            El = -74.0
-            Ee = 0.0
-        """,
-        equations="""
-            tau_m * dv/dt = El - v + g_exc *  (Ee - vr) : init = -60.0
-            tau_e * dg_exc/dt = - g_exc
-        """,
-        spike="""
-            v > vt
-        """,
-        reset="""
-            v = vr
-        """,
-    )
-
-    Input = ANNarchy.PoissonPopulation(name="Input", geometry=n_neurons, rates=50.0)
-    Output = ANNarchy.Population(name="Output", geometry=n_neurons, neuron=IF)
-    proj = ANNarchy.Projection(pre=Input, post=Output, target="exc", synapse=None)
-    proj.connect_all_to_all(weights=ANNarchy.Uniform(0.0, 1.0))
-
-    ANNarchy.compile()
-
-    t1 = t()
-
-    ANNarchy.simulate(duration=time)
-
-    return t() - t0, t() - t1
-
-
-def main(start=100, stop=1000, step=100, time=1000, interval=100, plot=False):
-    times = {"ANNarchy_gpu": [], "ANNarchy_gpu (w/ comp.)": []}
-
-    f = os.path.join(
-        benchmark_path, "benchmark_{start}_{stop}_{step}_{time}.csv".format(**locals())
-    )
-    if not os.path.isfile(f):
-        raise Exception("{0} not found.".format(f))
-
-    for n_neurons in range(start, stop + step, step):
-        print("\nRunning benchmark with {0} neurons.".format(n_neurons))
-        for framework in times.keys():
-            if "comp" in framework:
-                continue
-
-            print("- {0}:".format(framework), end=" ")
-
-            fn = globals()[framework]
-            total, sim = fn(n_neurons=n_neurons, time=time)
-            times[framework].append(sim)
-            times[framework + " (w/ comp.)"].append(total)
-
-            print("(total, sim: {0:.4f}, {1:.4f})".format(total, sim))
-
-    df = pd.read_csv(f, index_col=0)
-
-    for framework in times.keys():
-        print(pd.Series(times[framework]))
-        df[framework] = times[framework]
-
-    print()
-    print(df)
-    print()
-
-    df.to_csv(f)
-
-
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--start", type=int, default=100)
-    parser.add_argument("--stop", type=int, default=1000)
-    parser.add_argument("--step", type=int, default=100)
-    parser.add_argument("--time", type=int, default=1000)
-    parser.add_argument("--interval", type=int, default=100)
-    parser.add_argument("--plot", dest="plot", action="store_true")
-    parser.set_defaults(plot=False)
-    args = parser.parse_args()
-
-    print(args)
-
-    main(
-        start=args.start,
-        stop=args.stop,
-        step=args.step,
-        time=args.time,
-        interval=args.interval,
-        plot=args.plot,
-    )
diff --git a/examples/benchmark/plot_benchmark.py b/examples/benchmark/plot_benchmark.py
index 6975c645a..524ee850e 100755
--- a/examples/benchmark/plot_benchmark.py
+++ b/examples/benchmark/plot_benchmark.py
@@ -1,69 +1,36 @@
+"""Plots simulation time against n from the CSV that benchmark.py writes.
+
+python examples/benchmark/plot_benchmark.py benchmark.csv
+"""
+
 import argparse
 import os
 
 import matplotlib.pyplot as plt
 import pandas as pd
-from experiments import ROOT_DIR
-
-benchmark_path = os.path.join(ROOT_DIR, "benchmark")
-figure_path = os.path.join(ROOT_DIR, "figures")
-
-if not os.path.isdir(benchmark_path):
-    os.makedirs(benchmark_path)
 
 
-def main(start=100, stop=1000, step=100, time=1000, interval=100, plot=False):
-    name = f"benchmark_{start}_{stop}_{step}_{time}"
-    f = os.path.join(benchmark_path, name + ".csv")
-    df = pd.read_csv(f, index_col=0)
-
-    plt.plot(df["BindsNET_cpu"], label="BindsNET (CPU)", linestyle="-", color="b")
-    plt.plot(df["BindsNET_gpu"], label="BindsNET (GPU)", linestyle="-", color="g")
-    plt.plot(df["BRIAN2"], label="BRIAN2", linestyle="--", color="r")
-    plt.plot(df["BRIAN2GENN"], label="brian2genn", linestyle="--", color="c")
-    plt.plot(df["BRIAN2GENN comp."], label="brian2genn comp.", linestyle=":", color="c")
-    plt.plot(df["PyNEST"], label="PyNEST", linestyle="--", color="y")
-    plt.plot(df["ANNarchy_cpu"], label="ANNarchy (CPU)", linestyle="--", color="m")
-    plt.plot(df["ANNarchy_gpu"], label="ANNarchy (GPU)", linestyle="--", color="k")
-    plt.plot(
-        df["ANNarchy_gpu comp."], label="ANNarchy (GPU) comp.", linestyle=":", color="k"
-    )
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("csv", type=str)
+    parser.add_argument("--plot", dest="plot", action="store_true")
+    args = parser.parse_args()
 
-    # for c in df.columns:
-    #     if 'BindsNET' in c:
-    #         plt.plot(df[c], label=c, linestyle='-')
-    #     else:
-    #         plt.plot(df[c], label=c, linestyle='--')
+    df = pd.read_csv(args.csv)
+    for name, group in df.groupby("simulator", sort=False):
+        style = "-" if name.startswith("bindsnet") else "--"
+        plt.plot(group["n"], group["sim_s"], label=name, linestyle=style)
 
     plt.title("Benchmark comparison of SNN simulation libraries")
-    plt.xticks(range(0, stop + interval, interval))
     plt.xlabel("Number of input / output neurons")
     plt.ylabel("Simulation time (seconds)")
-    plt.legend(loc=1, prop={"size": 5})
+    plt.legend(loc="upper left", prop={"size": 7})
     plt.yscale("log")
+    plt.savefig(os.path.splitext(args.csv)[0] + ".png")
 
-    plt.savefig(os.path.join(figure_path, name + ".png"))
-
-    if plot:
+    if args.plot:
         plt.show()
 
 
 if __name__ == "__main__":
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--start", type=int, default=100)
-    parser.add_argument("--stop", type=int, default=1000)
-    parser.add_argument("--step", type=int, default=100)
-    parser.add_argument("--time", type=int, default=1000)
-    parser.add_argument("--interval", type=int, default=1000)
-    parser.add_argument("--plot", dest="plot", action="store_true")
-    parser.set_defaults(plot=False)
-    args = parser.parse_args()
-
-    main(
-        start=args.start,
-        stop=args.stop,
-        step=args.step,
-        time=args.time,
-        interval=args.interval,
-        plot=args.plot,
-    )
+    main()
diff --git a/examples/benchmark/results_2026-10/README.md b/examples/benchmark/results_2026-10/README.md
new file mode 100644
index 000000000..a1adf1936
--- /dev/null
+++ b/examples/benchmark/results_2026-10/README.md
@@ -0,0 +1,15 @@
+# Scaling benchmark results, 2026-10-18
+
+Measured on an AMD Ryzen Threadripper 2990WX (32 cores, 64 threads) with an NVIDIA
+GeForce RTX 2070 (driver 615.71.09). Simulation time for 1,000 ms at dt = 1.0 ms; the
+setup is in `../benchmark.py`. Summary and caveats: `REPRODUCING.md`, section
+"Scaling benchmark".
+
+| File | Produced by | Environment |
+|------|-------------|-------------|
+| `bindsnet_scaling.csv` | `python examples/benchmark/bindsnet_scaling.py` at commit `f392a65` | project environment: Python 3.11.9, PyTorch 2.14.0 (CUDA 13.0 build) |
+| `benchmark.csv` | `python examples/benchmark/benchmark.py --simulators bindsnet_cpu brian2 nest annarchy_cpu brian2genn annarchy_gpu --start 1000 --stop 10000 --step 1000` at commit `f392a65` (brian2genn and annarchy_gpu skipped: not installed / no CUDA compiler) | separate environment: BRIAN2 2.9.0, NEST 3.10.0, ANNarchy 5.0, PyTorch 2.14.1 CPU build, numpy 2.3.5 |
+| `annarchy_cpu_64_threads.csv` | same as `benchmark.csv`, `--simulators annarchy_cpu`, with ANNarchy set to 64 threads (a test; the committed `benchmark.py` uses ANNarchy's default of one thread) | as above |
+
+Columns: `sim_s` is the simulation call alone; `total_s` adds building the network and
+encoding the input (and, for ANNarchy, compiling it).
diff --git a/examples/benchmark/results_2026-10/annarchy_cpu_64_threads.csv b/examples/benchmark/results_2026-10/annarchy_cpu_64_threads.csv
new file mode 100644
index 000000000..aebdc3c20
--- /dev/null
+++ b/examples/benchmark/results_2026-10/annarchy_cpu_64_threads.csv
@@ -0,0 +1,11 @@
+n,simulator,sim_s,total_s
+1000,annarchy_cpu,5.331004,19.416848
+2000,annarchy_cpu,21.783171,37.006954
+3000,annarchy_cpu,102.750315,120.078945
+4000,annarchy_cpu,148.758482,170.783141
+5000,annarchy_cpu,149.296447,175.132875
+6000,annarchy_cpu,57.676196,87.164966
+7000,annarchy_cpu,181.293607,214.151751
+8000,annarchy_cpu,221.925256,263.761149
+9000,annarchy_cpu,267.979866,318.947949
+10000,annarchy_cpu,341.376518,398.998859
diff --git a/examples/benchmark/results_2026-10/benchmark.csv b/examples/benchmark/results_2026-10/benchmark.csv
new file mode 100644
index 000000000..5c6cd3710
--- /dev/null
+++ b/examples/benchmark/results_2026-10/benchmark.csv
@@ -0,0 +1,41 @@
+n,simulator,sim_s,total_s
+1000,bindsnet_cpu,0.293839,0.469469
+2000,bindsnet_cpu,1.217571,1.569120
+3000,bindsnet_cpu,2.118855,2.599682
+4000,bindsnet_cpu,3.346503,3.998598
+5000,bindsnet_cpu,5.017276,5.889279
+6000,bindsnet_cpu,7.054901,8.161672
+7000,bindsnet_cpu,8.928544,10.264814
+8000,bindsnet_cpu,11.902077,13.535372
+9000,bindsnet_cpu,15.019956,16.941857
+10000,bindsnet_cpu,18.540328,20.765416
+1000,brian2,0.414605,0.488819
+2000,brian2,1.550524,1.746375
+3000,brian2,2.907986,3.365023
+4000,brian2,4.499285,5.303062
+5000,brian2,6.685783,7.910205
+6000,brian2,11.051691,12.853773
+7000,brian2,16.616983,19.110844
+8000,brian2,21.976808,24.980678
+9000,brian2,30.919449,35.072334
+10000,brian2,33.514666,38.403776
+1000,nest,3.708811,4.285300
+2000,nest,5.457313,6.050769
+3000,nest,8.584372,9.368215
+4000,nest,10.410835,11.386339
+5000,nest,12.539116,13.673739
+6000,nest,14.460769,16.140607
+7000,nest,16.797107,18.531616
+8000,nest,18.638320,20.655046
+9000,nest,20.654695,23.099050
+10000,nest,22.821320,26.233835
+1000,annarchy_cpu,1.329493,15.798165
+2000,annarchy_cpu,6.152741,22.293137
+3000,annarchy_cpu,15.254862,33.528304
+4000,annarchy_cpu,20.634525,41.392665
+5000,annarchy_cpu,36.603881,62.004369
+6000,annarchy_cpu,52.836956,83.841308
+7000,annarchy_cpu,81.942182,120.693081
+8000,annarchy_cpu,116.262238,161.237963
+9000,annarchy_cpu,115.068515,166.466087
+10000,annarchy_cpu,163.542757,218.398253
diff --git a/examples/benchmark/results_2026-10/bindsnet_scaling.csv b/examples/benchmark/results_2026-10/bindsnet_scaling.csv
new file mode 100644
index 000000000..973da5281
--- /dev/null
+++ b/examples/benchmark/results_2026-10/bindsnet_scaling.csv
@@ -0,0 +1,81 @@
+n,device,sim_s,total_s
+250,cpu,0.178521,0.239002
+250,cuda,0.436185,0.440056
+500,cpu,0.197655,0.290911
+500,cuda,0.436627,0.442940
+750,cpu,0.240028,0.358870
+750,cuda,0.430795,0.441740
+1000,cpu,0.316971,0.478215
+1000,cuda,0.430571,0.444149
+1250,cpu,0.456749,0.644573
+1250,cuda,0.437293,0.463377
+1500,cpu,0.761022,0.964571
+1500,cuda,0.434337,0.459912
+1750,cpu,1.083341,1.322225
+1750,cuda,0.441847,0.488018
+2000,cpu,1.274368,1.562101
+2000,cuda,0.451153,0.500341
+2250,cpu,1.451540,1.792263
+2250,cuda,0.433404,0.487974
+2500,cpu,1.849701,2.210877
+2500,cuda,0.443770,0.504379
+2750,cpu,2.035141,2.458061
+2750,cuda,0.435267,0.527659
+3000,cpu,2.220441,2.678465
+3000,cuda,0.452950,0.560965
+3250,cpu,2.386480,2.903333
+3250,cuda,0.440806,0.565100
+3500,cpu,2.567510,3.123183
+3500,cuda,0.456886,0.600162
+3750,cpu,3.250204,3.851561
+3750,cuda,0.453197,0.580405
+4000,cpu,3.463382,4.111532
+4000,cuda,0.462108,0.684807
+4250,cpu,3.824384,4.517341
+4250,cuda,0.483596,0.690664
+4500,cpu,4.283830,5.034801
+4500,cuda,0.496825,0.730477
+4750,cpu,4.679010,5.478856
+4750,cuda,0.519962,0.779276
+5000,cpu,6.099752,6.953475
+5000,cuda,0.514456,0.794089
+5250,cpu,5.356654,6.266590
+5250,cuda,0.556842,0.880527
+5500,cpu,6.060059,7.173798
+5500,cuda,0.567678,0.910960
+5750,cpu,7.705449,8.748658
+5750,cuda,0.593767,0.967777
+6000,cpu,7.152771,8.253771
+6000,cuda,0.618340,1.020741
+6250,cpu,7.200378,8.363510
+6250,cuda,0.651328,1.088285
+6500,cpu,8.490247,9.725330
+6500,cuda,0.947270,1.418414
+6750,cpu,8.925832,10.225321
+6750,cuda,0.972485,1.603113
+7000,cpu,9.609491,10.976390
+7000,cuda,0.751130,1.290064
+7250,cpu,9.217090,10.648212
+7250,cuda,0.786495,1.368512
+7500,cpu,11.251869,12.744944
+7500,cuda,0.824260,1.447550
+7750,cpu,10.550629,12.107363
+7750,cuda,0.870547,1.533530
+8000,cpu,11.624152,13.286975
+8000,cuda,0.889405,1.591630
+8250,cpu,12.216137,14.095973
+8250,cuda,1.047930,1.797953
+8500,cpu,13.229201,15.058885
+8500,cuda,1.074910,1.891053
+8750,cpu,12.890625,14.812768
+8750,cuda,1.139088,1.985097
+9000,cpu,14.160773,16.147641
+9000,cuda,1.124698,2.057778
+9250,cpu,15.251417,17.317167
+9250,cuda,1.169328,2.149866
+9500,cpu,15.637174,17.734699
+9500,cuda,1.229931,2.512387
+9750,cpu,17.440537,19.653995
+9750,cuda,1.276026,2.359804
+10000,cpu,17.243479,19.557353
+10000,cuda,1.266757,2.409778
diff --git a/examples/breakout/README.md b/examples/breakout/README.md
index 17dac5a7d..cd0b8531b 100644
--- a/examples/breakout/README.md
+++ b/examples/breakout/README.md
@@ -17,7 +17,8 @@ see [../../DATA.md](../../DATA.md)).
 | Architecture | `nn.Linear(6400, 1000)` → `ReLU` → `nn.Linear(1000, 4)` (class `Net` in `play_breakout_from_ANN.py`) |
 | Input | 6400 features = a flattened 80×80 preprocessed Breakout frame |
 | Output | 4 units = the Breakout discrete action space |
-| Consumed by | `play_breakout_from_ANN.py:55` (`torch.load("trained_shallow_ANN.pt")`) |
+| Format | A `state_dict` (`fc1.weight`, `fc1.bias`, `fc2.weight`, `fc2.bias`), loaded with `weights_only=True`, PyTorch's safe default since 2.6 (it runs no pickled code). Re-saved in 2026-10 from the original pickled `Net` object, with identical tensors. |
+| Consumed by | `play_breakout_from_ANN.py` (`Net().load_state_dict(torch.load(..., weights_only=True))`) |
 | How it is used | Its `fc1`/`fc2` weights are transposed, scaled (`layer1scale=57.68`, `layer2scale=77.48`), and transplanted into a spiking network `Input(6400) → LIFNodes(1000) → LIFNodes(4)`, which is then run on Breakout through an `EnvironmentPipeline` with Poisson encoding — an ANN→SNN conversion demo. |
 
 ### Regeneration
@@ -25,5 +26,5 @@ see [../../DATA.md](../../DATA.md)).
 **The training script that produced `trained_shallow_ANN.pt` is not included in this
 repository.** The file is shipped as a pretrained weight blob. To regenerate it you would
 need to train a network with the `Net` architecture above (input 6400, hidden 1000,
-output 4) as a Breakout Q-network and save it with `torch.save`. If you reproduce or
+output 4) as a Breakout Q-network and save its `state_dict()` with `torch.save`. If you reproduce or
 replace this artifact, please document the training data, hyperparameters, and seed here.
diff --git a/examples/breakout/play_breakout_from_ANN.py b/examples/breakout/play_breakout_from_ANN.py
index 5b6be3e7d..1fb5cd4b2 100644
--- a/examples/breakout/play_breakout_from_ANN.py
+++ b/examples/breakout/play_breakout_from_ANN.py
@@ -51,8 +51,11 @@ def forward(self, x):
         return x
 
 
-# load ANN
-dqn_network = torch.load("trained_shallow_ANN.pt", map_location=device)
+# load ANN (a plain state dict, so weights_only loading runs no pickled code)
+dqn_network = Net()
+dqn_network.load_state_dict(
+    torch.load("trained_shallow_ANN.pt", map_location=device, weights_only=True)
+)
 
 # Build Spiking network.
 network = Network(dt=dt).to(device)
diff --git a/examples/breakout/random_network_baseline.py b/examples/breakout/random_network_baseline.py
index 31fe87706..6f12b6998 100644
--- a/examples/breakout/random_network_baseline.py
+++ b/examples/breakout/random_network_baseline.py
@@ -34,11 +34,10 @@
 print_interval = args.print_interval
 gpu = args.gpu
 
+# Seeds the CPU generator and every GPU.
+torch.manual_seed(seed)
 if gpu:
     torch.set_default_tensor_type("torch.cuda.FloatTensor")
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
 
 # Build network.
 network = Network(dt=dt)
diff --git a/examples/breakout/trained_shallow_ANN.pt b/examples/breakout/trained_shallow_ANN.pt
index 6ea56558c..69f12f796 100644
Binary files a/examples/breakout/trained_shallow_ANN.pt and b/examples/breakout/trained_shallow_ANN.pt differ
diff --git a/examples/dotTracing/dot_tracing.py b/examples/dotTracing/dot_tracing.py
index 94695a273..ec07e5eca 100644
--- a/examples/dotTracing/dot_tracing.py
+++ b/examples/dotTracing/dot_tracing.py
@@ -24,7 +24,9 @@
 # Handle arguments for dot tracing params.
 parser = argparse.ArgumentParser()
 parser.add_argument("--steps", type=int, default=100)
-parser.add_argument("--seed", type=int, default=0)
+parser.add_argument(
+    "--seed", type=int, default=1, help="0 seeds the dot movement from the clock"
+)
 parser.add_argument("--dim", type=int, default=28)
 parser.add_argument("--granularity", type=int, default=100)
 parser.add_argument("--neurons", type=int, default=100)
@@ -46,6 +48,10 @@
 
 args = parser.parse_args()
 
+# Seeds torch (weights, Poisson encoding, random actions) on CPU and GPU; the dot
+# simulator gets the same seed below.
+torch.manual_seed(args.seed)
+
 steps = args.steps  # timesteps in which the dot is moving
 dim = args.dim  # 28x28 square
 granularity = args.granularity  # granularity (or precision) of spike trains
@@ -263,6 +269,7 @@ def main():
         allow_stay=allow_stay,
         pandas=pandas,
         fpath=OUT_FILE_PATH,
+        seed=args.seed,
     )
     environment.reset()
 
diff --git a/examples/mnist/SOM_LM-SNNs.py b/examples/mnist/SOM_LM-SNNs.py
index c830ecba0..043ec4e9a 100644
--- a/examples/mnist/SOM_LM-SNNs.py
+++ b/examples/mnist/SOM_LM-SNNs.py
@@ -40,7 +40,7 @@
 parser.add_argument("--plot_interval", type=int, default=250)
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
-parser.set_defaults(plot=True, gpu=True)
+parser.set_defaults(plot=False, gpu=True)
 
 args = parser.parse_args()
 
@@ -62,14 +62,13 @@
 update_inhibation_weights = args.update_inhibation_weights
 
 # Sets up Gpu use
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling) and every GPU.
+torch.manual_seed(seed)
+
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
diff --git a/examples/mnist/batch_eth_mnist.py b/examples/mnist/batch_eth_mnist.py
index eb1374537..bb5ca9177 100644
--- a/examples/mnist/batch_eth_mnist.py
+++ b/examples/mnist/batch_eth_mnist.py
@@ -78,16 +78,14 @@
 update_steps = int(n_train / batch_size / n_updates)
 update_interval = update_steps * batch_size
 
-device = "cpu"
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling) and every GPU.
+torch.manual_seed(seed)
+
 # Sets up Gpu use
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
@@ -204,7 +202,7 @@
 
     pbar_training = tqdm(total=n_train)
     for step, batch in enumerate(train_dataloader):
-        if step * batch_size > n_train:
+        if step * batch_size >= n_train:
             break
 
         # Assign labels to excitatory neurons.
@@ -357,7 +355,7 @@
 pbar.set_description_str("Test progress: ")
 
 for step, batch in enumerate(test_dataloader):
-    if step * batch_size > n_test:
+    if step * batch_size >= n_test:
         break
     # Get next input sample.
     inputs = {"X": batch["encoded_image"]}
diff --git a/examples/mnist/conv1d_MNIST.py b/examples/mnist/conv1d_MNIST.py
index ca0446430..443874ff9 100644
--- a/examples/mnist/conv1d_MNIST.py
+++ b/examples/mnist/conv1d_MNIST.py
@@ -37,7 +37,7 @@
 parser.add_argument("--test", dest="train", action="store_false")
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
-parser.set_defaults(plot=True, gpu=True, train=True)
+parser.set_defaults(plot=False, gpu=True, train=True)
 
 args = parser.parse_args()
 
@@ -59,14 +59,13 @@
 plot = args.plot
 gpu = args.gpu
 
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling) and every GPU.
+torch.manual_seed(seed)
+
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
diff --git a/examples/mnist/conv3d_MNIST.py b/examples/mnist/conv3d_MNIST.py
index bfdb36844..a82b6adf7 100644
--- a/examples/mnist/conv3d_MNIST.py
+++ b/examples/mnist/conv3d_MNIST.py
@@ -38,7 +38,7 @@
 parser.add_argument("--test", dest="train", action="store_false")
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
-parser.set_defaults(plot=True, gpu=True, train=True)
+parser.set_defaults(plot=False, gpu=True, train=True)
 
 args = parser.parse_args()
 
@@ -60,14 +60,13 @@
 plot = args.plot
 gpu = args.gpu
 
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling) and every GPU.
+torch.manual_seed(seed)
+
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
diff --git a/examples/mnist/conv_mnist.py b/examples/mnist/conv_mnist.py
index e98fa4348..b169f2706 100644
--- a/examples/mnist/conv_mnist.py
+++ b/examples/mnist/conv_mnist.py
@@ -15,6 +15,7 @@
 )
 from bindsnet.datasets import MNIST
 from bindsnet.encoding import PoissonEncoder
+from bindsnet.evaluation import all_activity, assign_labels, proportion_weighting
 from bindsnet.learning import PostPre
 from bindsnet.network import Network
 from bindsnet.network.monitors import Monitor
@@ -42,7 +43,7 @@
 parser.add_argument("--test", dest="train", action="store_false")
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
-parser.set_defaults(plot=True, gpu=True, train=True)
+parser.set_defaults(plot=False, gpu=True, train=True)
 
 args = parser.parse_args()
 
@@ -64,14 +65,13 @@
 plot = args.plot
 gpu = args.gpu
 
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling) and every GPU.
+torch.manual_seed(seed)
+
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
@@ -80,14 +80,15 @@
     update_interval = n_test
 
 conv_size = int((28 - kernel_size + 2 * padding) / stride) + 1
-per_class = int((n_filters * conv_size * conv_size) / 10)
+n_neurons = n_filters * conv_size * conv_size
+n_classes = 10
 
 # Build network.
 network = Network()
 input_layer = Input(n=784, shape=(1, 28, 28), traces=True)
 
 conv_layer = DiehlAndCookNodes(
-    n=n_filters * conv_size * conv_size,
+    n=n_neurons,
     shape=(n_filters, conv_size, conv_size),
     traces=True,
 )
@@ -148,6 +149,14 @@
     voltages[layer] = Monitor(network.layers[layer], state_vars=["v"], time=time)
     network.add_monitor(voltages[layer], name="%s_voltages" % layer)
 
+# Label each output neuron by the class it fires most for, from the last
+# update_interval training samples (as in eth_mnist.py).
+spike_record = torch.zeros((update_interval, int(time / dt), n_neurons), device=device)
+labels = []
+assignments = -torch.ones(n_neurons, device=device)
+proportions = torch.zeros((n_neurons, n_classes), device=device)
+rates = torch.zeros((n_neurons, n_classes), device=device)
+
 # Train the network.
 print("Begin training.\n")
 start = t()
@@ -175,7 +184,7 @@
 
     for step, batch in enumerate(tqdm(train_dataloader)):
         # Get next input sample.
-        if step > n_train:
+        if step * batch_size >= n_train:
             break
         inputs = {"X": batch["encoded_image"].view(time, batch_size, 1, 28, 28)}
         if gpu:
@@ -185,6 +194,20 @@
         # Run the network on the input.
         network.run(inputs=inputs, time=time)
 
+        # Record output spikes per sample; reassign labels when the window is full.
+        s = spikes["Y"].get("s").view(int(time / dt), -1, n_neurons)
+        for b in range(s.size(1)):
+            if len(labels) == update_interval:
+                assignments, proportions, rates = assign_labels(
+                    spikes=spike_record,
+                    labels=torch.tensor(labels, device=device),
+                    n_labels=n_classes,
+                    rates=rates,
+                )
+                labels = []
+            spike_record[len(labels)] = s[:, b]
+            labels.append(int(label[b]))
+
         # Optionally plot various simulation information.
         if plot and batch_size == 1:
             image = batch["image"].view(28, 28)
@@ -212,3 +235,54 @@
 
 print("Progress: %d / %d (%.4f seconds)\n" % (n_epochs, n_epochs, t() - start))
 print("Training complete.\n")
+
+# Load MNIST test data.
+test_dataset = MNIST(
+    PoissonEncoder(time=time, dt=dt),
+    None,
+    "../../data/MNIST",
+    download=True,
+    train=False,
+    transform=transforms.Compose(
+        [transforms.ToTensor(), transforms.Lambda(lambda x: x * intensity)]
+    ),
+)
+test_dataloader = torch.utils.data.DataLoader(
+    test_dataset, batch_size=batch_size, shuffle=False, num_workers=0, pin_memory=gpu
+)
+
+# Test the network with learning switched off.
+print("Begin testing.\n")
+network.train(mode=False)
+accuracy = {"all": 0, "proportion": 0}
+n_seen = 0
+
+for step, batch in enumerate(tqdm(test_dataloader)):
+    if step * batch_size >= n_test:
+        break
+    inputs = {"X": batch["encoded_image"].view(time, -1, 1, 28, 28)}
+    if gpu:
+        inputs = {k: v.cuda() for k, v in inputs.items()}
+
+    network.run(inputs=inputs, time=time)
+
+    s = spikes["Y"].get("s").view(int(time / dt), -1, n_neurons).permute(1, 0, 2)
+    s = s.to(device)
+    label_tensor = batch["label"].to(device)
+    all_activity_pred = all_activity(
+        spikes=s, assignments=assignments, n_labels=n_classes
+    )
+    proportion_pred = proportion_weighting(
+        spikes=s, assignments=assignments, proportions=proportions, n_labels=n_classes
+    )
+    accuracy["all"] += float(torch.sum(label_tensor.long() == all_activity_pred).item())
+    accuracy["proportion"] += float(
+        torch.sum(label_tensor.long() == proportion_pred).item()
+    )
+    n_seen += label_tensor.numel()
+
+    network.reset_state_variables()
+
+print("\nAll activity accuracy: %.2f" % (accuracy["all"] / n_seen))
+print("Proportion weighting accuracy: %.2f \n" % (accuracy["proportion"] / n_seen))
+print("Testing complete.\n")
diff --git a/examples/mnist/eth_mnist.py b/examples/mnist/eth_mnist.py
index f5e6184b9..d4e6e8957 100644
--- a/examples/mnist/eth_mnist.py
+++ b/examples/mnist/eth_mnist.py
@@ -42,7 +42,7 @@
 parser.add_argument("--test", dest="train", action="store_false")
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
-parser.set_defaults(plot=True, gpu=True)
+parser.set_defaults(plot=False, gpu=True)
 
 args = parser.parse_args()
 
@@ -64,15 +64,14 @@
 plot = args.plot
 gpu = args.gpu
 
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling) and every GPU.
+torch.manual_seed(seed)
+
 # Sets up Gpu use
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
@@ -175,7 +174,7 @@
     )
 
     for step, batch in enumerate(tqdm(dataloader)):
-        if step > n_train:
+        if step >= n_train:
             break
         # Get next input sample.
         inputs = {"X": batch["encoded_image"].view(int(time / dt), 1, 1, 28, 28)}
diff --git a/examples/mnist/reservoir.py b/examples/mnist/reservoir.py
index 9946fc32d..8b4311e2f 100644
--- a/examples/mnist/reservoir.py
+++ b/examples/mnist/reservoir.py
@@ -37,7 +37,7 @@
 parser.add_argument("--update_interval", type=int, default=250)
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
-parser.set_defaults(plot=True, gpu=False, train=True)
+parser.set_defaults(plot=False, gpu=False, train=True)
 
 args = parser.parse_args()
 
@@ -55,19 +55,15 @@
 plot = args.plot
 gpu = args.gpu
 
+# Seeds numpy (thresholds), the CPU generator and every GPU.
 np.random.seed(seed)
-torch.cuda.manual_seed_all(seed)
 torch.manual_seed(seed)
 
 # Sets up Gpu use
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
 
@@ -128,7 +124,7 @@
 training_pairs = []
 pbar = tqdm(enumerate(dataloader))
 for i, dataPoint in pbar:
-    if i > n_iters:
+    if i >= n_iters:
         break
 
     # Extract & resize the MNIST samples image data for training
@@ -228,13 +224,27 @@ def forward(self, x):
         % (epoch + 1, n_epochs, avg_loss / len(training_pairs))
     )
 
-# Run same simulation on reservoir with testing data instead of training data
+# Run same simulation on reservoir with the held-out MNIST test set
 # (see training section for intuition)
+test_dataset = MNIST(
+    PoissonEncoder(time=time, dt=dt),
+    None,
+    root=os.path.join("..", "..", "data", "MNIST"),
+    download=True,
+    train=False,
+    transform=transforms.Compose(
+        [transforms.ToTensor(), transforms.Lambda(lambda x: x * intensity)]
+    ),
+)
+test_dataloader = torch.utils.data.DataLoader(
+    test_dataset, batch_size=1, shuffle=True, num_workers=0, pin_memory=gpu
+)
+
 n_iters = examples
 test_pairs = []
-pbar = tqdm(enumerate(dataloader))
+pbar = tqdm(enumerate(test_dataloader))
 for i, dataPoint in pbar:
-    if i > n_iters:
+    if i >= n_iters:
         break
     datum = dataPoint["encoded_image"].view(int(time / dt), 1, 1, 28, 28).to(device)
     label = dataPoint["label"]
@@ -279,5 +289,5 @@ def forward(self, x):
 
 print(
     "\n Accuracy of the model on %d test images: %.2f %%"
-    % (n_iters, 100 * correct / total)
+    % (total, 100 * correct / total)
 )
diff --git a/examples/mnist/supervised_mnist.py b/examples/mnist/supervised_mnist.py
index b394efa7a..2b8cf53a9 100644
--- a/examples/mnist/supervised_mnist.py
+++ b/examples/mnist/supervised_mnist.py
@@ -41,7 +41,7 @@
 parser.add_argument("--plot", dest="plot", action="store_true")
 parser.add_argument("--gpu", dest="gpu", action="store_true")
 parser.add_argument("--device_id", type=int, default=0)
-parser.set_defaults(plot=True, gpu=True, train=True)
+parser.set_defaults(plot=False, gpu=True, train=True)
 
 args = parser.parse_args()
 
@@ -63,15 +63,16 @@
 gpu = args.gpu
 device_id = args.device_id
 
+# Seeds the CPU generator (weight init, Poisson encoding, shuffling), every GPU,
+# and numpy (choice of the clamped neuron).
+torch.manual_seed(seed)
+np.random.seed(seed)
+
 # Sets up Gpu use
 device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-if gpu and torch.cuda.is_available():
-    torch.cuda.manual_seed_all(seed)
-else:
-    torch.manual_seed(seed)
+if not (gpu and torch.cuda.is_available()):
     device = "cpu"
-    if gpu:
-        gpu = False
+    gpu = False
 
 torch.set_num_threads(os.cpu_count() - 1)
 print("Running on Device = ", device)
@@ -155,7 +156,7 @@
 
 pbar = tqdm(total=n_train)
 for i, datum in enumerate(dataloader):
-    if i > n_train:
+    if i >= n_train:
         break
 
     image = datum["encoded_image"]
@@ -216,7 +217,9 @@
     # Optionally plot various simulation information.
     if plot:
         inpt = inputs["X"].view(time, 784).sum(0).view(28, 28)
-        input_exc_weights = network.connections[("X", "Ae")].w
+        input_exc_weights = (
+            network.connections[("X", "Ae")].feature_index["weight"].value
+        )
         square_weights = get_square_weights(
             input_exc_weights.view(784, n_neurons), n_sqrt, 28
         )
@@ -273,7 +276,7 @@
 
 pbar = tqdm(total=n_test)
 for step, batch in enumerate(test_dataset):
-    if step > n_test:
+    if step >= n_test:
         break
     # Get next input sample.
     inputs = {"X": batch["encoded_image"].view(int(time / dt), 1, 1, 28, 28)}
diff --git a/pyproject.toml b/pyproject.toml
index 964857ed8..1f65736ce 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -95,3 +95,7 @@ exclude = '''
   | logs
   )/
 '''
+
+[tool.pytest.ini_options]
+# Only the test suite; ignored folders such as logs/ can hold git worktrees.
+testpaths = ["test"]
diff --git a/test/environment/test_gym_environment.py b/test/environment/test_gym_environment.py
new file mode 100644
index 000000000..2840dec55
--- /dev/null
+++ b/test/environment/test_gym_environment.py
@@ -0,0 +1,15 @@
+from bindsnet.environment import GymEnvironment
+
+
+def test_breakout_deterministic_v4_still_resolves():
+    # ale-py 0.9+ no longer registers this id; the examples and the Breakout
+    # preprocessing branch still use it, so GymEnvironment maps it to
+    # Breakout-v4 with the fixed frameskip of 4 it used to mean.
+    env = GymEnvironment("BreakoutDeterministic-v4")
+    try:
+        assert env.env.unwrapped._frameskip == 4
+        env.reset()
+        obs, _, _, _ = env.step(1)
+        assert obs.shape[-2:] == (80, 80)
+    finally:
+        env.close()
diff --git a/test/repro/test_packaging_consistency.py b/test/repro/test_packaging_consistency.py
new file mode 100644
index 000000000..172da7320
--- /dev/null
+++ b/test/repro/test_packaging_consistency.py
@@ -0,0 +1,92 @@
+"""Static checks that the reproducibility artifacts agree with each other.
+
+Each check here guards a mismatch that once shipped: the Docker image on a
+Python that pyproject.toml does not support (#810), and examples whose
+``--seed`` left the CPU generator unseeded on GPU runs.
+"""
+
+import ast
+import re
+import tomllib
+from pathlib import Path
+
+from packaging.specifiers import SpecifierSet
+
+ROOT = Path(__file__).resolve().parents[2]
+
+
+def _requires_python():
+    with open(ROOT / "pyproject.toml", "rb") as f:
+        return SpecifierSet(tomllib.load(f)["project"]["requires-python"])
+
+
+def _supported_minors():
+    spec = _requires_python()
+    return {f"3.{m}" for m in range(0, 30) if f"3.{m}.0" in spec}
+
+
+def test_docker_python_is_supported():
+    dockerfile = (ROOT / "Dockerfile").read_text()
+    match = re.search(r"^FROM python:(\d+\.\d+)", dockerfile, re.M)
+    assert match, "Dockerfile base image is not python:X.Y"
+    assert f"{match.group(1)}.0" in _requires_python(), (
+        f"Dockerfile uses Python {match.group(1)}, outside requires-python "
+        f"{_requires_python()}"
+    )
+
+
+def test_ci_covers_every_supported_python():
+    workflow = (ROOT / ".github" / "workflows" / "python-app.yml").read_text()
+    tested = set(re.findall(r'"(3\.\d+)"', workflow))
+    assert tested == _supported_minors(), (
+        f"CI tests {sorted(tested)}, requires-python allows "
+        f"{sorted(_supported_minors())}"
+    )
+
+
+def test_citation_version_matches_package():
+    with open(ROOT / "pyproject.toml", "rb") as f:
+        version = tomllib.load(f)["project"]["version"]
+    citation = (ROOT / "CITATION.cff").read_text()
+    assert re.search(
+        rf"^version: {re.escape(version)}$", citation, re.M
+    ), f"CITATION.cff version does not match pyproject.toml ({version})"
+
+
+def _seeds_unconditionally(tree):
+    # True if some torch.manual_seed / torch.random.manual_seed call is not
+    # inside an if/else, a try, or a loop.
+    def visit(node, guarded):
+        if isinstance(node, ast.Call):
+            name = ast.unparse(node.func)
+            if (
+                name in ("torch.manual_seed", "torch.random.manual_seed")
+                and not guarded
+            ):
+                return True
+        guards = (ast.If, ast.Try, ast.For, ast.While, ast.IfExp)
+        child_guarded = guarded or isinstance(node, guards)
+        return any(visit(c, child_guarded) for c in ast.iter_child_nodes(node))
+
+    return visit(tree, False)
+
+
+def test_examples_seed_cpu_generator_unconditionally():
+    # torch.manual_seed seeds the CPU generator and every GPU. Called only in a
+    # CPU branch, GPU runs draw weights, Poisson spikes and shuffling unseeded.
+    missing = []
+    for path in sorted((ROOT / "examples").rglob("*.py")):
+        source = path.read_text()
+        if '"--seed"' not in source:
+            continue
+        if not _seeds_unconditionally(ast.parse(source)):
+            missing.append(str(path.relative_to(ROOT)))
+    assert not missing, f"--seed does not seed torch unconditionally in: {missing}"
+
+
+def test_publish_workflow_reads_project_version():
+    # The version lives in [project]; [tool.poetry] has none since 0.3.4, so a
+    # release's tag check that reads tool.poetry.version fails with KeyError.
+    workflow = (ROOT / ".github" / "workflows" / "publish.yml").read_text()
+    assert "['project']['version']" in workflow
+    assert "['tool']['poetry']['version']" not in workflow