P9 complete: compiled kernels validated (byte-identical, 78-80 GSA samples/min/core), scaling report

Graviton3 c7g.8xlarge, 360 GSE165845 arrays: 360/360 GTCs identical between the Rust
kernels and the pure-Python paths; 1 worker 78-80 samples/min (pure Python 29.6).
Scaling 1->8 workers efficiency 0.87-0.90; above that the 125 MB/s gp3 disk caps the
run (~620 samples/min) and, on a RAM disk, a fixed per-worker start-up cost limits a
360-sample batch to efficiency 0.48 at 32 workers (K-025). Worker state is pickled once.
docs/validation/P9.md, result JSONs, guide/README/CHANGELOG updates.

CHANGELOG.md

@@ -1,5 +1,13 @@
 # Changelog
 
+## Unreleased (P9)
+
+Optional compiled kernels `rust/simurg-core` (PyO3): fused GenCall calling and BAF/LRR, libm maps with CPython
+semantics, normalization nearest-neighbour searches, BPM parser. Byte-identical outputs to the pure-Python paths
+(`SIMURG_PURE_PYTHON=1`); 78–80 GSA samples/min per Graviton3 core (pure Python 30), manifest loading 4× faster.
+Address coverage memoised per IDAT layout; worker state pickled once. `bench/profile_genotype.py`,
+`bench/scaling.py`. Details and scaling limits: `docs/validation/P9.md`.
+
 ## 1.0.0 — 2026-09-28
 
 First release: the v1.0 goals of the master plan (§1.1), built and validated phase by phase (P0–P8).

README.md

@@ -22,8 +22,9 @@ differential-testing harness. Validation (`docs/validation/`):
 
 Commands: `simurg inspect`, `ingest`, `normalize`, `genotype`, `qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`,
 `gtc-to-vcf`, `export plink|bedgraph|finalreport|locus-summary`, `compare`. **User guide: `docs/user_guide.md`**;
-worked example: `examples/gsa_demo.sh`. Next (plan §18): P9 performance (compiled core), P10 packaging, then v2
+worked example: `examples/gsa_demo.sh`. P9 (performance) is done: optional Rust kernels, byte-identical outputs,
-modules (cluster training, CN/LOH, PGx).
+~80 GSA samples/min per server core (`docs/validation/P9.md`). Next (plan §18): P10 packaging, then v2 modules
+(cluster training, CN/LOH, PGx).
 
 ## Quick start
 

docs/known_unknowns.md

@@ -4,6 +4,7 @@ Plan §20 lists the open questions; this file records what the implementation ha
 
 | ID | Date | Topic | Finding | Status |
 |---|---|---|---|---|
+| K-025 | 2026-09-28 | Scaling limits (P9) | With the Rust kernels a Graviton3 core genotypes 78–80 GSA samples/min. At > 8 workers a 125 MB/s gp3 volume caps the run at ~620 samples/min (GTC ≈ 13 MB each); on a RAM disk efficiency is 0.77 at 16 and 0.48 at 32 workers for 360 samples, explained by a fixed start-up cost (≈ 2.5 s + 0.2 s per worker). Linear scaling for large batches is inferred from that model, not measured. | `docs/validation/P9.md`; measure a ≥ 5000-sample batch when one is at hand. |
 | K-024 | 2026-09-28 | GSE165845 (P8) | GSA-MD v1 chips have no public cluster file; the GSA v1 C1 EGT with the C2 build-38 BPM covers 99.996 % of bead addresses and gives median LogR deviation 0.246 (cluster mismatch, not noise). GSM5051309 (NA20362 rep 2) is a failed array (call rate 0.762, concordance 0.336 with its own replicates). 27 510 array records have no NYGC 30x record and are scored as hom-ref; 1444 of them are systematically called alt. | Recorded in `docs/validation/P8.md`; results reported with and without these. |
 | K-023 | 2026-09-28 | P6 optional outputs | Implemented and validated the same day: `--info-fields` (with `--egt`), `--extra-format`, `--merge-snv-indel`, `--auxiliary-loci`, PLINK 2 `.pgen`. Two decisions: OUT-035 rule (d) uses the MNV's last plus-strand base for minus-strand designs (spec amended); `--extra-format` writes NORMX/NORMY missing where the raw intensities are both 0, and THETA/R missing where both normalized intensities are 0 (148 values on 3 GSA samples; gtc2vcf writes 0). | resolved |
 | K-022 | 2026-09-28 | VCF oracles' limits (GSA v4 PGx, GRCh38) | gtc2vcf leaves 144 indel probes unresolved (`REF=I`/`D`) that SimurgArray resolves from `SourceSeq`; both leave 49 indels unresolved. GTCtoVCF (python-2.7 era) runs under Python 3 only with two local patches (pyvcf3 `Format` field; skip instead of crash on 3 CYP2D6 indels whose `SourceSeq` has the dash on the right, `[ALLELE/-]`, which SimurgArray handles), and needs the alt-contig loci filtered when the FASTA has no alts. Neither oracle applies ploidy (males' X/Y/MT hets stay diploid). | documented |

docs/validation/P9.md

@@ -0,0 +1,82 @@
+# Phase P9 validation — performance (compiled kernels)
+
+Date: 2026-09-28. Development and identity checks on the Raspberry Pi 5; scaling on AWS c7g.8xlarge (32 Graviton3
+cores, 64 GB) with GSE165845 (360 GSA-MD arrays, GSA-24v1-0 C2 manifest). Cloud cost ≈ $0.25 (three short spot runs,
+one of which stopped after 2 min because the image lacked a C linker).
+
+## What was built
+
+Profiling first (`bench/profile_genotype.py`): one GSA sample took 3.4 s on a Pi core, 63 % of it in calling — the
+bit-exact port evaluates `exp`/`log10`/`pow` element-wise through CPython's `math` (≈ 6.8 M scalar calls per sample) —
+and 24 % in the normalization fit's nearest-neighbour searches. Reading a manifest cost 13 s (GSA) to 53 s (Omni2.5)
+per command. The IDAT parser was already negligible (0.02 s), so it was not ported.
+
+`rust/simurg-core` (PyO3 + maturin, optional; imported as `simurg_core`, used through `simurg._core`):
+
+| Kernel | Replaces | Notes |
+|---|---|---|
+| `call_genotypes` | `genotype.gencall.call_genotypes` | one fused pass per locus: flair SDs, the six likelihood cases, posteriors, sub-scores, score map; every intermediate in the reference's f32/f64 type and order |
+| `baf_lrr` | `genotype.baf_lrr.baf_lrr` | canonical NaN 0xFFC00000 kept |
+| `exp`, `log10`, `log2`, `pow` | `gencall._libm` | glibc libm (the functions CPython's `math` calls) with CPython's error semantics (NaN passthrough, +NaN on domain error, `log10(0)` = −inf) |
+| `closest_binned`, `closest_exact` | the two nearest-neighbour searches of `normalization.fit` | float32 distances, first minimum wins, NaN wins as in `numpy.argmin` |
+| `bpm_locus_entries` | `formats.bpm._parse_locus_entries` | same columns, string interning and error messages |
+
+No floating-point contraction: Rust never fuses `a*b+c` (cf. K-012). Without the extension, or with
+`SIMURG_PURE_PYTHON=1`, the pure paths run unchanged. Also: address coverage is memoised per IDAT layout (it was
+recomputed per sample), and the worker state is pickled once instead of once per worker.
+
+## Acceptance (plan §18 P9)
+
+### Outputs identical to the Python reference
+
+| Check | Result |
+|---|---|
+| `tests/unit/test_core_kernels.py` (19 tests): every kernel vs its pure path on edge values (±0, ±inf, NaN payloads and signs, subnormals, overflow, ties, empty bins, NaN clusters, zero SDs, empty clusters, missing thresholds), bit patterns compared | pass (Pi; on Graviton3 the run log does not show the count, the 360-GTC identity below covers that machine) |
+| 36 GSA v4 demo samples (Pi) | 36/36 GTCs byte-identical; `sample_summary.tsv`, `controls_qc.tsv`, `relatedness.tsv` identical |
+| 6 Omni2.5 samples, `--preset autocall-1.6.2` (Pi) | 6/6 GTCs and all tables identical |
+| 360 GSE165845 samples (Graviton3) | **360/360 GTCs byte-identical**, sample summary identical |
+| 5 manifests (GSA v1 C1/C2, GSA v4, GDA v1, Omni2.5) parsed by both BPM parsers | every column identical |
+| Oracle tests P1, P3, P4 (vs gtc2vcf/vendor) | pass with the kernels |
+
+### Throughput: ≥ 20 samples/min/core (GSA-size arrays)
+
+| Machine, 1 worker | Pure Python | Kernels |
+|---|---|---|
+| Graviton3 (360 samples, per-sample phase) | 29.6 samples/min | **78–80 samples/min** (0.75 s per sample) |
+| Raspberry Pi 5 (`genotype_sample` only) | 3.4 s per sample | 1.1 s per sample (≈ 55/min) |
+
+Manifest loading: GSA 13 s → 3 s, Omni2.5 53 s → 13 s, GDA 40 s → 10 s (Pi). A 36-sample GSA run on the Pi
+(3 workers, with controls and cohort QC) went from 131 s to 45 s.
+
+### Scaling to 32 threads
+
+360 samples, `simurg genotype --no-cohort-qc`, per-sample phase from the run log (`bench/scaling.py`):
+
+| Workers | GTCs on the gp3 disk (125 MB/s) | | GTCs on a RAM disk | |
+|---|---|---|---|---|
+| | samples/min | efficiency | samples/min | efficiency |
+| 1 | 78 | 1.00 | 80 | 1.00 |
+| 4 | 290 | 0.93 | 305 | 0.95 |
+| 8 | 546 | 0.87 | 578 | 0.90 |
+| 16 | 618 | 0.49 | 979 | 0.77 |
+| 24 | — | — | 1186 | 0.62 |
+| 32 | 619 | 0.25 | 1226 | 0.48 |
+
+* On the default gp3 disk the run is **write-bound** above 8 workers: one GSA GTC is ~13 MB, so ~620 samples/min is
+  ~140 MB/s, the volume's throughput. The pure-Python path hit the same ceiling (550 samples/min at 32 workers).
+  Faster storage (gp3 with more MB/s, local NVMe, tmpfs) or fewer outputs is needed to go beyond it.
+* On a RAM disk the remaining gap is a **fixed start-up cost that grows with the worker count**: measured phase minus
+  ideal (0.75 s × 360 / workers) is 3.4 s at 4 workers, 3.6 s at 8, 5.2 s at 16, 7.0 s at 24 and 9.2 s at 32
+  (≈ 2.5 s + 0.2 s per worker: process spawn, imports, manifest context). It does not depend on the number of samples,
+  so with 360 samples (11 per worker at 32 workers) it dominates. **Inferred, not measured:** for a 10 000-sample batch
+  at 32 workers the same overhead gives an efficiency of ≈ 0.96. The plan's "linear scaling to 32 threads" is
+  therefore met for large batches on fast storage, and not for short batches or on a 125 MB/s disk.
+* Pickling the worker state once (instead of per worker) did not change start-up on the Pi (4–8 workers: 1.5–2.6 s
+  either way); it is kept because it removes serial work from the parent.
+
+## Notes
+
+* Remaining per-sample time with the kernels (Pi, 1.1 s): normalization fit 0.53 s (many small robust fits per bin),
+  calling 0.24 s, normalize/polar 0.08 s, raw X/Y 0.06 s. A full Rust normalization fit is the next step if needed.
+* Pi thread scaling is poor (36 samples, 3 workers ≈ 1.7× one core): the Pi 5's memory bandwidth, not the code.
+* CI (`.github/workflows/ci.yml`) builds the extension and runs the unit tests with and without it.

docs/validation/P9_identity_360.json

@@ -0,0 +1 @@
+{"gtcs": 360, "differ": 0, "sample_summary": "identical"}

docs/validation/P9_scaling_gp3.json

@@ -0,0 +1,74 @@
+{
+ "samples": 360,
+ "backend": "rust 1.0.0",
+ "cpu_count": 32,
+ "machine": "aarch64",
+ "runs": [
+  {
+   "threads": 1,
+   "rc": 0,
+   "wall_s": 299.1,
+   "setup_s": 6.9,
+   "per_sample_phase_s": 276.63,
+   "samples_per_min": 78.1,
+   "samples_per_min_per_core": 78.1,
+   "speedup_vs_first": 1.0,
+   "efficiency": 1.0
+  },
+  {
+   "threads": 2,
+   "rc": 0,
+   "wall_s": 166.5,
+   "setup_s": 7.1,
+   "per_sample_phase_s": 143.78,
+   "samples_per_min": 150.2,
+   "samples_per_min_per_core": 75.1,
+   "speedup_vs_first": 1.92,
+   "efficiency": 0.96
+  },
+  {
+   "threads": 4,
+   "rc": 0,
+   "wall_s": 96.9,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 74.58,
+   "samples_per_min": 289.6,
+   "samples_per_min_per_core": 72.4,
+   "speedup_vs_first": 3.71,
+   "efficiency": 0.927
+  },
+  {
+   "threads": 8,
+   "rc": 0,
+   "wall_s": 61.9,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 39.56,
+   "samples_per_min": 546.0,
+   "samples_per_min_per_core": 68.3,
+   "speedup_vs_first": 6.99,
+   "efficiency": 0.874
+  },
+  {
+   "threads": 16,
+   "rc": 0,
+   "wall_s": 57.4,
+   "setup_s": 6.9,
+   "per_sample_phase_s": 34.98,
+   "samples_per_min": 617.5,
+   "samples_per_min_per_core": 38.6,
+   "speedup_vs_first": 7.91,
+   "efficiency": 0.494
+  },
+  {
+   "threads": 32,
+   "rc": 0,
+   "wall_s": 57.4,
+   "setup_s": 6.9,
+   "per_sample_phase_s": 34.91,
+   "samples_per_min": 618.8,
+   "samples_per_min_per_core": 19.3,
+   "speedup_vs_first": 7.92,
+   "efficiency": 0.247
+  }
+ ]
+}

docs/validation/P9_scaling_ramdisk.json

@@ -0,0 +1,74 @@
+{
+ "samples": 360,
+ "backend": "rust 1.0.0",
+ "cpu_count": 32,
+ "machine": "aarch64",
+ "runs": [
+  {
+   "threads": 1,
+   "rc": 0,
+   "wall_s": 292.6,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 270.27,
+   "samples_per_min": 79.9,
+   "samples_per_min_per_core": 79.9,
+   "speedup_vs_first": 1.0,
+   "efficiency": 1.0
+  },
+  {
+   "threads": 4,
+   "rc": 0,
+   "wall_s": 93.3,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 70.94,
+   "samples_per_min": 304.5,
+   "samples_per_min_per_core": 76.1,
+   "speedup_vs_first": 3.81,
+   "efficiency": 0.953
+  },
+  {
+   "threads": 8,
+   "rc": 0,
+   "wall_s": 59.7,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 37.37,
+   "samples_per_min": 578.0,
+   "samples_per_min_per_core": 72.2,
+   "speedup_vs_first": 7.23,
+   "efficiency": 0.904
+  },
+  {
+   "threads": 16,
+   "rc": 0,
+   "wall_s": 44.4,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 22.06,
+   "samples_per_min": 979.3,
+   "samples_per_min_per_core": 61.2,
+   "speedup_vs_first": 12.25,
+   "efficiency": 0.766
+  },
+  {
+   "threads": 24,
+   "rc": 0,
+   "wall_s": 40.6,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 18.22,
+   "samples_per_min": 1185.6,
+   "samples_per_min_per_core": 49.4,
+   "speedup_vs_first": 14.83,
+   "efficiency": 0.618
+  },
+  {
+   "threads": 32,
+   "rc": 0,
+   "wall_s": 40.0,
+   "setup_s": 6.8,
+   "per_sample_phase_s": 17.62,
+   "samples_per_min": 1225.7,
+   "samples_per_min_per_core": 38.3,
+   "speedup_vs_first": 15.34,
+   "efficiency": 0.479
+  }
+ ]
+}

src/simurg/pipeline/genotype.py

@@ -9,6 +9,8 @@ import json
 import math
 import multiprocessing as mp
 import os
+import pickle
+import tempfile
 from concurrent.futures import ProcessPoolExecutor
 from dataclasses import dataclass, field, replace
 from pathlib import Path
@@ -147,6 +149,13 @@ def _init(
     )
 
 
+def _init_from_file(path: str) -> None:
+    """Worker initializer: the shared state is pickled once by the parent (not once per worker, which serialised
+    worker start-up: ~0.25 s per worker for a GSA manifest context) and read by every worker in parallel."""
+    with open(path, "rb") as fh:
+        _init(*pickle.load(fh))
+
+
 def _process(job: Job) -> dict[str, Any]:
     ctx: ManifestContext = _W["ctx"]
     params: GenotypeParams = _W["params"]
@@ -338,13 +347,17 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
             _init(ctx, params, auto_1_22, probes, started, beads, th.bg_offset)
             results = [_process(j) for j in jobs]
         else:
+            opts.out.mkdir(parents=True, exist_ok=True)
+            fd, state = tempfile.mkstemp(prefix=".simurg-workers-", suffix=".pkl", dir=opts.out)
+            try:
+                with os.fdopen(fd, "wb") as fh:
+                    pickle.dump((ctx, params, auto_1_22, probes, started, beads, th.bg_offset), fh, protocol=5)
                 with ProcessPoolExecutor(
-                workers,
+                    workers, mp_context=mp.get_context("spawn"), initializer=_init_from_file, initargs=(state,)
-                mp_context=mp.get_context("spawn"),
-                initializer=_init,
-                initargs=(ctx, params, auto_1_22, probes, started, beads, th.bg_offset),
                 ) as ex:
                     results = list(ex.map(_process, jobs))
+            finally:
+                Path(state).unlink(missing_ok=True)
     except KeyboardInterrupt:
         rm.status = "interrupted"
         rm.write(opts.out / "run_manifest.json")