Release 1.1.0: compiled kernels (P9) and packaging (P10)

P10 complete: packaging verified on fresh x86_64 and aarch64 VMs
Wheels (pure + manylinux_2_28 kernels), sdist, Docker image and conda packages built and installed on fresh Ubuntu 24.04 machines; 225 unit tests pass against the installed package; the tutorial runs with pip (with and without kernels); GTC/VCF content digests are identical across pip, Docker and conda on both architectures and equal to the Pi's run. docs/validation/P10.md; guide: install from wheels, Docker, conda; exit code 141.
2026-09-28 23:56:39 +01:00 · 2026-09-28 23:53:59 +01:00 · 2026-09-28 23:14:47 +01:00 · 2026-09-28 22:41:56 +01:00 · 2026-09-28 22:32:48 +01:00 · 2026-09-28 22:04:10 +01:00
47 changed files with 2580 additions and 25 deletions
@@ -0,0 +1,14 @@
+.git
+.venv
+data
+refs
+dist
+results
+_phase0_archive
+**/__pycache__
+**/target
+*.idat
+*.idat.gz
+*.bpm
+*.egt
+*.gtc
@@ -23,3 +23,9 @@ jobs:
        run: uv run mypy
      - name: Unit tests (no downloaded data needed)
        run: uv run pytest tests/unit -q
+      - name: Build the compiled kernels (rust/simurg-core)
+        run: |
+          curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal
+          PATH="$HOME/.cargo/bin:$PATH" uv pip install ./rust/simurg-core
+      - name: Unit tests with the compiled kernels (must equal the pure-Python results)
+        run: uv run pytest tests/unit -q
@@ -1,5 +1,20 @@
 # Changelog

+## 1.1.0 — 2026-09-28 (P9, P10)
+
+**Packaging (P10).** `scripts/build_dist.sh` (sdist, pure wheel, manylinux_2_28 kernel wheels for x86_64 and
+aarch64), `Dockerfile`, conda recipes (`packaging/conda/`). Verified on fresh x86_64 and aarch64 VMs: pip, Docker
+and conda installs give identical outputs (`docs/validation/P10.md`). A closed stdout pipe now exits 141 silently
+instead of reporting an internal error.
+
+**Performance (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).
@@ -0,0 +1,28 @@
+# SimurgArray container: the `simurg` CLI with the compiled kernels.
+#   docker build -t simurg:1.1.0 .
+#   docker run --rm -v "$PWD":/work simurg:1.1.0 genotype --bpm M.bpm --egt C.egt --idats idats --out run/gtc
+# Inputs and outputs live in the mounted /work folder; the image holds no data.
+ARG PYTHON=3.12
+
+FROM python:${PYTHON}-slim-bookworm AS build
+RUN apt-get update -q && apt-get install -y -q --no-install-recommends curl gcc libc6-dev \
+    && rm -rf /var/lib/apt/lists/*
+RUN curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal
+ENV PATH=/root/.cargo/bin:$PATH
+WORKDIR /src
+COPY pyproject.toml README.md LICENSE THIRD_PARTY_NOTICES.md ./
+COPY src ./src
+COPY rust/simurg-core ./rust/simurg-core
+RUN pip wheel --no-cache-dir --no-deps -w /wheels . ./rust/simurg-core
+
+FROM python:${PYTHON}-slim-bookworm
+LABEL org.opencontainers.image.title="SimurgArray" \
+      org.opencontainers.image.description="Research toolkit for Illumina Infinium genotyping BeadChip data" \
+      org.opencontainers.image.licenses="MIT"
+COPY --from=build /wheels /wheels
+RUN pip install --no-cache-dir /wheels/*.whl && rm -rf /wheels \
+    && useradd --create-home --uid 1000 simurg
+USER simurg
+WORKDIR /work
+ENTRYPOINT ["simurg"]
+CMD ["--help"]
@@ -9,7 +9,7 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data

 ## Status

-**v1.0.0** (2026-09-28) — see `CHANGELOG.md`. IDAT → normalization → GenCall-equivalent genotype calling → GTC v5 →
+**v1.1.0** (2026-09-28) — see `CHANGELOG.md`. IDAT → normalization → GenCall-equivalent genotype calling → GTC v5 →
 sample/control QC with an HTML report → plus-strand VCF/BCF → PLINK 1/2, bedGraph, final report; plus a
 differential-testing harness. Validation (`docs/validation/`):

@@ -22,8 +22,10 @@ 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`). P10 (packaging) is done: wheels, Docker image and
+conda recipes, verified on fresh x86_64 and aarch64 machines (`docs/validation/P10.md`). Next (plan §18): v2 modules
+(cluster training, CN/LOH, PGx).

 ## Quick start

@@ -12,6 +12,7 @@ used only as external subprocess oracles and are not copied or linked.
 | [klib `ksort.h`](https://github.com/attractivechaos/klib) as shipped in htslib — Copyright (c) 2008, 2012-2013, 2017-2019 Genome Research Ltd | MIT | Quickselect `ks_ksmall` partition order, ported to `simurg/qc/contamination.py` (`ks_ksmall`) to reproduce BAFregress medians |
 | [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) — Copyright Illumina, Inc. | Apache-2.0 | Indel source-sequence handling ported to `simurg/harmonize/sites.py` (design-strand orientation, left shift, IUPAC flank matching, deletion/insertion context score); multi-probe consistency idea in `simurg/harmonize/combine.py`; test oracle (run under Python 3 with two local compatibility patches, never shipped) |
 | [Picard `picard.arrays`](https://github.com/broadinstitute/picard) | MIT | Second format reference |
+| [PyO3](https://github.com/PyO3/pyo3), [rust-numpy](https://github.com/PyO3/rust-numpy), [ndarray](https://github.com/rust-ndarray/ndarray) and their Rust dependencies | MIT / Apache-2.0 / BSD-2-Clause | Compiled into the optional `simurg_core` extension (`rust/simurg-core`); list with `cargo tree` |
 | [illuminaio](https://github.com/HenrikBengtsson/illuminaio) | GPL-2 | Documentation of IDAT field codes only; oracle via subprocess |

 ## BSD-2-Clause (BeadArrayFiles, verbatim from its README)
@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+"""Per-stage timing of `simurg genotype` for one sample (plan §19 measurements).
+
+python bench/profile_genotype.py BPM EGT GRN_IDAT [--repeat N]
+"""
+
+from __future__ import annotations
+
+import argparse
+import datetime
+import time
+from dataclasses import replace
+from pathlib import Path
+
+import numpy as np
+
+from simurg.formats.bpm import read_bpm
+from simurg.formats.egt import read_egt
+from simurg.formats.gtc import write_gtc
+from simurg.formats.idat import read_idat
+from simurg.genotype.baf_lrr import baf_lrr, logr_dev
+from simurg.genotype.clusters import align_clusters
+from simurg.genotype.gencall import call_genotypes
+from simurg.genotype.sample import genotype_sample, manifest_context, params_for_preset
+from simurg.genotype.sex import estimate_sex
+from simurg.signal.normalization.apply import normalize, polar
+from simurg.signal.normalization.fit import fit_transforms
+from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("bpm")
+    ap.add_argument("egt")
+    ap.add_argument("grn")
+    ap.add_argument("--repeat", type=int, default=3)
+    ap.add_argument("--out", default="/tmp/profile.gtc")
+    a = ap.parse_args()
+    m = read_bpm(a.bpm)
+    ctx = manifest_context(m, align_clusters(read_egt(a.egt), m))
+    lite = m.lite()
+    ctx = replace(ctx, manifest=lite, cache=LocusMapCache(lite))
+    params = params_for_preset("current")
+    red_path = a.grn.replace("_Grn", "_Red")
+    t: dict[str, list[float]] = {}
+
+    def lap(name: str, t0: float) -> float:
+        now = time.perf_counter()
+        t.setdefault(name, []).append(now - t0)
+        return now
+
+    for _ in range(a.repeat):
+        t0 = time.perf_counter()
+        grn, red = read_idat(a.grn), read_idat(red_path)
+        t0 = lap("read_idat x2", t0)
+        raw = assemble_raw_xy(grn, red, lite, cache=ctx.cache)
+        t0 = lap("assemble_raw_xy", t0)
+        xf = fit_transforms(raw.raw_x, raw.raw_y, ctx.norm_bins, params.norm_version)
+        t0 = lap("fit_transforms", t0)
+        nx, ny = normalize(raw.raw_x, raw.raw_y, ctx.lookup, xf)
+        th, r = polar(nx, ny)
+        t0 = lap("normalize+polar", t0)
+        z = (raw.raw_x == 0) & (raw.raw_y == 0)
+        thc = np.where(z, np.float32(np.nan), th).astype(np.float32)
+        rc = np.where(z, np.float32(np.nan), r).astype(np.float32)
+        calls = call_genotypes(thc, rc, z, ctx.clusters, lite.intensity_only, params.cutoff)
+        t0 = lap("call_genotypes", t0)
+        estimate_sex(calls.genotypes, r, ctx.clusters.gentrain, ctx.clusters.present, ctx.sex_classes, params.sex)
+        t0 = lap("estimate_sex", t0)
+        _, lrr = baf_lrr(thc, rc, z, ctx.clusters, lite.intensity_only)
+        logr_dev(lrr, ctx.autosomal, lite.intensity_only, ctx.clusters.present)
+        t0 = lap("baf_lrr+logr_dev", t0)
+        g = genotype_sample(grn, red, ctx, params, datetime.datetime(2026, 1, 1))
+        t0 = lap("genotype_sample (whole)", t0)
+        write_gtc(g, Path(a.out))
+        t0 = lap("write_gtc", t0)
+    print(f"{'stage':28s} {'median s':>9s}")
+    for k, v in t.items():
+        print(f"{k:28s} {sorted(v)[len(v) // 2]:9.3f}")
+
+
+if __name__ == "__main__":
+    main()
@@ -0,0 +1,114 @@
+#!/usr/bin/env python3
+"""Thread scaling of ``simurg genotype`` (plan §18 P9 gate: ≥ 20 samples/min/core on GSA-size arrays, linear
+scaling to 32 threads).
+
+For every worker count the same N samples are genotyped (``--no-cohort-qc``); the per-sample phase is timed from
+the run log ("genotyping … sample(s)" → "done"), so manifest loading is reported separately and does not blur
+the scaling. Writes a JSON summary.
+
+    python bench/scaling.py --bpm M.bpm --egt C.egt --idats DIR --out OUT --threads 1,2,4,8 [--samples N]
+"""
+
+from __future__ import annotations
+
+import argparse
+import datetime as dt
+import json
+import os
+import shutil
+import subprocess
+import sys
+import time
+from pathlib import Path
+
+
+def phase_seconds(log: Path) -> tuple[float, float]:
+    """(manifest+setup seconds, per-sample phase seconds) from a genotype run log."""
+    ts: dict[str, dt.datetime] = {}
+    for line in log.read_text().splitlines():
+        rec = json.loads(line)
+        msg = rec.get("msg") or ""
+        t = dt.datetime.fromisoformat(rec["ts"].replace("Z", "+00:00"))
+        ts.setdefault("first", t)
+        if msg.startswith("genotyping "):
+            ts["start"] = t
+        elif msg.startswith("done:"):
+            ts["done"] = t
+    return (ts["start"] - ts["first"]).total_seconds(), (ts["done"] - ts["start"]).total_seconds()
+
+
+def subset(idats: Path, n: int, dest: Path) -> Path:
+    grn = sorted(p for p in idats.rglob("*_Grn.idat*"))
+    if n <= 0 or n >= len(grn):
+        return idats
+    dest.mkdir(parents=True, exist_ok=True)
+    for g in grn[:n]:
+        d = dest / g.parent.name
+        d.mkdir(exist_ok=True)
+        for p in (g, g.with_name(g.name.replace("_Grn.", "_Red."))):
+            link = d / p.name
+            if not link.exists():
+                link.symlink_to(p.resolve())
+    return dest
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--bpm", required=True)
+    ap.add_argument("--egt", required=True)
+    ap.add_argument("--idats", required=True, type=Path)
+    ap.add_argument("--out", required=True, type=Path)
+    ap.add_argument("--threads", default="1,2,4")
+    ap.add_argument("--samples", type=int, default=0, help="use the first N samples (0 = all)")
+    ap.add_argument("--simurg", default=str(Path(sys.executable).with_name("simurg")))
+    ap.add_argument("--pure-python", action="store_true", help="run with SIMURG_PURE_PYTHON=1")
+    ap.add_argument("--work", type=Path, help="folder for the GTC outputs (e.g. a tmpfs); default: --out")
+    ap.add_argument("--clean", action="store_true", help="delete each run's GTCs after timing it")
+    a = ap.parse_args()
+    a.out.mkdir(parents=True, exist_ok=True)
+    idats = subset(a.idats, a.samples, a.out / "idats")
+    n = len(list(idats.rglob("*_Grn.idat*")))
+    env = {**os.environ, "SIMURG_PURE_PYTHON": "1" if a.pure_python else "0"}
+    runs = []
+    for t in (int(x) for x in a.threads.split(",")):
+        out = (a.work or a.out) / f"t{t}"
+        t0 = time.time()
+        cmd = [a.simurg, "genotype", "--bpm", a.bpm, "--egt", a.egt, "--idats", str(idats), "--out", str(out),
+               "--threads", str(t), "--no-cohort-qc", "--autocall-date", "1/1/2026 12:00:00 AM",
+               "--quiet", "--force"]  # fmt: skip
+        rc = subprocess.run(cmd, env=env).returncode
+        wall = time.time() - t0
+        setup, phase = phase_seconds(out / "logs" / "genotype.log")
+        if a.clean:
+            shutil.rmtree(out / "samples", ignore_errors=True)
+        runs.append(
+            {
+                "threads": t,
+                "rc": rc,
+                "wall_s": round(wall, 1),
+                "setup_s": round(setup, 1),
+                "per_sample_phase_s": round(phase, 2),
+                "samples_per_min": round(60 * n / phase, 1) if phase else None,
+                "samples_per_min_per_core": round(60 * n / phase / t, 1) if phase else None,
+            }
+        )
+        print(json.dumps(runs[-1]), flush=True)
+    base = runs[0]["per_sample_phase_s"] * runs[0]["threads"]
+    for r in runs:
+        r["speedup_vs_first"] = round(base / r["per_sample_phase_s"], 2) if r["per_sample_phase_s"] else None
+        r["efficiency"] = round(r["speedup_vs_first"] / r["threads"], 3) if r["speedup_vs_first"] else None
+    from simurg._core import backend
+
+    summary = {
+        "samples": n,
+        "backend": "python" if a.pure_python else backend(),
+        "cpu_count": os.cpu_count(),
+        "machine": os.uname().machine,
+        "runs": runs,
+    }
+    (a.out / "scaling.json").write_text(json.dumps(summary, indent=1))
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
@@ -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 |
@@ -1,4 +1,4 @@
-# SimurgArray user guide (v1.0)
+# SimurgArray user guide (v1.1)

 SimurgArray turns Illumina Infinium genotyping BeadChip scans (IDAT files) into genotype calls, QC tables and
 reports, VCF/BCF and PLINK files. It is a personal, non-commercial **research** tool — not for clinical or
@@ -19,6 +19,40 @@ uv venv && uv pip install -e .          # add ".[dev]" for tests and linters
 .venv/bin/simurg --version
 ```

+**Optional compiled kernels** (≈ 3–4× faster genotyping, 4× faster manifest loading; results are byte-identical):
+
+```bash
+curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal   # Rust (needs a C linker, e.g. gcc)
+uv pip install ./rust/simurg-core
+.venv/bin/python -c "from simurg._core import backend; print(backend())"                   # rust 1.1.0
+```
+
+Without them the same code runs in pure Python/NumPy; `SIMURG_PURE_PYTHON=1` forces that path.
+
+**From release files** (`bash scripts/build_dist.sh` writes them to `dist/`: a pure wheel, an sdist and
+`manylinux_2_28` kernel wheels for x86_64 and aarch64, glibc ≥ 2.28):
+
+```bash
+python3 -m venv simurg-env && simurg-env/bin/pip install simurg-1.1.0-py3-none-any.whl \
+    simurg_core-1.1.0-cp312-abi3-manylinux_2_28_$(uname -m).whl        # the second wheel is optional
+```
+
+**Docker** (image ≈ 550 MB, kernels included, runs as a non-root user; data stay in the mounted folder):
+
+```bash
+docker build -t simurg:1.1.0 .
+docker run --rm --user "$(id -u):$(id -g)" -v "$PWD":/work simurg:1.1.0 \
+    genotype --bpm M.bpm --egt C.egt --idats idats --out run/gtc
+```
+
+**conda** (recipes in `packaging/conda/`):
+
+```bash
+conda build -c conda-forge --output-folder conda-bld packaging/conda/simurg-core packaging/conda/simurg
+conda create -n simurg -c ./conda-bld -c conda-forge simurg simurg-core
+```
+
+All four installs give identical results on x86_64 and aarch64 (`docs/validation/P10.md`).
 Only pure-Python dependencies are needed (numpy, pyarrow, typer, rich, pyyaml, openpyxl). VCF/BCF, tabix/CSI
 indexes and PLINK files are written by SimurgArray itself; bcftools/plink2 are not required (they are used only
 as test oracles).
@@ -148,8 +182,9 @@ A nightly differential run (`scripts/nightly.py`) re-checks these against the or
  diplotyping are not in v1.0.
 * The largest cohort tested is 2141 samples (plan goal ≥ 10 000). Memory grows with samples × loci in the cohort
  QC steps; on an 8 GB machine keep `--threads` low (≈ 2–3 GB per worker for GSA/Omni-size arrays).
-* Speed is that of NumPy: genotyping a GSA sample takes ~3 CPU-seconds on a server core (360 samples in 2 min on
-  16 cores) and ~25 s on a Raspberry Pi 5 core; a compiled core is planned (P9).
+* Speed: with the compiled kernels a GSA sample takes ~1.1 s on a Raspberry Pi 5 core (pure Python ~3.4 s); see
+  `docs/validation/P9.md` for server numbers and thread scaling. Loading a manifest/cluster file is a fixed cost
+  per run (GSA ~5 s, Omni2.5 ~15 s on the Pi).
 * Methylation arrays: control-probe QC only.

 ## Troubleshooting
@@ -157,5 +192,6 @@ A nightly differential run (`scripts/nightly.py`) re-checks these against the or
 Errors and warnings carry a stable code (`E_…`, `W_…`) — all are listed in `docs/messages.md` with their spec item.
 Add `--debug` for a traceback or `--log-json` for machine-readable logs. Exit codes (CLI-014): 0 success,
 1 internal error, 2 usage or pre-flight validation error (nothing processed), 3 some samples failed (or
-`compare --fail-on-discrepancy` found discrepancies), 4 no sample succeeded, 130/143 interrupted. Findings about vendor behaviour that SimurgArray reproduces
+`compare --fail-on-discrepancy` found discrepancies), 4 no sample succeeded, 130/143 interrupted, 141 the reader of
+stdout went away (e.g. `simurg inspect … | head`; nothing is reported). Findings about vendor behaviour that SimurgArray reproduces
 (or deliberately does not) are in `docs/known_unknowns.md`.
@@ -0,0 +1,43 @@
+# Phase P10 validation — packaging and fresh-machine acceptance
+
+Date: 2026-09-28. Artifacts built and tested on fresh AWS Ubuntu 24.04 VMs (4 vCPU): aarch64 c7g.xlarge and x86_64
+m6a.xlarge / c6a.xlarge (one x86 spot instance was reclaimed by AWS after 4 min and relaunched). Cloud cost ≈ $0.10.
+
+## Artifacts
+
+| Artifact | How | Notes |
+|---|---|---|
+| `simurg-1.0.0-py3-none-any.whl`, `simurg-1.0.0.tar.gz` | `uv build` (hatchling) | QC report assets, profiles, `LICENSE` and `THIRD_PARTY_NOTICES.md` inside |
+| `simurg_core-1.0.0-cp312-abi3-manylinux_2_28_{aarch64,x86_64}.whl` | `maturin --zig --compatibility manylinux_2_28` (both built on either machine) | one wheel per architecture for every Python ≥ 3.12; glibc ≥ 2.28 |
+| Docker image (`Dockerfile`) | multi-stage, `python:3.12-slim-bookworm`, non-root user | 545 MB (x86_64), 551 MB (aarch64) |
+| conda packages (`packaging/conda/`) | `conda build -c conda-forge` | `simurg` (noarch) + `simurg-core` (Rust compiler from conda-forge) |
+
+`scripts/build_dist.sh` builds the wheels and sdist and writes `dist/SHA256SUMS`. A wheel built on the Pi is not
+portable (glibc 2.41 symbols such as `exp@GLIBC_2.29`), hence the zig-linked manylinux build.
+
+## Acceptance (plan §18 P10: clean install + tutorial works on a fresh VM)
+
+`scripts/cloud/tasks/p10.sh` (and `p10t.sh`, the same without Docker/conda) on each architecture:
+
+| Step | aarch64 | x86_64 |
+|---|---|---|
+| Build sdist, pure wheel, both kernel wheels | ok | ok |
+| `pip install` into a new venv of the system Python 3.12.3 (with / without kernels) | ok / ok | ok / ok |
+| Unit tests against the installed package (225) | pass | pass |
+| Tutorial `examples/gsa_demo.sh` on 6 demo samples, pip install with kernels | ok (506 s) | ok (521 s) |
+| Tutorial, pure-Python install | ok (539 s) | ok (494 s) |
+| Docker build + `genotype` + `gtc-to-vcf` | ok | ok |
+| conda build + install + `genotype` | ok | ok |
+
+Content digests (`scripts/p10/digest.py`: every numeric GTC field and sample statistic with NaNs canonicalised, all
+VCF data lines) of the 6 GTCs and 6 VCFs are **identical across all of them**: pip with kernels, pip pure, Docker and
+conda, on aarch64 and on x86_64, and equal to the Raspberry Pi's own run of the same tutorial. The x86_64 results come
+from AMD EPYC (m6a/c6a) machines; bit-exactness had so far only been shown on ARM.
+
+## Found and fixed
+
+* `simurg inspect … | head` made the CLI report `E_INTERNAL: BrokenPipeError` and exit 1, which stopped the tutorial
+  under `set -o pipefail` (it passed on the Pi only by timing luck). A closed stdout now exits with 141 (like a tool
+  killed by SIGPIPE) and prints nothing; the tutorial reads with `sed` instead of `head`; unit test added.
+* The `fast` extra (`simurg[fast]` → `simurg-core`) was dropped: `uv` must resolve every extra, and `simurg-core` is
+  not on a package index. The kernels are installed as their own wheel.
@@ -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.
@@ -0,0 +1 @@
+{"gtcs": 360, "differ": 0, "sample_summary": "identical"}
@@ -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
+  }
+ ]
+}
@@ -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
+  }
+ ]
+}
@@ -30,7 +30,7 @@ mkdir -p "$OUT"
 # Optional subset: symlink the first N IDAT pairs into $OUT/idats (same <barcode>/<barcode>_<pos>_<ch>.idat layout).
 if [[ $N -gt 0 ]]; then
  rm -rf "$OUT/idats"
-  find "$IDATS" -name '*_Grn.idat' | sort | head -n "$N" | while read -r g; do
+  find "$IDATS" -name '*_Grn.idat' | sort | sed -n "1,${N}p" | while read -r g; do
    b=$(basename "$(dirname "$g")"); mkdir -p "$OUT/idats/$b"
    ln -s "$(realpath "$g")" "$OUT/idats/$b/"; ln -s "$(realpath "${g%_Grn.idat}_Red.idat")" "$OUT/idats/$b/"
  done
@@ -38,12 +38,12 @@ if [[ $N -gt 0 ]]; then
 fi

 echo "== 1. inspect the inputs"
-$SIMURG inspect "$BPM" "$EGT" | head -20
+$SIMURG inspect "$BPM" "$EGT" | sed -n 1,20p   # sed reads to the end: no SIGPIPE under pipefail

 echo "== 2. genotype: GTC v5 per sample, sample_summary.tsv, control-probe QC, contamination, relatedness"
 $SIMURG genotype --bpm "$BPM" --egt "$EGT" --csv "$CSV" --idats "$IDATS" --out "$OUT/gtc" \
  --qc-profile genotyping --threads "$THREADS" --force
-column -t -s$'\t' "$OUT/gtc/sample_summary.tsv" | cut -c1-150 | head -8
+column -t -s$'\t' "$OUT/gtc/sample_summary.tsv" | cut -c1-150 | sed -n 1,8p

 echo "== 3. HTML QC report + QC table"
 $SIMURG qc report --data "$OUT/gtc" --out "$OUT/qc_report" --output-format both
@@ -0,0 +1,40 @@
+# conda recipe for the optional compiled kernels of SimurgArray (Rust, PyO3; abi3 for Python >= 3.12).
+#   conda build packaging/conda/simurg-core -c conda-forge
+{% set version = load_file_regex(load_file="Cargo.toml", regex_pattern='version = "([^"]+)"', from_recipe_dir=False).group(1) %}
+
+package:
+  name: simurg-core
+  version: {{ version }}
+
+source:
+  path: ../../../rust/simurg-core
+
+build:
+  number: 0
+  skip: true  # [win]
+  script: {{ PYTHON }} -m pip install . --no-deps -vv
+
+requirements:
+  build:
+    - {{ compiler('c') }}
+    - {{ compiler('rust') }}
+  host:
+    - python >=3.12
+    - pip
+    - maturin >=1.7,<2
+  run:
+    - python >=3.12
+    - numpy >=1.26
+
+test:
+  imports:
+    - simurg_core
+  commands:
+    - python -c "import simurg_core, numpy as np; assert simurg_core.log10(np.array([0.0]))[0] == -np.inf"
+
+about:
+  summary: Optional compiled kernels of SimurgArray (results identical to the pure-Python paths)
+  license: MIT
+  license_file:
+    - LICENSE
+    - NOTICE
@@ -0,0 +1,48 @@
+# conda recipe for SimurgArray (pure Python). Build from a checkout:
+#   conda build packaging/conda/simurg -c conda-forge
+{% set version = load_file_regex(load_file="src/simurg/__init__.py", regex_pattern='__version__ = "([^"]+)"', from_recipe_dir=False).group(1) %}
+
+package:
+  name: simurg
+  version: {{ version }}
+
+source:
+  path: ../../..
+
+build:
+  noarch: python
+  number: 0
+  script: {{ PYTHON }} -m pip install . --no-deps --no-build-isolation -vv
+  entry_points:
+    - simurg = simurg.cli.main:app
+
+requirements:
+  host:
+    - python >=3.12
+    - pip
+    - hatchling >=1.27
+  run:
+    - python >=3.12
+    - numpy >=1.26
+    - typer >=0.12
+    - rich >=13
+    - pyyaml >=6
+    - pyarrow >=15
+    - openpyxl >=3.1
+  run_constrained:
+    - simurg-core =={{ version }}
+
+test:
+  imports:
+    - simurg
+    - simurg.cli.main
+  commands:
+    - simurg --version
+    - simurg --help
+
+about:
+  summary: Research toolkit for Illumina Infinium genotyping BeadChip data (IDAT to GTC, QC, VCF, PLINK)
+  license: MIT
+  license_file:
+    - LICENSE
+    - THIRD_PARTY_NOTICES.md
@@ -1,14 +1,25 @@
 [build-system]
-requires = ["hatchling>=1.24"]
+requires = ["hatchling>=1.27"]
 build-backend = "hatchling.build"

 [project]
 name = "simurg"
 description = "SimurgArray: open research toolkit for Illumina Infinium genotyping BeadChip data"
 readme = "README.md"
-license = { text = "MIT" }
+license = "MIT"
+license-files = ["LICENSE", "THIRD_PARTY_NOTICES.md"]
 requires-python = ">=3.12"
 authors = [{ name = "Kemal Yaylali" }]
+keywords = ["genotyping", "microarray", "infinium", "idat", "gtc", "vcf", "bioinformatics"]
+classifiers = [
+    "Development Status :: 5 - Production/Stable",
+    "Intended Audience :: Science/Research",
+    "Operating System :: POSIX :: Linux",
+    "Operating System :: MacOS",
+    "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+    "Topic :: Scientific/Engineering :: Bio-Informatics",
+]
 dynamic = ["version"]
 dependencies = [
    "numpy>=1.26",   # all numerics
@@ -0,0 +1 @@
+target/
@@ -0,0 +1,228 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "autocfg"
+version = "1.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53"
+
+[[package]]
+name = "heck"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
+
+[[package]]
+name = "libc"
+version = "0.2.189"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2"
+
+[[package]]
+name = "matrixmultiply"
+version = "0.3.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3f607c237553f086e7043417a51df26b2eb899d3caff94e6a67592ff992fedc7"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
+[[package]]
+name = "ndarray"
+version = "0.17.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "520080814a7a6b4a6e9070823bb24b4531daac8c4627e08ba5de8c5ef2f2752d"
+dependencies = [
+ "matrixmultiply",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "portable-atomic",
+ "portable-atomic-util",
+ "rawpointer",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.4.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.47"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7ce2d95d4b3734dc35aa2f45e1aa22cd416814592a4f9d9205e11affd5b8e10b"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "numpy"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6a5b15d63a5ff39e378daed0e1340d3a5964703ea9712eb09a0dc66fade996f4"
+dependencies = [
+ "libc",
+ "ndarray",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "pyo3",
+ "pyo3-build-config",
+ "rustc-hash",
+]
+
+[[package]]
+name = "once_cell"
+version = "1.21.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50"
+
+[[package]]
+name = "portable-atomic"
+version = "1.15.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "05c8b63e8d9609db387f0324918f81d68fe27748f084ef092fb35954d0539a85"
+
+[[package]]
+name = "portable-atomic-util"
+version = "0.2.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "10ab3eb7f3becc3a1cbc4f2c6f20267996cfc1a6467a873763411b136a122715"
+dependencies = [
+ "portable-atomic",
+]
+
+[[package]]
+name = "proc-macro2"
+version = "1.0.107"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "pyo3"
+version = "0.29.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4688ddedf473e32662b9b067670129a8afb8c18e351482c70d62ba4a88171e8b"
+dependencies = [
+ "libc",
+ "once_cell",
+ "portable-atomic",
+ "pyo3-build-config",
+ "pyo3-ffi",
+ "pyo3-macros",
+]
+
+[[package]]
+name = "pyo3-build-config"
+version = "0.29.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f41027e41b4bd03f6e60f9f417fe24a6341a6bb744edd62b6f709f2a52ea30e9"
+dependencies = [
+ "target-lexicon",
+]
+
+[[package]]
+name = "pyo3-ffi"
+version = "0.29.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e591a95526fead067432c3b3a33fc74770b87b1e04e73671090d9c2055a2b327"
+dependencies = [
+ "libc",
+ "pyo3-build-config",
+]
+
+[[package]]
+name = "pyo3-macros"
+version = "0.29.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73225868fc1cd84eef2c3c230ddb91273bf1de46aeb8a4248da76d32a0924a1c"
+dependencies = [
+ "proc-macro2",
+ "pyo3-macros-backend",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "pyo3-macros-backend"
+version = "0.29.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "571575aa3749fa6216757dd47d2a3e7ef360f329a40f0666a9fbd14889024952"
+dependencies = [
+ "heck",
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "quote"
+version = "1.0.47"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001"
+dependencies = [
+ "proc-macro2",
+]
+
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
+[[package]]
+name = "rustc-hash"
+version = "2.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6b1e7f9a428571be2dc5bc0505c13fb6bf936822b894ec87abf8a08a4e51742d"
+
+[[package]]
+name = "simurg-core"
+version = "1.1.0"
+dependencies = [
+ "numpy",
+ "pyo3",
+]
+
+[[package]]
+name = "syn"
+version = "2.0.119"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "target-lexicon"
+version = "0.13.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "adb6935a6f5c20170eeceb1a3835a49e12e19d792f6dd344ccc76a985ca5a6ca"
+
+[[package]]
+name = "unicode-ident"
+version = "1.0.26"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d245f478577f809a851594d02313b640fb437e0bb33866753cff937863096954"
@@ -0,0 +1,20 @@
+[package]
+name = "simurg-core"
+version = "1.1.0"
+edition = "2021"
+license = "MIT"
+description = "Compiled kernels for SimurgArray (plan §18 P9)"
+publish = false
+
+[lib]
+name = "simurg_core"
+crate-type = ["cdylib"]
+
+[dependencies]
+numpy = "0.29.0"
+pyo3 = { version = "0.29", features = ["extension-module", "abi3-py312"] }
+
+[profile.release]
+opt-level = 3
+lto = true
+codegen-units = 1
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2026 Kemal Yaylali
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
@@ -0,0 +1,6 @@
+simurg-core: optional compiled kernels of SimurgArray (MIT, see LICENSE).
+
+The kernels port numerics that SimurgArray's pure-Python modules port from gtc2vcf (MIT, (C) Giulio Genovese) and
+Illumina BeadArrayFiles (BSD-2-Clause, (C) 2020 Illumina); full notices: THIRD_PARTY_NOTICES.md of SimurgArray.
+The binary statically links PyO3, rust-numpy, ndarray and their dependencies (MIT / Apache-2.0 / BSD-2-Clause);
+`cargo tree` lists them and each crate's licence ships with its source on crates.io.
@@ -0,0 +1,15 @@
+[build-system]
+requires = ["maturin>=1.7,<2"]
+build-backend = "maturin"
+
+[project]
+name = "simurg-core"
+description = "Compiled kernels for SimurgArray (optional; the pure-Python paths give identical results)"
+requires-python = ">=3.12"
+license = "MIT"
+license-files = ["LICENSE", "NOTICE"]
+dependencies = ["numpy>=1.26"]
+dynamic = ["version"]
+
+[tool.maturin]
+features = ["pyo3/extension-module"]
@@ -0,0 +1,275 @@
+//! BPM LocusEntry records — port of `simurg.formats.bpm._parse_locus_entries` (format: Illumina BeadArrayFiles
+//! `BeadPoolManifest.py`, BSD-2). Same field order, checks and error messages; strings of ≤ 64 bytes are shared
+//! Python objects, as in the Python parser.
+
+use std::collections::HashMap;
+
+use numpy::{IntoPyArray, PyArray1, PyArray2};
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+use pyo3::types::{PyList, PyString};
+
+pub const STR_COLS: usize = 16;
+// output column order (matches the Python `str_cols` list)
+const ILMN: usize = 0;
+const SNP: usize = 1;
+const CHR: usize = 2;
+const ISTR: usize = 3;
+const SSTR: usize = 4;
+const RSTR: usize = 5;
+const PLOIDY: usize = 6;
+const SPECIES: usize = 7;
+const BUILD: usize = 8;
+const SOURCE: usize = 9;
+const SVER: usize = 10;
+const CSTR: usize = 11;
+const SSEQ: usize = 12;
+const TOP: usize = 13;
+const PA: usize = 14;
+const PB: usize = 15;
+
+struct Cursor<'a> {
+    b: &'a [u8],
+    pos: usize,
+}
+
+fn err(msg: impl Into<String>) -> PyErr {
+    PyValueError::new_err(msg.into())
+}
+
+impl<'a> Cursor<'a> {
+    fn string(&mut self) -> PyResult<&'a [u8]> {
+        let size = self.b.len();
+        let mut p = self.pos;
+        if p >= size {
+            return Err(err("unexpected end of file in locus entry"));
+        }
+        let mut byte = self.b[p];
+        p += 1;
+        let mut ln = (byte & 0x7f) as usize;
+        let mut shift = 7;
+        while byte & 0x80 != 0 {
+            if p >= size || shift > 28 {
+                return Err(err("bad string length in locus entry"));
+            }
+            byte = self.b[p];
+            p += 1;
+            ln |= ((byte & 0x7f) as usize) << shift;
+            shift += 7;
+        }
+        let end = p + ln;
+        if end > size {
+            return Err(err("string runs past end of file in locus entry"));
+        }
+        self.pos = end;
+        Ok(&self.b[p..end])
+    }
+    fn i32(&mut self, what: &str) -> PyResult<i32> {
+        if self.pos + 4 > self.b.len() {
+            return Err(err(format!("unexpected end of file at {what}")));
+        }
+        let v = i32::from_le_bytes(self.b[self.pos..self.pos + 4].try_into().unwrap());
+        self.pos += 4;
+        Ok(v)
+    }
+}
+
+type Parsed<'py> = (
+    Vec<Bound<'py, PyList>>,
+    Bound<'py, PyArray1<i64>>,
+    Bound<'py, PyArray1<i32>>,
+    Bound<'py, PyArray1<i32>>,
+    Bound<'py, PyArray1<i16>>,
+    Bound<'py, PyArray1<i8>>,
+    Bound<'py, PyArray1<bool>>,
+    Bound<'py, PyArray1<i8>>,
+    Option<Bound<'py, PyArray2<f32>>>,
+    usize,
+);
+
+/// Parse `n` LocusEntry records starting at byte 0 of `buf`. Returns the 16 string columns, map_infos,
+/// address_a, address_b, assay types (−1 = infer from sequences, entry version < 6), expected clusters,
+/// intensity-only flags, entry versions, base fractions (None unless some entry is v7+) and the bytes consumed.
+/// Errors are ValueError with the Python parser's message.
+#[pyfunction]
+pub fn bpm_locus_entries<'py>(py: Python<'py>, buf: &[u8], n: usize, names: Vec<String>) -> PyResult<Parsed<'py>> {
+    let name_to_idx: HashMap<&str, usize> = names.iter().enumerate().map(|(i, s)| (s.as_str(), i)).collect();
+    let mut cols: Vec<Vec<Option<Py<PyString>>>> = (0..STR_COLS).map(|_| (0..n).map(|_| None).collect()).collect();
+    let mut intern: HashMap<&[u8], Py<PyString>> = HashMap::new();
+    let mut map_infos = vec![0i64; n];
+    let mut address_a = vec![0i32; n];
+    let mut address_b = vec![0i32; n];
+    let mut assay = vec![0i16; n];
+    let mut exp_clusters = vec![0i8; n];
+    let mut int_only = vec![false; n];
+    let mut versions = vec![0i8; n];
+    let mut fracs = vec![0f32; 4 * n];
+    let mut seen = vec![false; n];
+    let mut any_v7 = false;
+    let mut c = Cursor { b: buf, pos: 0 };
+
+    let to_py = |raw: &'_ [u8], intern: &mut HashMap<&[u8], Py<PyString>>| -> Py<PyString> {
+        if let Some(s) = intern.get(raw) {
+            return s.clone_ref(py);
+        }
+        let s = PyString::new(py, &String::from_utf8_lossy(raw)).unbind();
+        s
+    };
+
+    for _ in 0..n {
+        let ver = c.i32("locus entry")?;
+        if ![4, 6, 7, 8].contains(&ver) {
+            return Err(err(format!("unsupported locus entry version {ver}")));
+        }
+        let ilmn_id = c.string()?;
+        let name = c.string()?;
+        for _ in 0..3 {
+            c.string()?;
+        }
+        c.i32("locus index")?;
+        c.string()?;
+        let ilmn_strand = c.string()?;
+        let snp = c.string()?;
+        let chrom = c.string()?;
+        let ploidy = c.string()?;
+        let species = c.string()?;
+        let map_info = c.string()?;
+        let top_seq = c.string()?;
+        let cust_strand = c.string()?;
+        if c.pos + 8 > buf.len() {
+            return Err(err("unexpected end of file at addresses"));
+        }
+        let a_addr = c.i32("addresses")?;
+        let b_addr = c.i32("addresses")?;
+        let probe_a = c.string()?;
+        let probe_b = c.string()?;
+        let build = c.string()?;
+        let source = c.string()?;
+        let source_ver = c.string()?;
+        let source_strand = c.string()?;
+        let source_seq = c.string()?;
+        let (mut at, mut exp_cl, mut io) = (-1i16, 0i8, 0u8);
+        if ver >= 6 {
+            if c.pos + 4 > buf.len() {
+                return Err(err("unexpected end of file at assay type"));
+            }
+            exp_cl = buf[c.pos + 1] as i8;
+            io = buf[c.pos + 2];
+            at = buf[c.pos + 3] as i16;
+            c.pos += 4;
+            if at > 2 {
+                return Err(err(format!("invalid assay type {at}")));
+            }
+            if (b_addr == 0) != (at == 0) {
+                return Err(err("assay type is inconsistent with address B"));
+            }
+        }
+        let mut f4: Option<[f32; 4]> = None;
+        if ver >= 7 {
+            if c.pos + 16 > buf.len() {
+                return Err(err("unexpected end of file at base fractions"));
+            }
+            let mut f = [0f32; 4];
+            for (k, v) in f.iter_mut().enumerate() {
+                *v = f32::from_le_bytes(buf[c.pos + 4 * k..c.pos + 4 * k + 4].try_into().unwrap());
+            }
+            f4 = Some(f);
+            c.pos += 16;
+        }
+        let ref_strand: &[u8] = if ver >= 8 { c.string()? } else { b"" };
+
+        let name_s = String::from_utf8_lossy(name);
+        let idx = match name_to_idx.get(name_s.as_ref()) {
+            Some(&i) => i,
+            None => return Err(err(format!("locus entry {} not in the names table", py_repr(&name_s)))),
+        };
+        if seen[idx] {
+            return Err(err(format!("locus entry {} appears twice", py_repr(&name_s))));
+        }
+        seen[idx] = true;
+        for (col, raw) in [
+            (ILMN, ilmn_id),
+            (SNP, snp),
+            (CHR, chrom),
+            (ISTR, ilmn_strand),
+            (SSTR, source_strand),
+            (RSTR, ref_strand),
+            (PLOIDY, ploidy),
+            (SPECIES, species),
+            (BUILD, build),
+            (SOURCE, source),
+            (SVER, source_ver),
+            (CSTR, cust_strand),
+            (SSEQ, source_seq),
+            (TOP, top_seq),
+            (PA, probe_a),
+            (PB, probe_b),
+        ] {
+            let s = to_py(raw, &mut intern);
+            if raw.len() <= 64 && !intern.contains_key(raw) {
+                intern.insert(raw, s.clone_ref(py));
+            }
+            cols[col][idx] = Some(s);
+        }
+        map_infos[idx] = if map_info.is_empty() {
+            0
+        } else {
+            let t = String::from_utf8_lossy(map_info);
+            match t.trim().parse::<i64>() {
+                Ok(v) => v,
+                Err(_) => {
+                    return Err(err(format!(
+                        "non-integer MapInfo {} for {}",
+                        py_repr(&t),
+                        py_repr(&name_s)
+                    )))
+                }
+            }
+        };
+        address_a[idx] = a_addr;
+        address_b[idx] = b_addr;
+        assay[idx] = at;
+        exp_clusters[idx] = exp_cl;
+        int_only[idx] = io != 0;
+        versions[idx] = ver as i8;
+        if let Some(f) = f4 {
+            fracs[4 * idx..4 * idx + 4].copy_from_slice(&f);
+            any_v7 = true;
+        }
+    }
+    let missing = seen.iter().filter(|s| !**s).count();
+    if missing > 0 {
+        return Err(err(format!("{missing} loci in the names table have no locus entry")));
+    }
+    let empty = PyString::new(py, "").unbind();
+    let lists = cols
+        .into_iter()
+        .map(|col| PyList::new(py, col.into_iter().map(|s| s.unwrap_or_else(|| empty.clone_ref(py)))))
+        .collect::<PyResult<Vec<_>>>()?;
+    let fr = if any_v7 {
+        Some(numpy::ndarray::Array2::from_shape_vec((n, 4), fracs).map_err(|e| err(e.to_string()))?.into_pyarray(py))
+    } else {
+        None
+    };
+    Ok((
+        lists,
+        map_infos.into_pyarray(py),
+        address_a.into_pyarray(py),
+        address_b.into_pyarray(py),
+        assay.into_pyarray(py),
+        exp_clusters.into_pyarray(py),
+        int_only.into_pyarray(py),
+        versions.into_pyarray(py),
+        fr,
+        c.pos,
+    ))
+}
+
+/// Python `repr()` of a str for error messages (single quotes unless the text contains one).
+fn py_repr(s: &str) -> String {
+    if s.contains('\'') && !s.contains('"') {
+        format!("\"{s}\"")
+    } else {
+        format!("'{}'", s.replace('\\', "\\\\").replace('\'', "\\'"))
+    }
+}
@@ -0,0 +1,357 @@
+//! Fused per-locus ports of `simurg.genotype.gencall.call_genotypes` and `simurg.genotype.baf_lrr.baf_lrr`
+//! (which port gtc2vcf `gen_std_flair`/`modilik`/`compute_score_call_prelim`/`gencall_score_map`/`make_calls`
+//! and `compute_baf_lrr`, MIT, © Giulio Genovese). Every intermediate keeps the NumPy reference's type (f32 or
+//! f64) and evaluation order; NumPy's `fmin`/`fmax` are `f32::min`/`max`; libm calls use CPython `math`
+//! semantics. One pass per locus instead of dozens of full-array temporaries.
+
+use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+
+use crate::{math_exp, math_log, math_pow};
+
+const M_2_SQRTPI: f64 = 1.12837916709551257390;
+const M_SQRT1_2: f64 = 0.707106781186547524401;
+const NEG_NAN32: f32 = f32::from_bits(0xffc0_0000);
+
+#[inline]
+fn smf(x: f32, a: f32, b: f32) -> f32 {
+    let mid = (a + b) / 2.0;
+    if a >= b {
+        return if x >= mid { 1.0 } else { 0.0 };
+    }
+    let u = (x - a) / (b - a);
+    let v = (x - b) / (b - a);
+    if x <= a {
+        0.0
+    } else if a < x && x <= mid {
+        2.0 * (u * u)
+    } else if mid < x && x <= b {
+        1.0 - 2.0 * (v * v)
+    } else {
+        1.0
+    }
+}
+
+#[inline]
+fn zmf(x: f32, a: f32, b: f32) -> f32 {
+    let mid = (a + b) / 2.0;
+    let u = (x - a) / (b - a);
+    let v = (x - b) / (b - a);
+    if x <= a {
+        1.0
+    } else if a < x && x <= mid {
+        1.0 - 2.0 * (u * u)
+    } else if mid < x && x <= b {
+        2.0 * (v * v)
+    } else {
+        0.0
+    }
+}
+
+/// Normal pdf pieces of `_pdf`: (full, peak, bad) for value x, mean mu, sd sigma.
+#[inline]
+fn pdf(x: f32, mu: f32, sigma: f32) -> (f64, f64, bool) {
+    let z = ((x - mu) / sigma) as f64;
+    let e = math_exp(-0.5 * (z * z));
+    let sd = sigma as f64;
+    let full = e * 0.5 * M_2_SQRTPI * M_SQRT1_2 / sd;
+    let peak = (0.5 * M_2_SQRTPI * M_SQRT1_2) / sd;
+    (full, peak, sigma <= 0.0)
+}
+
+#[derive(Clone, Copy)]
+enum Side {
+    Left,
+    Right,
+    Full,
+}
+
+#[inline]
+fn get(p: (f64, f64, bool), x: f32, mu: f32, side: Side) -> f64 {
+    let (full, peak, bad) = p;
+    if bad {
+        return f64::NAN;
+    }
+    match side {
+        Side::Left => {
+            if x < mu {
+                peak
+            } else {
+                full
+            }
+        }
+        Side::Right => {
+            if x > mu {
+                peak
+            } else {
+                full
+            }
+        }
+        Side::Full => full,
+    }
+}
+
+fn check(len: usize, n: usize, what: &str) -> PyResult<()> {
+    if len != n {
+        return Err(PyValueError::new_err(format!("gencall: {what} has length {len}, expected {n}")));
+    }
+    Ok(())
+}
+
+/// Genotypes (0 NC, 1 AA, 2 AB, 3 BB) and GenCall scores; (n, 3) cluster arrays are passed flattened (C order).
+#[pyfunction]
+#[allow(clippy::too_many_arguments)]
+pub fn call_genotypes<'py>(
+    py: Python<'py>,
+    theta: PyReadonlyArray1<'py, f32>,
+    r: PyReadonlyArray1<'py, f32>,
+    raw_zero: PyReadonlyArray1<'py, bool>,
+    present: PyReadonlyArray1<'py, bool>,
+    n_cl: PyReadonlyArray1<'py, i32>,
+    theta_mean: PyReadonlyArray1<'py, f32>,
+    theta_dev: PyReadonlyArray1<'py, f32>,
+    r_mean: PyReadonlyArray1<'py, f32>,
+    r_dev: PyReadonlyArray1<'py, f32>,
+    intensity_threshold: PyReadonlyArray1<'py, f32>,
+    gentrain: PyReadonlyArray1<'py, f32>,
+    intensity_only: PyReadonlyArray1<'py, bool>,
+    cutoff: f32,
+) -> PyResult<(Bound<'py, PyArray1<u8>>, Bound<'py, PyArray1<f32>>)> {
+    let (th, rr, rz, pres) = (theta.as_slice()?, r.as_slice()?, raw_zero.as_slice()?, present.as_slice()?);
+    let (ncl, tmv, tdv, rmv, rdv) =
+        (n_cl.as_slice()?, theta_mean.as_slice()?, theta_dev.as_slice()?, r_mean.as_slice()?, r_dev.as_slice()?);
+    let (ithr, gt, io) = (intensity_threshold.as_slice()?, gentrain.as_slice()?, intensity_only.as_slice()?);
+    let n = th.len();
+    for (l, w) in [(rr.len(), "r"), (rz.len(), "raw_zero"), (pres.len(), "present"), (ithr.len(), "threshold"),
+                   (gt.len(), "gentrain"), (io.len(), "intensity_only")] {
+        check(l, n, w)?;
+    }
+    for (l, w) in [(ncl.len(), "n"), (tmv.len(), "theta_mean"), (tdv.len(), "theta_dev"), (rmv.len(), "r_mean"),
+                   (rdv.len(), "r_dev")] {
+        check(l, 3 * n, w)?;
+    }
+    let tiny_d = f64::MIN_POSITIVE * f64::EPSILON;
+    let third: f32 = 1.0 / 3.0;
+    let mut geno = vec![0u8; n];
+    let mut scores = vec![0f32; n];
+
+    for i in 0..n {
+        let t = th[i];
+        let r = rr[i];
+        let tm = [tmv[3 * i], tmv[3 * i + 1], tmv[3 * i + 2]];
+        let td = [tdv[3 * i], tdv[3 * i + 1], tdv[3 * i + 2]];
+        let rm = [rmv[3 * i], rmv[3 * i + 1], rmv[3 * i + 2]];
+        let rd = [rdv[3 * i], rdv[3 * i + 1], rdv[3 * i + 2]];
+        let nc = [ncl[3 * i], ncl[3 * i + 1], ncl[3 * i + 2]];
+
+        // gen_std_flair
+        let z1 = 0.5f32 * (tm[1] - tm[0]) / (td[1] + td[0]);
+        let z2 = 0.5f32 * (tm[2] - tm[1]) / (td[2] + td[1]);
+        let mz = z1.min(z2);
+        let alpha = (third * mz).max(1.0);
+        let beta = (third * mz).min(1.0);
+        let eta = alpha * beta;
+        let mut ftd = [0f32; 3];
+        for k in 0..3 {
+            let v = td[k] * eta;
+            ftd[k] = if v < 0.02 { 0.02 } else { v };
+        }
+        let med = td[0].min(td[1]).max(td[0].max(td[1]).min(td[2]));
+        let min_r = if 0.2f32 < med { med } else { 0.2 };
+        let mut frd = [0f32; 3];
+        for k in 0..3 {
+            let mrk0 = rm[k] / 7.0;
+            let mrk = if mrk0 < min_r { min_r } else { mrk0 };
+            frd[k] = if rd[k] < mrk { mrk } else { rd[k] };
+        }
+
+        // modilik
+        let d: [f64; 3] = std::array::from_fn(|k| 100.0 * ((t - tm[k]).abs() as f64) + ((r - rm[k]).abs() as f64));
+        let (daa, dab, dbb) = (d[0], d[1], d[2]);
+        let ok_t = !t.is_nan();
+        let pt: [(f64, f64, bool); 3] = std::array::from_fn(|k| pdf(t, tm[k], ftd[k]));
+        let pr: [(f64, f64, bool); 3] = std::array::from_fn(|k| pdf(r, rm[k], frd[k]));
+        let pdf2 = |k: usize, st: Side, sr: Side| get(pt[k], t, tm[k], st) * get(pr[k], r, rm[k], sr);
+        use Side::{Full, Left, Right};
+        let mut lik = [0f64; 3];
+        let mut covered = false;
+        if ok_t {
+            let cases: [(bool, [(usize, Side, Side); 2]); 6] = [
+                (daa <= dbb && dab <= dbb && rm[0] <= rm[1], [(0, Left, Left), (1, Right, Right)]),
+                (daa <= dab && dbb <= dab && rm[0] <= rm[2], [(0, Left, Left), (2, Right, Right)]),
+                (dab <= daa && dbb <= daa && rm[1] <= rm[2], [(1, Left, Left), (2, Right, Right)]),
+                (daa <= dbb && dab <= dbb && rm[0] > rm[1], [(0, Left, Right), (1, Right, Left)]),
+                (dab <= daa && dbb <= daa && rm[1] > rm[2], [(1, Left, Right), (2, Right, Left)]),
+                (daa <= dab && dbb <= dab && rm[0] > rm[2], [(0, Left, Right), (2, Right, Left)]),
+            ];
+            for (cond, assigns) in cases {
+                if cond {
+                    for (k, st, sr) in assigns {
+                        lik[k] = pdf2(k, st, sr);
+                    }
+                    covered = true;
+                }
+            }
+        }
+        if !covered {
+            for (k, l) in lik.iter_mut().enumerate() {
+                *l = pdf2(k, Full, Full);
+            }
+        }
+        for l in lik.iter_mut() {
+            if l.is_nan() {
+                *l = 0.0;
+            }
+        }
+        if !ok_t {
+            lik = [f64::NAN; 3];
+        }
+
+        // priors and posteriors
+        let ntot = nc[0] as i64 + nc[1] as i64 + nc[2] as i64;
+        let p: [f64; 3] = std::array::from_fn(|k| (((nc[k] as f32) / (ntot as f32)) as f64 + 1.0) / 4.0);
+        let lp: [f64; 3] = std::array::from_fn(|k| lik[k] * p[k]);
+        let ev = (lp[0] + lp[1]) + lp[2];
+        let (aa, ab, bb) = (lp[0] / ev, lp[1] / ev, lp[2] / ev);
+        let mut iap: u8 = if aa >= ab && aa >= bb {
+            1
+        } else if ab >= aa && ab >= bb {
+            2
+        } else if bb >= aa && bb >= ab {
+            3
+        } else {
+            0
+        };
+        let mut mx = if aa > ab { aa } else { ab };
+        let mut scnd = if aa > ab { ab } else { aa };
+        let b_gt_mx = bb > mx;
+        scnd = if b_gt_mx { mx } else if bb > scnd { bb } else { scnd };
+        mx = if b_gt_mx { bb } else { mx };
+        let ap_lod = math_log(mx / (tiny_d + scnd), true);
+        let score_ap = smf(ap_lod as f32, 0.0, 2.0);
+
+        // sub-scores
+        let score_r1 = smf(r, 0.0, 0.1);
+        let nclus = nc.iter().filter(|&&v| v > 0).count();
+        let misclust: f32 = if nclus == 1 && nc[1] == 0 {
+            0.7
+        } else if nclus != 3 {
+            0.95
+        } else {
+            1.0
+        };
+        let g = (iap as i64 - 1).clamp(0, 2) as usize;
+        let (gm_t, gd_t, gm_r, gd_r) = (tm[g], ftd[g], rm[g], frd[g]);
+        let st_aa = zmf(t, gm_t + 2.0 * gd_t, gm_t + 6.0 * gd_t);
+        let st_ab = zmf(((t - gm_t) / gd_t).abs(), 2.0, 6.0);
+        let st_bb = smf(t, gm_t - 6.0 * gd_t, gm_t - 2.0 * gd_t);
+        let mut score_t = match iap {
+            1 => st_aa,
+            2 => st_ab,
+            3 => st_bb,
+            _ => 1.0,
+        };
+        let r2 = smf(r, 0.0, gm_r / 10.0);
+        let r3 = smf(r, gm_r - 6.0 * gd_r, gm_r - 2.0 * gd_r);
+        let r4 = zmf(r, gm_r + 6.0 * gd_r, gm_r + 12.0 * gd_r);
+        let mut score_r = if iap > 0 { ((score_r1 * r2) * r3) * r4 } else { 0.0 };
+        let mut score_n: f32 = 1.0;
+        if gt[i].is_nan() {
+            score_t = f32::NAN;
+            score_r = f32::NAN;
+        }
+        if !ok_t {
+            score_t = f32::NAN;
+            score_r = f32::NAN;
+            score_n = f32::NAN;
+        }
+        let mut prelim = ((((score_ap * gt[i]) * score_t) * score_r) * score_n) * misclust;
+        if r < ithr[i] {
+            prelim = f32::NAN;
+            iap = 0;
+        }
+
+        // gencall_score_map
+        let x = prelim as f64;
+        let q = (x - 1.08) / 0.504;
+        let gbell = 1.0 / (1.0 + math_pow(q * q, 1.71));
+        let mut score = (math_pow(x, 0.35) * gbell) as f32;
+        if score.is_nan() {
+            score = 0.0;
+        }
+        if rz[i] || !pres[i] {
+            score = 0.0;
+            iap = 0;
+        }
+        geno[i] = if score < cutoff || io[i] { 0 } else { iap };
+        scores[i] = score;
+    }
+    Ok((geno.into_pyarray(py), scores.into_pyarray(py)))
+}
+
+/// BAF and LRR (float32, canonical NaN 0xFFC00000 where undefined).
+#[pyfunction]
+#[allow(clippy::too_many_arguments)]
+pub fn baf_lrr<'py>(
+    py: Python<'py>,
+    theta: PyReadonlyArray1<'py, f32>,
+    r: PyReadonlyArray1<'py, f32>,
+    raw_zero: PyReadonlyArray1<'py, bool>,
+    present: PyReadonlyArray1<'py, bool>,
+    theta_mean: PyReadonlyArray1<'py, f32>,
+    r_mean: PyReadonlyArray1<'py, f32>,
+    mean_r: PyReadonlyArray1<'py, f32>,
+    intensity_only: PyReadonlyArray1<'py, bool>,
+) -> PyResult<(Bound<'py, PyArray1<f32>>, Bound<'py, PyArray1<f32>>)> {
+    let (th, rr, rz, pres) = (theta.as_slice()?, r.as_slice()?, raw_zero.as_slice()?, present.as_slice()?);
+    let (tmv, rmv, mr, io) = (theta_mean.as_slice()?, r_mean.as_slice()?, mean_r.as_slice()?, intensity_only.as_slice()?);
+    let n = th.len();
+    for (l, w) in [(rr.len(), "r"), (rz.len(), "raw_zero"), (pres.len(), "present"), (mr.len(), "mean_r"), (io.len(), "intensity_only")] {
+        check(l, n, w)?;
+    }
+    check(tmv.len(), 3 * n, "theta_mean")?;
+    check(rmv.len(), 3 * n, "r_mean")?;
+    let mut baf = vec![0f32; n];
+    let mut lrr = vec![0f32; n];
+    for i in 0..n {
+        let (t, r) = (th[i], rr[i]);
+        let (at, bt, ct) = (tmv[3 * i], tmv[3 * i + 1], tmv[3 * i + 2]);
+        let (ar, br, cr) = (rmv[3 * i], rmv[3 * i + 1], rmv[3 * i + 2]);
+        let (eq, lo, hi) = (t == bt, t < bt, t > bt);
+        let (r_ref, mut b) = if eq {
+            (br, 0.5f32)
+        } else if lo {
+            (ar + ((t - at) * (ar - br)) / (at - bt), ((t - at) / (bt - at)) * 0.5)
+        } else if hi {
+            (br + ((t - bt) * (cr - br)) / (ct - bt), 0.5 + ((t - bt) / (ct - bt)) * 0.5)
+        } else {
+            (f32::NAN, NEG_NAN32)
+        };
+        let denom = if io[i] { mr[i] } else { r_ref };
+        let ratio = (r / denom) as f64;
+        let mut l = math_log(ratio, false) as f32;
+        if r == 0.0 {
+            l = -f32::MAX;
+        }
+        if !(eq || lo || hi) {
+            l = NEG_NAN32;
+        }
+        b = if b < 0.0 {
+            0.0
+        } else if b > 1.0 {
+            1.0
+        } else {
+            b
+        };
+        if rz[i] || !pres[i] || b.is_nan() {
+            b = NEG_NAN32;
+        }
+        if rz[i] || !pres[i] || l.is_nan() {
+            l = NEG_NAN32;
+        }
+        baf[i] = b;
+        lrr[i] = l;
+    }
+    Ok((baf.into_pyarray(py), lrr.into_pyarray(py)))
+}
@@ -0,0 +1,226 @@
+//! Compiled kernels for SimurgArray (plan §18 P9, §19).
+//!
+//! Every kernel is a drop-in replacement for a pure-Python/NumPy reference in the `simurg` package and returns
+//! bit-identical results (checked by `tests/unit/test_core_kernels.py`):
+//!
+//! * `exp`, `log10`, `log2`, `pow` — element-wise glibc libm calls with the error semantics of CPython's `math`
+//!   module as used by `simurg.genotype.gencall._c_call` (NaN input returned as is, a domain error gives the
+//!   canonical +NaN, `log10(0)` = −inf, overflow = +inf). The libm symbols are the ones CPython's `math` calls, so
+//!   results equal the C reference implementations (gtc2vcf, MIT) that the Python port follows.
+//! * `closest_binned` — the per-bin candidate scan of the vendor's binned nearest-neighbour search
+//!   (`simurg.signal.normalization.fit._closest_binned`): float32 squared distances, first minimum wins.
+//!
+//! * `call_genotypes`, `baf_lrr` — fused per-locus GenCall calling/scoring and BAF/LRR (`gencall.rs`).
+//! * `bpm_locus_entries` — the BPM LocusEntry parser (`simurg.formats.bpm._parse_locus_entries`).
+//!
+//! No floating-point contraction (FMA) happens: Rust never fuses `a*b+c` implicitly.
+
+mod bpm;
+mod gencall;
+
+use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+
+#[link(name = "m")]
+unsafe extern "C" {
+    #[link_name = "exp"]
+    fn c_exp(x: f64) -> f64;
+    #[link_name = "log10"]
+    fn c_log10(x: f64) -> f64;
+    #[link_name = "log2"]
+    fn c_log2(x: f64) -> f64;
+    #[link_name = "pow"]
+    fn c_pow(x: f64, y: f64) -> f64;
+}
+
+/// CPython's `float('nan')`.
+const PY_NAN: f64 = f64::from_bits(0x7ff8_0000_0000_0000);
+
+#[inline]
+pub(crate) fn math_exp(v: f64) -> f64 {
+    if v.is_nan() {
+        return v;
+    }
+    let r = unsafe { c_exp(v) };
+    if r.is_nan() { PY_NAN } else if r.is_infinite() && v.is_finite() { f64::INFINITY } else { r }
+}
+
+/// CPython `m_log10`/`m_log2` + `_c_call`: x > 0 → libm; 0 → −inf (log10) or NaN (log2); x < 0, −inf → NaN.
+#[inline]
+pub(crate) fn math_log(v: f64, ten: bool) -> f64 {
+    if v.is_nan() {
+        return v;
+    }
+    if v.is_infinite() {
+        return if v > 0.0 { v } else { PY_NAN };
+    }
+    if v > 0.0 {
+        return unsafe { if ten { c_log10(v) } else { c_log2(v) } };
+    }
+    if v == 0.0 && ten { f64::NEG_INFINITY } else { PY_NAN }
+}
+
+/// CPython `math.pow(v, e)` for a finite exponent e > 0 that is not an integer, + `_c_call`.
+#[inline]
+pub(crate) fn math_pow(v: f64, e: f64) -> f64 {
+    if v.is_nan() {
+        return v;
+    }
+    if v.is_infinite() {
+        return f64::INFINITY;
+    }
+    if v == 0.0 {
+        return 0.0;
+    }
+    let r = unsafe { c_pow(v, e) };
+    if r.is_nan() { PY_NAN } else if r.is_infinite() { f64::INFINITY } else { r }
+}
+
+fn map_f64<'py>(py: Python<'py>, x: PyReadonlyArray1<'py, f64>, f: impl Fn(f64) -> f64) -> PyResult<Bound<'py, PyArray1<f64>>> {
+    let s = x.as_slice()?;
+    let out: Vec<f64> = s.iter().map(|&v| f(v)).collect();
+    Ok(out.into_pyarray(py))
+}
+
+#[pyfunction]
+fn exp<'py>(py: Python<'py>, x: PyReadonlyArray1<'py, f64>) -> PyResult<Bound<'py, PyArray1<f64>>> {
+    map_f64(py, x, math_exp)
+}
+
+#[pyfunction]
+fn log10<'py>(py: Python<'py>, x: PyReadonlyArray1<'py, f64>) -> PyResult<Bound<'py, PyArray1<f64>>> {
+    map_f64(py, x, |v| math_log(v, true))
+}
+
+#[pyfunction]
+fn log2<'py>(py: Python<'py>, x: PyReadonlyArray1<'py, f64>) -> PyResult<Bound<'py, PyArray1<f64>>> {
+    map_f64(py, x, |v| math_log(v, false))
+}
+
+#[pyfunction]
+fn pow<'py>(py: Python<'py>, x: PyReadonlyArray1<'py, f64>, e: f64) -> PyResult<Bound<'py, PyArray1<f64>>> {
+    if !(e.is_finite() && e > 0.0 && e.fract() != 0.0) {
+        return Err(PyValueError::new_err("pow kernel: exponent must be finite, positive and not an integer"));
+    }
+    map_f64(py, x, |v| math_pow(v, e))
+}
+
+/// For each sweep point `p` searching bin `choice[p]` (0..11 regular, 12..22 shifted; any other value: none), the
+/// index of the closest candidate among `eligible[i]` with `b1[i] == choice` (regular) or `b2[i] == choice − 12`
+/// (shifted), in ascending order; squared distance in float32; −1 if no candidate or the best is not < 1e20.
+#[pyfunction]
+#[allow(clippy::too_many_arguments)]
+fn closest_binned<'py>(
+    py: Python<'py>,
+    raw_x: PyReadonlyArray1<'py, f32>,
+    raw_y: PyReadonlyArray1<'py, f32>,
+    axis_x: PyReadonlyArray1<'py, f32>,
+    axis_y: PyReadonlyArray1<'py, f32>,
+    choice: PyReadonlyArray1<'py, i64>,
+    eligible: PyReadonlyArray1<'py, i64>,
+    b1: PyReadonlyArray1<'py, i64>,
+    b2: PyReadonlyArray1<'py, i64>,
+) -> PyResult<Bound<'py, PyArray1<i64>>> {
+    let (rx, ry) = (raw_x.as_slice()?, raw_y.as_slice()?);
+    let (ax, ay, ch) = (axis_x.as_slice()?, axis_y.as_slice()?, choice.as_slice()?);
+    let (el, b1, b2) = (eligible.as_slice()?, b1.as_slice()?, b2.as_slice()?);
+    if ax.len() != ay.len() || ax.len() != ch.len() || el.len() != b1.len() || el.len() != b2.len() || rx.len() != ry.len() {
+        return Err(PyValueError::new_err("closest_binned: length mismatch"));
+    }
+    let mut groups: Vec<Vec<usize>> = vec![Vec::new(); 23];
+    for (i, &e) in el.iter().enumerate() {
+        let e = usize::try_from(e).map_err(|_| PyValueError::new_err("closest_binned: negative index"))?;
+        if e >= rx.len() {
+            return Err(PyValueError::new_err("closest_binned: index out of range"));
+        }
+        if (0..12).contains(&b1[i]) {
+            groups[b1[i] as usize].push(e);
+        }
+        if (0..11).contains(&b2[i]) {
+            groups[12 + b2[i] as usize].push(e);
+        }
+    }
+    let out: Vec<i64> = (0..ax.len())
+        .map(|p| {
+            let c = ch[p];
+            if !(0..23).contains(&c) {
+                return -1;
+            }
+            let cand = &groups[c as usize];
+            let mut best_j = usize::MAX;
+            let mut best = f32::INFINITY;
+            for (j, &e) in cand.iter().enumerate() {
+                let dx = rx[e] - ax[p];
+                let dy = ry[e] - ay[p];
+                let d = dx * dx + dy * dy;
+                // numpy argmin: the first NaN wins, otherwise the first strict minimum
+                if d.is_nan() {
+                    best_j = j;
+                    best = d;
+                    break;
+                }
+                if best_j == usize::MAX || d < best {
+                    best_j = j;
+                    best = d;
+                }
+            }
+            if best_j == usize::MAX || !((best as f64) < 1e20) { -1 } else { cand[best_j] as i64 }
+        })
+        .collect();
+    Ok(out.into_pyarray(py))
+}
+
+/// Exact nearest neighbour (`simurg.signal.normalization.fit._closest_slow`): for each sweep point (x, y) the
+/// index of the reference point with the smallest float32 squared distance `(x-xr)²+(y-yr)²`; first minimum wins,
+/// a NaN distance wins as in `numpy.argmin`. `xref` must not be empty.
+#[pyfunction]
+fn closest_exact<'py>(
+    py: Python<'py>,
+    xref: PyReadonlyArray1<'py, f32>,
+    yref: PyReadonlyArray1<'py, f32>,
+    x: PyReadonlyArray1<'py, f32>,
+    y: PyReadonlyArray1<'py, f32>,
+) -> PyResult<Bound<'py, PyArray1<i64>>> {
+    let (xr, yr, xs, ys) = (xref.as_slice()?, yref.as_slice()?, x.as_slice()?, y.as_slice()?);
+    if xr.is_empty() || xr.len() != yr.len() || xs.len() != ys.len() {
+        return Err(PyValueError::new_err("closest_exact: empty reference or length mismatch"));
+    }
+    let out: Vec<i64> = xs
+        .iter()
+        .zip(ys)
+        .map(|(&px, &py_)| {
+            let mut best_j = 0usize;
+            let mut best = f32::INFINITY;
+            for (j, (&rx, &ry)) in xr.iter().zip(yr).enumerate() {
+                let dx = px - rx;
+                let dy = py_ - ry;
+                let d = dx * dx + dy * dy;
+                if d.is_nan() {
+                    return j as i64;
+                }
+                if j == 0 || d < best {
+                    best_j = j;
+                    best = d;
+                }
+            }
+            best_j as i64
+        })
+        .collect();
+    Ok(out.into_pyarray(py))
+}
+
+#[pymodule]
+fn simurg_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add("__version__", env!("CARGO_PKG_VERSION"))?;
+    m.add_function(wrap_pyfunction!(exp, m)?)?;
+    m.add_function(wrap_pyfunction!(log10, m)?)?;
+    m.add_function(wrap_pyfunction!(log2, m)?)?;
+    m.add_function(wrap_pyfunction!(pow, m)?)?;
+    m.add_function(wrap_pyfunction!(closest_binned, m)?)?;
+    m.add_function(wrap_pyfunction!(closest_exact, m)?)?;
+    m.add_function(wrap_pyfunction!(bpm::bpm_locus_entries, m)?)?;
+    m.add_function(wrap_pyfunction!(gencall::call_genotypes, m)?)?;
+    m.add_function(wrap_pyfunction!(gencall::baf_lrr, m)?)?;
+    Ok(())
+}
@@ -0,0 +1,22 @@
+#!/usr/bin/env bash
+# Build the release artifacts into dist/:
+#   simurg-<v>-py3-none-any.whl + simurg-<v>.tar.gz           (pure Python)
+#   simurg_core-<v>-cp312-abi3-manylinux_2_28_{aarch64,x86_64}.whl   (compiled kernels, Python >= 3.12)
+# Portable Linux wheels are linked against glibc 2.28 with `maturin --zig` (cross-compiles the other architecture).
+# Needs: uv, a Rust toolchain (rustup), ~60 MB for zig on first use.
+#   bash scripts/build_dist.sh [aarch64|x86_64|all]      (default: all)
+set -euo pipefail
+cd "$(dirname "$0")/.."
+export PATH="$HOME/.cargo/bin:$PATH"
+WHICH=${1:-all}
+mkdir -p dist
+uv build -o dist                                   # sdist + pure wheel of simurg
+MATURIN="uvx --from maturin[zig]>=1.7,<2 maturin"
+for arch in aarch64 x86_64; do
+  [[ $WHICH == all || $WHICH == "$arch" ]] || continue
+  rustup target add "$arch-unknown-linux-gnu" > /dev/null
+  $MATURIN build --release --zig --compatibility manylinux_2_28 --target "$arch-unknown-linux-gnu" \
+    -m rust/simurg-core/Cargo.toml -o dist
+done
+( cd dist && sha256sum ./*.whl ./*.tar.gz > SHA256SUMS )
+ls -la dist
@@ -0,0 +1,125 @@
+#!/bin/bash
+# Phase P10 acceptance on a fresh Ubuntu 24.04 VM (plan §18: "clean install + tutorial works on fresh VM"):
+#   1. build the release artifacts (scripts/build_dist.sh: sdist, pure wheel, manylinux_2_28 kernel wheels for both
+#      architectures);
+#   2. install them with plain pip into new venvs of the system Python (with and without the kernels) and run the
+#      unit tests against the installed package;
+#   3. fetch the public demo data and run the tutorial (examples/gsa_demo.sh, 6 samples) with each install;
+#   4. build the Docker image and the conda packages, run genotype + gtc-to-vcf with them;
+#   5. content digests of every output (scripts/p10/digest.py) must agree across all installs.
+# Args: <unused data dir>  <results dir>. Run with OMNI=0.
+set -uo pipefail
+OUT=$2
+R=$OUT/p10; W=/work/p10; mkdir -p $R $W
+N=$(nproc)
+ARCH=$(uname -m)
+declare -A STEP
+step() { STEP[$1]=$2; echo "--- $1: $2"; }
+
+export DEBIAN_FRONTEND=noninteractive
+apt-get install -y -q gcc docker.io bsdextrautils > /dev/null
+curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1
+export PATH="$HOME/.cargo/bin:$PATH"
+
+echo "=== 1. build artifacts $(date -u +%T) ==="
+if bash scripts/build_dist.sh all > $R/build_dist.log 2>&1; then step build ok; else step build FAILED; tail -30 $R/build_dist.log; fi
+cp dist/*.whl dist/*.tar.gz dist/SHA256SUMS $R/ 2>/dev/null
+ls -la dist
+
+echo "=== 2. pip install into fresh venvs $(date -u +%T) ==="
+PY=/usr/bin/python3
+$PY --version
+CORE=$(ls dist/simurg_core-*manylinux_2_28_$ARCH.whl 2>/dev/null | head -1)
+$PY -m venv $W/fresh && $W/fresh/bin/pip install -q dist/simurg-*-py3-none-any.whl "$CORE" \
+  && step pip_install ok || step pip_install FAILED
+$PY -m venv $W/pure && $W/pure/bin/pip install -q dist/simurg-*-py3-none-any.whl \
+  && step pip_install_pure ok || step pip_install_pure FAILED
+$W/fresh/bin/python -c "from simurg._core import backend; print('fresh:', backend())"
+$W/pure/bin/python -c "from simurg._core import backend; print('pure:', backend())"
+$W/fresh/bin/simurg --version
+$W/fresh/bin/pip install -q pytest hypothesis
+# tests import the installed package (the checkout keeps its sources under src/, not on sys.path)
+( cd $W && cp -r /work/repo/tests /work/repo/docs . && $W/fresh/bin/python -m pytest tests/unit -q -p no:cacheprovider \
+    -k "not selftest_in_browser" > $R/unit_tests_installed.txt 2>&1 ) && step unit_tests_installed ok \
+  || { step unit_tests_installed FAILED; tail -20 $R/unit_tests_installed.txt; }
+tail -3 $R/unit_tests_installed.txt
+
+echo "=== 3. tutorial $(date -u +%T) ==="
+bash scripts/fetch_data.sh p0 > $W/fetch_p0.log 2>&1 && bash scripts/fetch_data.sh p6 > $W/fetch_p6.log 2>&1 \
+  && step fetch_data ok || { step fetch_data FAILED; tail -5 $W/fetch_p0.log $W/fetch_p6.log; }
+for v in fresh pure; do
+  t0=$SECONDS
+  SIMURG=$W/$v/bin/simurg THREADS=$N bash examples/gsa_demo.sh $W/tut_$v 6 > $R/tutorial_$v.log 2>&1 \
+    && step tutorial_$v "ok ($((SECONDS - t0)) s)" || { step tutorial_$v FAILED; tail -20 $R/tutorial_$v.log; }
+  $W/fresh/bin/python scripts/p10/digest.py $W/tut_$v/gtc --out $R/digest_gtc_$v.json
+  $W/fresh/bin/python scripts/p10/digest.py $W/tut_$v/vcf --out $R/digest_vcf_$v.json
+done
+
+echo "=== 4a. docker $(date -u +%T) ==="
+if [[ "${SKIP_DOCKER:-0}" != "1" ]]; then
+D=data/p0/gsa_v4_pgx
+BPM=$D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm
+CSV=$D/GSA-PGx-48v4-0_20079540_E2_csv/GSA-PGx-48v4-0_20079540_E2.csv
+EGT=$D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt
+FASTA=data/p6/GRCh38.primary.fa
+if docker build -q -t simurg:test . > $R/docker_build.log 2>&1; then
+  step docker_build ok
+  docker image ls simurg:test --format '{{.Size}}' | sed 's/^/image size: /'
+  DR="docker run --rm --user 0:0 -v /work:/work -w /work/repo simurg:test"
+  $DR --version
+  $DR genotype --bpm $BPM --egt $EGT --csv $CSV --idats $W/tut_fresh/idats --out $W/docker/gtc \
+      --qc-profile genotyping --threads $N --quiet --force \
+    && $DR gtc-to-vcf --bpm $BPM --csv $CSV --fasta $FASTA --gtcs $W/docker/gtc --out $W/docker/vcf --threads $N \
+      --quiet --force && step docker_run ok || step docker_run FAILED
+  $W/fresh/bin/python scripts/p10/digest.py $W/docker/gtc --out $R/digest_gtc_docker.json
+  $W/fresh/bin/python scripts/p10/digest.py $W/docker/vcf --out $R/digest_vcf_docker.json
+else
+  step docker_build FAILED; tail -30 $R/docker_build.log
+fi
+fi
+
+echo "=== 4b. conda $(date -u +%T) ==="
+if [[ "${SKIP_CONDA:-0}" != "1" ]]; then
+curl -sSfL -o $W/miniforge.sh "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-$ARCH.sh" \
+  && bash $W/miniforge.sh -b -p $W/conda > /dev/null
+C=$W/conda/bin
+if $C/conda install -y -q -n base conda-build > $R/conda.log 2>&1 \
+   && $C/conda build -c conda-forge --output-folder $W/condabld packaging/conda/simurg-core packaging/conda/simurg \
+      >> $R/conda.log 2>&1; then
+  step conda_build ok
+  $C/conda create -y -q -p $W/cenv -c $W/condabld -c conda-forge simurg simurg-core >> $R/conda.log 2>&1 \
+    && $W/cenv/bin/simurg --version \
+    && $W/cenv/bin/python -c "from simurg._core import backend; print('conda:', backend())" \
+    && $W/cenv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/tut_fresh/idats --out $W/conda_run/gtc \
+         --threads $N --quiet --force \
+    && step conda_run ok || { step conda_run FAILED; tail -20 $R/conda.log; }
+  $W/fresh/bin/python scripts/p10/digest.py $W/conda_run/gtc --out $R/digest_gtc_conda.json
+else
+  step conda_build FAILED; tail -40 $R/conda.log
+fi
+fi
+
+echo "=== 5. digests $(date -u +%T) ==="
+$W/fresh/bin/python - "$R" <<'PY' | tee $R/agreement.txt
+import json, sys
+from pathlib import Path
+r = Path(sys.argv[1])
+for kind in ("gtc", "vcf"):
+    runs = {p.stem.split("_", 2)[2]: json.loads(p.read_text())[kind] for p in sorted(r.glob(f"digest_{kind}_*.json"))}
+    if not runs:
+        continue
+    ref = runs.get("fresh", {})
+    if not ref:
+        print(f"{kind}: no reference (the pip-installed tutorial produced no outputs)")
+        continue
+    for name, d in runs.items():
+        if name != "fresh":
+            print(f"{kind}: {name} vs fresh: {'IDENTICAL' if d == ref else 'DIFFERENT'} ({len(d)} vs {len(ref)} files)")
+PY
+{
+  echo "{"; echo "  \"arch\": \"$ARCH\","
+  for k in "${!STEP[@]}"; do echo "  \"$k\": \"${STEP[$k]}\","; done
+  echo "  \"python\": \"$($PY --version)\""; echo "}"
+} > $R/steps.json
+cat $R/steps.json
+echo "=== task done $(date -u +%T) ==="
@@ -0,0 +1,3 @@
+#!/bin/bash
+# P10 acceptance without the Docker and conda parts (already verified): build, pip installs, tutorial, digests.
+SKIP_DOCKER=1 SKIP_CONDA=1 bash "$(dirname "$0")/p10.sh" "$@"
@@ -0,0 +1,54 @@
+#!/bin/bash
+# Phase P9 performance gate on a 32-vCPU Graviton instance, GSE165845 GSA arrays (360 samples):
+#   1. build the compiled kernels (rust/simurg-core) into the venv;
+#   2. thread scaling 1 → 32 with the kernels (bench/scaling.py, per-sample phase from the run log);
+#   3. pure-Python reference at 1 and 32 threads on a subset, and byte-identity of all 360 GTCs (kernels vs pure).
+# Args: <unused data dir>  <results dir>. Run with OMNI=0.
+set -uo pipefail
+OUT=$2
+W=/work/p9; R=$OUT/p9; mkdir -p $W $R
+N=$(nproc)
+
+echo "=== build simurg-core $(date -u +%T) ==="
+DEBIAN_FRONTEND=noninteractive apt-get install -y -q gcc > /dev/null  # Rust links through the C toolchain
+curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1
+export PATH="$HOME/.cargo/bin:$PATH"
+rustc --version
+VIRTUAL_ENV=$PWD/.venv uv pip install -q ./rust/simurg-core || { echo "!! simurg-core build failed"; exit 1; }
+.venv/bin/python -c "from simurg._core import backend; print('backend', backend())"
+.venv/bin/python -m pytest tests/unit/test_core_kernels.py -q -p no:cacheprovider 2>&1 | tail -2
+
+echo "=== data $(date -u +%T) ==="
+GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
+mkdir -p $W/gsa && cd $W/gsa
+for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
+         infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
+  curl -sSfL --retry 5 -O "$GSA/$z" && unzip -oq $z
+done
+cd - > /dev/null
+BPM=$W/gsa/GSA-24v1-0_C2.bpm; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
+.venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo --out $W/idats --workers 24 \
+  || exit 1
+
+echo "=== scaling with kernels $(date -u +%T) ==="
+.venv/bin/python bench/scaling.py --bpm $BPM --egt $EGT --idats $W/idats --out $W/rust \
+  --threads 1,2,4,8,16,$N 2>&1 | grep -v WARNING
+cp $W/rust/scaling.json $R/scaling_rust.json
+
+echo "=== pure Python reference $(date -u +%T) ==="
+.venv/bin/python bench/scaling.py --bpm $BPM --egt $EGT --idats $W/idats --out $W/python --threads 1 --samples 64 \
+  --pure-python 2>&1 | grep -v WARNING
+cp $W/python/scaling.json $R/scaling_python_t1.json
+.venv/bin/python bench/scaling.py --bpm $BPM --egt $EGT --idats $W/idats --out $W/python_all --threads $N \
+  --pure-python 2>&1 | grep -v WARNING
+cp $W/python_all/scaling.json $R/scaling_python_t$N.json
+
+echo "=== identity: 360 GTCs, kernels vs pure Python $(date -u +%T) ==="
+n=0; d=0
+for f in $(cd $W/rust/t$N && find samples -name '*.gtc'); do
+  n=$((n+1)); cmp -s $W/rust/t$N/$f $W/python_all/t$N/$f || { d=$((d+1)); echo "differs: $f"; }
+done
+s=$(cmp -s $W/rust/t$N/sample_summary.tsv $W/python_all/t$N/sample_summary.tsv && echo identical || echo DIFFERS)
+echo "{\"gtcs\": $n, \"differ\": $d, \"sample_summary\": \"$s\"}" | tee $R/identity.json
+lscpu | grep -i "model name\|^CPU(s)" > $R/cpu.txt; free -g >> $R/cpu.txt
+echo "=== task done $(date -u +%T) ==="
@@ -0,0 +1,29 @@
+#!/bin/bash
+# Phase P9 follow-up: thread scaling with the GTC outputs on a RAM disk (/dev/shm), to separate compute scaling
+# from the gp3 disk's write throughput (125 MB/s baseline; one GSA GTC is ~13 MB). Run with OMNI=0.
+set -uo pipefail
+OUT=$2
+W=/work/p9; R=$OUT/p9b; mkdir -p $W $R
+N=$(nproc)
+DEBIAN_FRONTEND=noninteractive apt-get install -y -q gcc > /dev/null
+curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1
+export PATH="$HOME/.cargo/bin:$PATH"
+VIRTUAL_ENV=$PWD/.venv uv pip install -q ./rust/simurg-core || { echo "!! simurg-core build failed"; exit 1; }
+.venv/bin/python -c "from simurg._core import backend; print('backend', backend())"
+GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
+mkdir -p $W/gsa && cd $W/gsa
+for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
+         infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
+  curl -sSfL --retry 5 -O "$GSA/$z" && unzip -oq $z
+done
+cd - > /dev/null
+BPM=$W/gsa/GSA-24v1-0_C2.bpm; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
+.venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl /dev/shm/geo --out /dev/shm/idats \
+  --workers 24 2>&1 | tail -2 || exit 1
+df -h /dev/shm | tail -1
+echo "=== scaling, outputs on /dev/shm $(date -u +%T) ==="
+.venv/bin/python bench/scaling.py --bpm $BPM --egt $EGT --idats /dev/shm/idats --out $W/shm --work /dev/shm/out \
+  --clean --threads 1,4,8,16,24,$N 2>&1 | grep -v WARNING
+cp $W/shm/scaling.json $R/scaling_rust_shm.json
+lscpu | grep -i "model name\|^CPU(s)" > $R/cpu.txt
+echo "=== task done $(date -u +%T) ==="
@@ -14,6 +14,8 @@ set -euo pipefail
 cd "$(dirname "$0")/../.."
 REGION=${REGION:-us-east-1}
 TYPE=${TYPE:-m7g.4xlarge}
+# Graviton families carry a "g" after the generation number (m7g, c7gn, r8g): arm64 image; everything else amd64
+if [[ $TYPE =~ ^[a-z]+[0-9]+g ]]; then ARCH=${ARCH:-arm64}; else ARCH=${ARCH:-amd64}; fi
 MAX_PRICE=${MAX_PRICE:-0.30}        # USD/hour spot cap (on-demand is 0.6528)
 MAX_MINUTES=${MAX_MINUTES:-360}     # hard cap: 6 hours
 DISK_GB=${DISK_GB:-250}
@@ -59,7 +61,7 @@ EOF
 launch)
  running=$($AWS ec2 describe-instances --filters Name=tag:Project,Values=$TAG \
    Name=instance-state-name,Values=pending,running --query 'length(Reservations[].Instances[])')
-  [[ "$running" == "0" ]] || { echo "a trial instance is already running; refusing to start another"; exit 1; }
+  [[ "$running" == "0" || "${PARALLEL:-0}" == "1" ]] || { echo "a trial instance is already running; refusing to start another (PARALLEL=1 overrides)"; exit 1; }
  RUN=run-$(date -u +%Y%m%dT%H%M%SZ)
  echo "-> packing repository at $(git rev-parse --short HEAD) -> s3://$BUCKET/$RUN/"
  git archive --format=tar.gz -o /tmp/simurg-repo.tar.gz HEAD
@@ -69,7 +71,7 @@ launch)
  sed -e "s/__BUCKET__/$BUCKET/" -e "s/__RUN__/$RUN/" -e "s/__MAX_MINUTES__/$MAX_MINUTES/" -e "s/__TASK__/$TASK/g" \
    -e "s/__OMNI__/$OMNI/" \
    scripts/cloud/userdata.sh > /tmp/simurg-userdata.sh
-  AMI=$($AWS ssm get-parameter --name /aws/service/canonical/ubuntu/server/24.04/stable/current/arm64/hvm/ebs-gp3/ami-id \
+  AMI=$($AWS ssm get-parameter --name /aws/service/canonical/ubuntu/server/24.04/stable/current/$ARCH/hvm/ebs-gp3/ami-id \
    --query Parameter.Value --output text)
  echo "-> task: $TASK"
  echo "-> launching $TYPE spot (cap \$$MAX_PRICE/h, hard limit ${MAX_MINUTES} min, ${DISK_GB} GB disk deleted on exit)"
@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+"""Content digests of SimurgArray outputs, for comparing runs across machines, architectures and installs (P10).
+
+GTC: every numeric field and sample statistic, NaNs canonicalised (their sign/payload is platform-specific, K-014);
+free-text metadata (dates, software version) is left out. VCF(.gz): all data lines, header lines excluded.
+TSV: the listed columns of sample_summary.tsv.
+
+    python scripts/p10/digest.py RUN_DIR [--out digest.json]
+"""
+
+from __future__ import annotations
+
+import argparse
+import gzip
+import hashlib
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+
+from simurg.formats.gtc import read_gtc
+
+GTC_FIELDS = ("normalization_transforms", "controls_x", "controls_y", "raw_x", "raw_y", "genotypes", "base_calls",
+              "genotype_scores", "b_allele_freqs", "logr_ratios", "percentiles_x", "percentiles_y", "sample_data",
+              "call_rate", "gender", "logr_dev", "gc10")  # fmt: skip
+SUMMARY_COLUMNS = ("sample_id", "status", "autosomal_call_rate", "call_rate", "logr_dev", "sex_est", "gc10",
+                   "gc50", "n_calls", "n_nocalls", "het_rate", "contamination", "functional_status", "control_status",
+                   "flags")  # fmt: skip
+
+
+def _canon(a: np.ndarray) -> bytes:
+    if a.dtype.names:  # structured: canonicalise float members
+        a = a.copy()
+        for n in a.dtype.names:
+            if np.issubdtype(a.dtype[n], np.floating):
+                v = a[n]
+                v[np.isnan(v)] = np.nan
+        return a.tobytes()
+    if np.issubdtype(a.dtype, np.floating):
+        a = np.where(np.isnan(a), np.array(np.nan, a.dtype), a)
+    return np.ascontiguousarray(a).tobytes()
+
+
+def gtc_digest(path: Path) -> dict[str, str]:
+    g = read_gtc(path)
+    out = {}
+    for f in GTC_FIELDS:
+        v = getattr(g, f)
+        b = _canon(np.asarray(v)) if not isinstance(v, bytes) else v
+        out[f] = hashlib.sha256(b).hexdigest()[:16]
+    return out
+
+
+def vcf_digest(path: Path) -> str:
+    h = hashlib.sha256()
+    with gzip.open(path, "rt") if path.suffix == ".gz" else open(path) as fh:
+        for line in fh:
+            if not line.startswith("#"):
+                h.update(line.encode())
+    return h.hexdigest()[:16]
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("run", type=Path)
+    ap.add_argument("--out", type=Path)
+    a = ap.parse_args()
+    d: dict[str, object] = {}
+    d["gtc"] = {p.stem: gtc_digest(p) for p in sorted(a.run.rglob("*.gtc"))}
+    d["vcf"] = {p.name: vcf_digest(p) for p in sorted(a.run.rglob("*.vcf.gz"))}
+    summ = {}
+    for p in sorted(a.run.rglob("sample_summary.tsv")):
+        rows = [ln.rstrip("\n").split("\t") for ln in p.open()]
+        idx = [rows[0].index(c) for c in SUMMARY_COLUMNS if c in rows[0]]
+        summ[str(p.relative_to(a.run))] = hashlib.sha256(
+            "\n".join("\t".join(r[i] for i in idx) for r in rows).encode()
+        ).hexdigest()[:16]
+    d["summary"] = summ
+    text = json.dumps(d, indent=1, sort_keys=True)
+    (a.out.write_text(text) if a.out else sys.stdout.write(text + "\n"))
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
@@ -4,5 +4,5 @@ Specification: SIMURGARRAY_MASTER_PLAN.md (algorithms, roadmap) and
 SIMURGARRAY_FUNCTIONAL_SPEC.md (externally visible behaviour, item IDs such as CLI-021).
 """

-__version__ = "1.0.0"
+__version__ = "1.1.0"
 PRODUCT_NAME = "SimurgArray"
@@ -0,0 +1,24 @@
+"""Optional compiled kernels (``rust/simurg-core``, plan §18 P9).
+
+``core`` is the ``simurg_core`` extension module when it is installed, else ``None``; every caller keeps a pure
+Python/NumPy path with bit-identical results. ``SIMURG_PURE_PYTHON=1`` forces the pure paths (tests, debugging).
+"""
+
+from __future__ import annotations
+
+import os
+from types import ModuleType
+
+core: ModuleType | None
+if os.environ.get("SIMURG_PURE_PYTHON", "") not in ("", "0"):
+    core = None
+else:
+    try:
+        import simurg_core as core  # type: ignore[no-redef]
+    except ImportError:  # pragma: no cover - depends on the installation
+        core = None
+
+
+def backend() -> str:
+    """``"rust <version>"`` or ``"python"`` (recorded in run manifests)."""
+    return f"rust {core.__version__}" if core is not None else "python"
@@ -7,6 +7,7 @@ Later phases add the other commands of CLI-001.
 from __future__ import annotations

 import json
+import os
 import sys
 from collections.abc import Callable
 from pathlib import Path
@@ -58,6 +59,11 @@ def run_guarded(reporter: Reporter, fn: Callable[[], int]) -> None:
    code = ExitCode.OK
    try:
        code = ExitCode(fn())
+        sys.stdout.flush()  # surface a closed pipe here, not at interpreter exit
+    except BrokenPipeError:
+        # the consumer of stdout stopped reading (`| head`): silence further writes, exit like a SIGPIPE'd tool
+        os.dup2(os.open(os.devnull, os.O_WRONLY), sys.stdout.fileno())
+        code = ExitCode.SIGPIPE
    except SimurgError as e:
        reporter.error(e.code, e.message + (f"\n  hint: {e.hint}" if e.hint else ""), exc=e, file=e.path)
        code = e.exit_code
@@ -23,6 +23,7 @@ class ExitCode(IntEnum):
    PARTIAL = 3  # >= 1 sample failed and >= 1 succeeded
    ALL_FAILED = 4  # run completed but no sample succeeded
    SIGINT = 130
+    SIGPIPE = 141  # the reader of stdout went away (e.g. `simurg inspect … | head`); nothing is reported
    SIGTERM = 143


@@ -32,6 +32,7 @@ from typing import Any
 import numpy as np
 import numpy.typing as npt

+from simurg._core import core
 from simurg.core.errors import FormatError
 from simurg.formats.binutil import BinaryReader, load_bytes

@@ -223,7 +224,9 @@ def read_bpm(path: str | Path) -> Manifest:


 def _parse_locus_entries(r: BinaryReader, n: int, name_to_idx: dict[str, int]) -> dict[str, Any]:
-    """Hot loop: parse ``n`` LocusEntry records (optimised pure Python; Rust in phase P9)."""
+    """Hot loop: parse ``n`` LocusEntry records (compiled in ``simurg_core`` when installed; same results)."""
+    if core is not None:
+        return _parse_locus_entries_core(r, n, name_to_idx)
    view = memoryview(r._buf)
    size = r.size
    pos = r.pos
@@ -416,6 +419,53 @@ def _parse_locus_entries(r: BinaryReader, n: int, name_to_idx: dict[str, int]) -
    }


+_STR_COLS = (
+    "ilmn_ids",
+    "snps",
+    "chroms",
+    "ilmn_strands",
+    "source_strands",
+    "ref_strands",
+    "ploidies",
+    "species",
+    "genome_builds",
+    "sources",
+    "source_versions",
+    "customer_strands",
+    "source_seqs",
+    "top_genomic_seqs",
+    "allele_a_probe_seqs",
+    "allele_b_probe_seqs",
+)
+
+
+def _parse_locus_entries_core(r: BinaryReader, n: int, name_to_idx: dict[str, int]) -> dict[str, Any]:
+    """``_parse_locus_entries`` through ``simurg_core.bpm_locus_entries`` (identical columns and error messages)."""
+    assert core is not None
+    pos = r.pos
+    try:
+        lists, map_infos, a, b, assay, exp_cl, int_only, versions, fracs, used = core.bpm_locus_entries(
+            bytes(r._view[pos:]), n, list(name_to_idx)
+        )
+    except ValueError as e:
+        raise r.fail(str(e)) from None
+    out: dict[str, Any] = dict(zip(_STR_COLS, lists, strict=True))
+    for i in np.flatnonzero(assay < 0):
+        assay[i] = infer_assay_type(out["allele_a_probe_seqs"][i], out["allele_b_probe_seqs"][i], out["source_seqs"][i])
+    r.pos = pos + used
+    return {
+        **out,
+        "map_infos": map_infos,
+        "address_a": a,
+        "address_b": b,
+        "assay_types": assay.astype(np.uint8),
+        "exp_clusters": exp_cl,
+        "intensity_only": int_only,
+        "entry_versions": versions,
+        "base_fractions": fracs,
+    }
+
+
 _IUPAC = {
    "A": 1,
    "C": 2,
@@ -11,6 +11,7 @@ import math
 import numpy as np
 import numpy.typing as npt

+from simurg._core import core
 from simurg.genotype.clusters import ClusterArrays
 from simurg.genotype.gencall import _libm

@@ -22,6 +23,17 @@ NEG_NAN = np.frombuffer(np.uint32(0xFFC00000).tobytes(), dtype=np.float32)[0]  #
 def baf_lrr(
    theta: F32Arr, r: F32Arr, raw_zero: npt.NDArray[np.bool_], c: ClusterArrays, intensity_only: npt.NDArray[np.bool_]
 ) -> tuple[F32Arr, F32Arr]:
+    if core is not None:
+        return core.baf_lrr(  # type: ignore[no-any-return]
+            np.ascontiguousarray(theta, np.float32),
+            np.ascontiguousarray(r, np.float32),
+            np.ascontiguousarray(raw_zero, np.bool_),
+            np.ascontiguousarray(c.present, np.bool_),
+            np.ascontiguousarray(c.theta_mean, np.float32).ravel(),
+            np.ascontiguousarray(c.r_mean, np.float32).ravel(),
+            np.ascontiguousarray(c.mean_r, np.float32),
+            np.ascontiguousarray(intensity_only, np.bool_),
+        )
    t = np.asarray(theta, np.float32)
    rr = np.asarray(r, np.float32)
    at, bt, ct = c.theta_mean[:, 0], c.theta_mean[:, 1], c.theta_mean[:, 2]
@@ -15,6 +15,7 @@ from dataclasses import dataclass
 import numpy as np
 import numpy.typing as npt

+from simurg._core import core
 from simurg.genotype.clusters import ClusterArrays

 F32 = np.float32
@@ -42,8 +43,17 @@ def _c_call(fn, v: float) -> float:  # type: ignore[no-untyped-def]
        return math.inf


+_CORE_FN = {math.exp: "exp", math.log10: "log10", math.log2: "log2"}
+
+
 def _libm(fn, x: F64Arr) -> F64Arr:  # type: ignore[no-untyped-def]
-    """Element-wise libm call (bit-identical to C)."""
+    """Element-wise libm call (bit-identical to C); compiled when ``simurg_core`` is installed."""
+    if core is not None:
+        xa = np.ascontiguousarray(x, dtype=np.float64).ravel()
+        if (e := getattr(fn, "exponent", None)) is not None:
+            return core.pow(xa, e).reshape(np.shape(x))  # type: ignore[no-any-return]
+        if (name := _CORE_FN.get(fn)) is not None:
+            return getattr(core, name)(xa).reshape(np.shape(x))  # type: ignore[no-any-return]
    return np.fromiter((_c_call(fn, v) for v in x.tolist()), dtype=np.float64, count=len(x))


@@ -51,6 +61,7 @@ def _pow(e: float):  # type: ignore[no-untyped-def]
    def f(v: float) -> float:
        return math.pow(v, e)

+    f.exponent = e  # type: ignore[attr-defined]
    return f


@@ -137,6 +148,23 @@ def call_genotypes(
    cutoff: float = DEFAULT_CUTOFF,
 ) -> Calls:
    """Bayesian calls and GenCall scores for one sample (plan §7.5)."""
+    if core is not None:
+        f32 = [np.ascontiguousarray(a, np.float32).ravel() for a in (c.theta_mean, c.theta_dev, c.r_mean, c.r_dev)]
+        geno, score = core.call_genotypes(
+            np.ascontiguousarray(theta, np.float32),
+            np.ascontiguousarray(r, np.float32),
+            np.ascontiguousarray(raw_zero, np.bool_),
+            np.ascontiguousarray(c.present, np.bool_),
+            np.ascontiguousarray(c.n, np.int32).ravel(),
+            *f32,
+            np.ascontiguousarray(c.intensity_threshold, np.float32),
+            np.ascontiguousarray(c.gentrain, np.float32),
+            np.ascontiguousarray(intensity_only, np.bool_),
+            F32(cutoff),
+        )
+        return Calls(
+            genotypes=geno, scores=score, num_calls=int((geno > 0).sum()), num_intensity_only=int(intensity_only.sum())
+        )
    n = len(theta)
    t = np.asarray(theta, np.float32)
    rr = np.asarray(r, np.float32)
@@ -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"]
@@ -160,7 +169,7 @@ def _process(job: Job) -> dict[str, Any]:
    try:
        grn, red = read_idat(job.grn), read_idat(job.red)
        check_pair(grn, red)
-        _, notes = check_coverage(grn, ctx.manifest)
+        _, notes = check_coverage(grn, ctx.manifest, ctx.cache)
        row["notes"].extend(notes)
        g = genotype_sample(grn, red, ctx, params, _W["started"])
        write_gtc(g, Path(job.gtc_path))
@@ -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,
-                mp_context=mp.get_context("spawn"),
-                initializer=_init,
-                initargs=(ctx, params, auto_1_22, probes, started, beads, th.bg_offset),
+                    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")
@@ -14,6 +14,7 @@ from dataclasses import dataclass
 import numpy as np
 import numpy.typing as npt

+from simurg._core import core
 from simurg.signal.normalization.numerics import (
    F32,
    TINY,
@@ -86,6 +87,8 @@ def closest_points(xref: FloatArr, yref: FloatArr, x: FloatArr, y: FloatArr) ->

 def _closest_slow(xref: FloatArr, yref: FloatArr, x: FloatArr, y: FloatArr) -> IdxArr:
    """Exact search: squared distances in float32; the first minimum wins (C strict ``<``)."""
+    if core is not None and len(xref) and all(a.dtype == np.float32 for a in (xref, yref, x, y)):
+        return core.closest_exact(*(np.ascontiguousarray(a) for a in (xref, yref, x, y)))  # type: ignore[no-any-return]
    dx = x[:, None] - xref[None, :]
    dy = y[:, None] - yref[None, :]
    return np.argmin(dx * dx + dy * dy, axis=1).astype(np.int64)
@@ -118,6 +121,12 @@ def _closest_binned(raw_x: FloatArr, raw_y: FloatArr, axis_x: FloatArr, axis_y:
    low = rem < F32(0.25)
    choice = np.where(mid, k, np.where(low, np.where(k == 0, k, 12 + k - 1), np.where(k == 11, k, 12 + k)))

+    if core is not None and all(a.dtype == np.float32 for a in (raw_x, raw_y, axis_x, axis_y)):
+        i64 = np.int64
+        return core.closest_binned(  # type: ignore[no-any-return]
+            *(np.ascontiguousarray(a, np.float32) for a in (raw_x, raw_y, axis_x, axis_y)),
+            *(np.ascontiguousarray(a, i64) for a in (choice, eligible, b1, b2)),
+        )
    out = np.full(len(axis_a), -1, dtype=np.int64)
    for c in np.unique(choice):
        cand = eligible[b1 == c] if c < 12 else eligible[b2 == c - 12]
@@ -90,9 +90,11 @@ def address_coverage(idat: Idat, manifest: Manifest) -> float:
    return float(found.mean())


-def check_coverage(idat: Idat, manifest: Manifest) -> tuple[float, list[tuple[str, str]]]:
+def check_coverage(
+    idat: Idat, manifest: Manifest, cache: LocusMapCache | None = None
+) -> tuple[float, list[tuple[str, str]]]:
    """IN-027: < 0.99 → E_IDAT_MANIFEST_MISMATCH (raised); < 1 → W_IDAT_LOW_ADDRESS_COVERAGE (returned)."""
-    cov = address_coverage(idat, manifest)
+    cov = cache.coverage(idat) if cache is not None and cache.manifest is manifest else address_coverage(idat, manifest)
    notes: list[tuple[str, str]] = []
    if cov < COVERAGE_ERROR:
        raise SimurgError(
@@ -149,10 +151,24 @@ class LocusMapCache:
        self.manifest = manifest
        self.caddr = control_addresses(manifest)
        self._maps: dict[bytes, LocusMap] = {}
+        self._coverage: dict[bytes, float] = {}
        self.max_entries = max_entries

+    @staticmethod
+    def _key(idat: Idat) -> bytes:
+        return hashlib.blake2b(np.ascontiguousarray(idat.address_ids).tobytes(), digest_size=16).digest()
+
+    def coverage(self, idat: Idat) -> float:
+        """``address_coverage`` memoised per address layout."""
+        key = self._key(idat)
+        if (cov := self._coverage.get(key)) is None:
+            if len(self._coverage) >= self.max_entries:
+                self._coverage.pop(next(iter(self._coverage)))
+            cov = self._coverage[key] = address_coverage(idat, self.manifest)
+        return cov
+
    def get(self, idat: Idat) -> LocusMap:
-        key = hashlib.blake2b(np.ascontiguousarray(idat.address_ids).tobytes(), digest_size=16).digest()
+        key = self._key(idat)
        lm = self._maps.get(key)
        if lm is None:
            if len(self._maps) >= self.max_entries:
@@ -93,3 +93,18 @@ def test_locale_independent_output(tmp_path) -> None:  # type: ignore[no-untyped
        assert r.returncode == 0, r.stderr
        outs.add(r.stdout)
    assert len(outs) == 1
+
+
+def test_closed_stdout_is_not_an_internal_error(tmp_path) -> None:  # type: ignore[no-untyped-def]
+    # `simurg inspect … | head`: the reader goes away; the CLI must exit like a SIGPIPE'd tool (141), silently
+    f = _files(tmp_path)
+    p = subprocess.Popen(
+        [sys.executable, "-m", "simurg.cli.main", "inspect", "--json", str(f["bpm"]), str(f["egt"]), str(f["idat"])],
+        stdout=subprocess.PIPE,
+        stderr=subprocess.PIPE,
+    )
+    assert p.stdout is not None and p.stderr is not None
+    p.stdout.close()  # nobody reads stdout
+    err = p.stderr.read().decode()
+    assert p.wait() == 141, err
+    assert "E_INTERNAL" not in err and "Traceback" not in err and "Exception ignored" not in err
@@ -0,0 +1,238 @@
+"""The compiled kernels (``rust/simurg-core``, plan §18 P9) must equal the pure-Python paths bit for bit."""
+
+from __future__ import annotations
+
+import datetime as dt
+import math
+from dataclasses import fields, replace
+
+import numpy as np
+import pytest
+from tests.conftest import GSA_BPM, GSA_EGT, GSA_IDATS
+
+import simurg.formats.bpm as bpm_mod
+import simurg.genotype.baf_lrr as baf_mod
+import simurg.genotype.gencall as gencall_mod
+import simurg.signal.normalization.fit as fit_mod
+from simurg._core import core
+from simurg.core.errors import FormatError
+from simurg.formats.bpm import read_bpm
+from simurg.synth.files import SynthLocus, encode_bpm
+
+pytestmark = pytest.mark.skipif(core is None, reason="simurg_core extension not installed")
+
+MODULES = (bpm_mod, gencall_mod, fit_mod, baf_mod)
+
+
+@pytest.fixture
+def pure(monkeypatch):  # type: ignore[no-untyped-def]
+    """Context switch: ``with pure(): ...`` runs the pure-Python paths."""
+
+    class _Pure:
+        def __enter__(self) -> None:
+            for m in MODULES:
+                monkeypatch.setattr(m, "core", None)
+
+        def __exit__(self, *exc: object) -> None:
+            for m in MODULES:
+                monkeypatch.setattr(m, "core", core)
+
+    return _Pure
+
+
+def _bits(a: np.ndarray) -> bytes:
+    return np.ascontiguousarray(a).tobytes()
+
+
+def _edge_values() -> np.ndarray:
+    neg_nan = np.frombuffer(np.uint64(0xFFF8000000000000).tobytes(), np.float64)[0]
+    payload_nan = np.frombuffer(np.uint64(0x7FF8000000001234).tobytes(), np.float64)[0]
+    tiny = np.finfo(np.float64).tiny
+    special = [0.0, -0.0, 1.0, -1.0, math.inf, -math.inf, math.nan, neg_nan, payload_nan, tiny, tiny / 4, -tiny,
+               1e-320, 709.78, 709.79, -745.2, -746.0, 1e300, -1e300, 2.0**-1074, 0.5, 10.0, 1e-5]  # fmt: skip
+    rng = np.random.default_rng(7)
+    rand = np.concatenate(
+        [rng.normal(0, 3, 20_000), rng.uniform(-800, 800, 5_000), np.exp(rng.uniform(-700, 700, 5_000))]
+    )
+    return np.concatenate([np.array(special, np.float64), rand, -rand[:1000]])
+
+
+@pytest.mark.parametrize("fn", [math.exp, math.log10, math.log2, gencall_mod._pow(1.71), gencall_mod._pow(0.35)])
+def test_libm_kernels_bit_identical(fn, pure) -> None:  # type: ignore[no-untyped-def]
+    x = _edge_values()
+    fast = gencall_mod._libm(fn, x)
+    with pure():
+        ref = gencall_mod._libm(fn, x)
+    assert fast.dtype == ref.dtype == np.float64
+    assert _bits(fast) == _bits(ref)
+
+
+def test_pow_kernel_rejects_integer_exponent() -> None:
+    assert core is not None
+    with pytest.raises(ValueError):
+        core.pow(np.ones(3), 2.0)
+
+
+def test_closest_binned_identical(pure) -> None:  # type: ignore[no-untyped-def]
+    rng = np.random.default_rng(3)
+    for trial in range(20):
+        n = int(rng.integers(0, 3000))
+        raw_x = rng.gamma(2.0, 400.0, n).astype(np.float32)
+        raw_y = rng.gamma(1.0, 60.0 if trial % 2 else 3000.0, n).astype(np.float32)
+        if trial == 5 and n:
+            raw_x[: n // 3] = raw_x[0]  # ties: first minimum must win
+            raw_y[: n // 3] = raw_y[0]
+        if trial == 7 and n:
+            raw_x[n // 2] = np.nan
+        m = int(rng.integers(1, 400))
+        axis = np.linspace(0, float(np.float32(rng.uniform(100, 5000))), m).astype(np.float32)
+        zeros = np.zeros(m, np.float32)
+        ax, ay = (axis, zeros) if trial % 3 else (zeros, axis)
+        fast = fit_mod._closest_binned(raw_x, raw_y, ax, ay)
+        with pure():
+            ref = fit_mod._closest_binned(raw_x, raw_y, ax, ay)
+        assert fast.dtype == ref.dtype
+        assert np.array_equal(fast, ref), trial
+
+
+def _loci() -> list[SynthLocus]:
+    return [
+        SynthLocus("rs1", "1", 100, "[A/G]", 11, 0, 0, norm_id=1, source_seq="ACGT[A/G]TTGA", probe_a="ACGTA"),
+        SynthLocus("rs2", "2", 200, "[A/T]", 12, 13, 1, ilmn_strand="BOT", ref_strand="-", probe_a="AC", probe_b="AT"),
+        SynthLocus("rs3", "X", 300, "[C/G]", 14, 15, 2, norm_id=99),
+        SynthLocus("rsü", "Y", 400, "[T/C]", 16, 0, 0, norm_id=2, intensity_only=1),
+        SynthLocus("rs5", "MT", 500, "[A/C]", 17, 0, 0, norm_id=99, source_seq="A" * 200 + "[A/C]"),
+    ]
+
+
+def _manifest_dict(m) -> dict:  # type: ignore[no-untyped-def]
+    out = {}
+    for f in fields(m):
+        v = getattr(m, f.name)
+        out[f.name] = (v.dtype.str, v.shape, v.tobytes()) if isinstance(v, np.ndarray) else v
+    return out
+
+
+@pytest.mark.parametrize("entry_version", [4, 6, 7, 8])
+@pytest.mark.parametrize("order", [None, [4, 3, 2, 1, 0]])
+def test_bpm_parser_identical(tmp_path, pure, entry_version, order) -> None:  # type: ignore[no-untyped-def]
+    p = tmp_path / "m.bpm"
+    p.write_bytes(encode_bpm(_loci(), entry_version=entry_version, entry_order=order))
+    fast = _manifest_dict(read_bpm(p))
+    with pure():
+        ref = _manifest_dict(read_bpm(p))
+    assert fast == ref
+
+
+def test_bpm_parser_errors_identical(tmp_path, pure) -> None:  # type: ignore[no-untyped-def]
+    good = encode_bpm(_loci())
+    cases = {
+        "truncated": good[: len(good) - 40],
+        "bad_assay": encode_bpm([SynthLocus("rs1", address_b=0, assay_type=1)]),
+        "bad_map": encode_bpm([SynthLocus("rs1", pos=12345)]).replace(b"\x0512345", b"\x0512x45"),
+    }
+    for label, data in cases.items():
+        p = tmp_path / f"{label}.bpm"
+        p.write_bytes(data)
+        with pytest.raises(FormatError) as fast:
+            read_bpm(p)
+        with pure(), pytest.raises(FormatError) as ref:
+            read_bpm(p)
+        assert (fast.value.code, str(fast.value)) == (ref.value.code, str(ref.value)), label
+
+
+@pytest.mark.data
+@pytest.mark.skipif(not GSA_BPM.exists(), reason="GSA demo data not fetched")
+def test_genotype_sample_identical(pure) -> None:  # type: ignore[no-untyped-def]
+    from simurg.formats.egt import read_egt
+    from simurg.formats.gtc import encode_gtc
+    from simurg.formats.idat import read_idat
+    from simurg.genotype.clusters import align_clusters
+    from simurg.genotype.sample import genotype_sample, manifest_context, params_for_preset
+    from simurg.signal.raw_xy import LocusMapCache
+
+    def run() -> bytes:
+        m = read_bpm(GSA_BPM)
+        ctx = manifest_context(m, align_clusters(read_egt(GSA_EGT), m))
+        ctx = replace(ctx, manifest=m.lite(), cache=LocusMapCache(m.lite()))
+        stem = next(iter(sorted(GSA_IDATS.rglob("*_Grn.idat")))).as_posix()[: -len("_Grn.idat")]
+        g = genotype_sample(read_idat(stem + "_Grn.idat"), read_idat(stem + "_Red.idat"), ctx,
+                            params_for_preset("current"), dt.datetime(2026, 1, 1))  # fmt: skip
+        return encode_gtc(g)
+
+    fast = run()
+    with pure():
+        ref = run()
+    assert fast == ref
+
+
+def test_closest_exact_identical(pure) -> None:  # type: ignore[no-untyped-def]
+    rng = np.random.default_rng(11)
+    for trial in range(10):
+        n = int(rng.integers(1, 1500))
+        xr = rng.normal(500, 200, n).astype(np.float32)
+        yr = rng.normal(300, 100, n).astype(np.float32)
+        if trial == 2:
+            xr[n // 2 :] = xr[0]
+            yr[n // 2 :] = yr[0]
+        if trial == 3:
+            xr[n // 3] = np.nan
+        m = int(rng.integers(1, 300))
+        x = rng.normal(500, 200, m).astype(np.float32)
+        y = rng.normal(300, 100, m).astype(np.float32)
+        fast = fit_mod._closest_slow(xr, yr, x, y)
+        with pure():
+            ref = fit_mod._closest_slow(xr, yr, x, y)
+        assert np.array_equal(fast, ref), trial
+
+
+def _synthetic_clusters(n: int, rng: np.random.Generator):  # type: ignore[no-untyped-def]
+    from simurg.genotype.clusters import ClusterArrays
+
+    tm = np.sort(rng.uniform(0, 1, (n, 3)), axis=1).astype(np.float32)
+    td = rng.uniform(0.0, 0.08, (n, 3)).astype(np.float32)
+    rm = rng.uniform(0.2, 3, (n, 3)).astype(np.float32)
+    rd = rng.uniform(0.0, 0.3, (n, 3)).astype(np.float32)
+    counts = rng.integers(0, 200, (n, 3)).astype(np.int32)
+    counts[rng.random((n, 3)) < 0.2] = 0
+    counts[: n // 50] = 0  # no samples at all
+    td[rng.random((n, 3)) < 0.02] = 0.0
+    tm[rng.random((n, 3)) < 0.01] = np.nan
+    rm[rng.random((n, 3)) < 0.01] = np.nan
+    thr = rng.uniform(0, 0.3, n).astype(np.float32)
+    thr[rng.random(n) < 0.3] = np.nan
+    gt = rng.uniform(0, 1, n).astype(np.float32)
+    gt[rng.random(n) < 0.02] = np.nan
+    return ClusterArrays(
+        present=rng.random(n) > 0.01, n=counts, theta_mean=tm, theta_dev=td, r_mean=rm, r_dev=rd,
+        intensity_threshold=thr, gentrain=gt, mean_r=rng.uniform(0.2, 3, n).astype(np.float32),
+    )  # fmt: skip
+
+
+def test_call_genotypes_and_baf_lrr_identical(pure) -> None:  # type: ignore[no-untyped-def]
+    from simurg.genotype.baf_lrr import baf_lrr
+    from simurg.genotype.gencall import call_genotypes
+
+    rng = np.random.default_rng(5)
+    n = 60_000
+    c = _synthetic_clusters(n, rng)
+    theta = rng.uniform(-0.05, 1.05, n).astype(np.float32)
+    theta[: n // 3] = c.theta_mean[: n // 3, rng.integers(0, 3)]  # exactly on a cluster centre
+    r = rng.gamma(2.0, 0.6, n).astype(np.float32)
+    raw_zero = rng.random(n) < 0.01
+    theta[raw_zero] = np.nan
+    r[raw_zero] = np.nan
+    r[rng.random(n) < 0.01] = 0.0
+    intensity_only = rng.random(n) < 0.02
+    for cutoff in (0.15, 0.0):
+        fast = call_genotypes(theta, r, raw_zero, c, intensity_only, cutoff)
+        with pure():
+            ref = call_genotypes(theta, r, raw_zero, c, intensity_only, cutoff)
+        assert _bits(fast.genotypes) == _bits(ref.genotypes)
+        assert _bits(fast.scores) == _bits(ref.scores)
+        assert (fast.num_calls, fast.num_intensity_only) == (ref.num_calls, ref.num_intensity_only)
+    fb, fl = baf_lrr(theta, r, raw_zero, c, intensity_only)
+    with pure():
+        rb, rl = baf_lrr(theta, r, raw_zero, c, intensity_only)
+    assert _bits(fb) == _bits(rb)
+    assert _bits(fl) == _bits(rl)
Author	SHA1	Message	Date
kemal	d235dc5419	Release 1.1.0: compiled kernels (P9) and packaging (P10) ci / test (3.13) (push) Failing after 12s Details ci / test (3.12) (push) Failing after 9s Details	2026-09-28 23:56:39 +01:00
kemal	0b05b0d346	P10 complete: packaging verified on fresh x86_64 and aarch64 VMs Wheels (pure + manylinux_2_28 kernels), sdist, Docker image and conda packages built and installed on fresh Ubuntu 24.04 machines; 225 unit tests pass against the installed package; the tutorial runs with pip (with and without kernels); GTC/VCF content digests are identical across pip, Docker and conda on both architectures and equal to the Pi's run. docs/validation/P10.md; guide: install from wheels, Docker, conda; exit code 141.	2026-09-28 23:53:59 +01:00
kemal	522429b731	CLI: a closed stdout pipe exits 141 silently (was E_INTERNAL); example avoids head under pipefail Found by the P10 fresh-VM run: `simurg inspect … \| head` raised BrokenPipeError, reported as an internal error, and stopped examples/gsa_demo.sh (set -o pipefail). run_guarded now flushes stdout inside the guard and maps BrokenPipeError to ExitCode.SIGPIPE (141) without a message. p10.sh: SKIP_DOCKER/SKIP_CONDA, clearer digest report; p10t.sh runs build+pip+tutorial only.	2026-09-28 23:14:47 +01:00
kemal	035e5ad66f	P10 (WIP): release build script, Dockerfile, conda recipes, fresh-VM acceptance task scripts/build_dist.sh: sdist + pure wheel + manylinux_2_28 kernel wheels (maturin --zig, both architectures). Dockerfile (multi-stage, kernels included, non-root user). packaging/conda/{simurg,simurg-core}. scripts/p10/digest.py: platform-independent content digests of GTC/VCF/summary outputs. Cloud task p10 installs the artifacts on a fresh Ubuntu VM and runs the tutorial with pip, Docker and conda installs. trial.sh picks the arm64/amd64 image from the instance type and can run instances in parallel (PARALLEL=1). Licence files ship in both wheels; Rust crate notices added.	2026-09-28 22:41:56 +01:00
kemal	4d5c9be4ae	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.	2026-09-28 22:32:48 +01:00
kemal	d40da1d670	P9: scaling on a RAM disk (bench --work/--clean, cloud task p9b)	2026-09-28 22:04:10 +01:00
kemal	f0d3739493	P9 cloud task: install gcc (linker for the Rust build); guide: kernel installation	2026-09-28 21:35:42 +01:00
kemal	587fadbaef	bench: format	2026-09-28 21:29:07 +01:00
kemal	95e2338897	P9 (WIP): compiled kernels (rust/simurg-core, PyO3) with bit-identical pure-Python fallback Rust ports of the hot paths found by profiling: fused per-locus GenCall calling and scoring, BAF/LRR, libm maps with CPython math semantics, both nearest-neighbour searches of the normalization fit, and the BPM LocusEntry parser. Address coverage is memoised per IDAT layout. SIMURG_PURE_PYTHON=1 forces the pure paths. On the Pi: genotype_sample 3.4 s -> 1.1 s, 36-sample GSA run 131 s -> 46 s, BPM read 4x faster; 36 GSA + 6 Omni2.5 GTCs and QC tables byte-identical to the pure path; tests/unit/test_core_kernels.py checks every kernel on edge cases. bench/profile_genotype.py, bench/scaling.py, cloud task p9 (thread scaling).	2026-09-28 21:28:55 +01:00
				`@@ -0,0 +1 @@`
				`{"gtcs": 360, "differ": 0, "sample_summary": "identical"}`