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 2583 additions and 28 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
+worked example: `examples/gsa_demo.sh`. P9 (performance) is done: optional Rust kernels, byte-identical outputs,
-modules (cluster training, CN/LOH, PGx).
+~80 GSA samples/min per server core (`docs/validation/P9.md`). 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
+* Speed: with the compiled kernels a GSA sample takes ~1.1 s on a Raspberry Pi 5 core (pure Python ~3.4 s); see
-  16 cores) and ~25 s on a Raspberry Pi 5 core; a compiled core is planned (P9).
+  `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:
-            with ProcessPoolExecutor(
+            opts.out.mkdir(parents=True, exist_ok=True)
-                workers,
+            fd, state = tempfile.mkstemp(prefix=".simurg-workers-", suffix=".pkl", dir=opts.out)
-                mp_context=mp.get_context("spawn"),
+            try:
-                initializer=_init,
+                with os.fdopen(fd, "wb") as fh:
-                initargs=(ctx, params, auto_1_22, probes, started, beads, th.bg_offset),
+                    pickle.dump((ctx, params, auto_1_22, probes, started, beads, th.bg_offset), fh, protocol=5)
-            ) as ex:
+                with ProcessPoolExecutor(
-                results = list(ex.map(_process, jobs))
+                    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"}`