From 4d5c9be4ae110d7ecc606b3310cc30b4d38926ab Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 22:32:48 +0100 Subject: [PATCH] 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 | 8 +++ README.md | 5 +- docs/known_unknowns.md | 1 + docs/validation/P9.md | 82 +++++++++++++++++++++++++ docs/validation/P9_identity_360.json | 1 + docs/validation/P9_scaling_gp3.json | 74 ++++++++++++++++++++++ docs/validation/P9_scaling_ramdisk.json | 74 ++++++++++++++++++++++ src/simurg/pipeline/genotype.py | 27 +++++--- 8 files changed, 263 insertions(+), 9 deletions(-) create mode 100644 docs/validation/P9.md create mode 100644 docs/validation/P9_identity_360.json create mode 100644 docs/validation/P9_scaling_gp3.json create mode 100644 docs/validation/P9_scaling_ramdisk.json diff --git a/CHANGELOG.md b/CHANGELOG.md index fe71354..62de98d 100644 --- a/CHANGELOG.md +++ b/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). diff --git a/README.md b/README.md index 134725c..12c743f 100644 --- a/README.md +++ b/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 -modules (cluster training, CN/LOH, PGx). +worked example: `examples/gsa_demo.sh`. P9 (performance) is done: optional Rust kernels, byte-identical outputs, +~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 diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 953c3ff..c69450b 100644 --- a/docs/known_unknowns.md +++ b/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 | diff --git a/docs/validation/P9.md b/docs/validation/P9.md new file mode 100644 index 0000000..652b5ef --- /dev/null +++ b/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. diff --git a/docs/validation/P9_identity_360.json b/docs/validation/P9_identity_360.json new file mode 100644 index 0000000..0be9878 --- /dev/null +++ b/docs/validation/P9_identity_360.json @@ -0,0 +1 @@ +{"gtcs": 360, "differ": 0, "sample_summary": "identical"} diff --git a/docs/validation/P9_scaling_gp3.json b/docs/validation/P9_scaling_gp3.json new file mode 100644 index 0000000..30330b0 --- /dev/null +++ b/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 + } + ] +} \ No newline at end of file diff --git a/docs/validation/P9_scaling_ramdisk.json b/docs/validation/P9_scaling_ramdisk.json new file mode 100644 index 0000000..d3aefdf --- /dev/null +++ b/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 + } + ] +} \ No newline at end of file diff --git a/src/simurg/pipeline/genotype.py b/src/simurg/pipeline/genotype.py index f9b2a4c..5a11be3 100644 --- a/src/simurg/pipeline/genotype.py +++ b/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: - with ProcessPoolExecutor( - workers, - 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)) + 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, mp_context=mp.get_context("spawn"), initializer=_init_from_file, initargs=(state,) + ) 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")