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.