SimurgArray/docs/validation/P9.md

# 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.