From 17bd0e8f78e3000a2eb1fb1f585b0ac7fa20b0e9 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Sun, 27 Sep 2026 20:57:42 +0100 Subject: [PATCH] P2 validation: oracle test and large-scale raw X/Y check - tests/oracle/test_p2_raw_xy.py: IDAT-assembled raw X/Y and controls equal the arrays in oracle GTCs - validate_p1_large.py: phase-P2 comparison against vendor GTCs with diagnostics (channel swap, bead-count disagreements) - cloud trial: also runs `simurg ingest` on all samples (time/memory), disk 250 GB --- scripts/cloud/trial.sh | 2 +- scripts/cloud/userdata.sh | 12 ++++- scripts/validate_p1_large.py | 96 ++++++++++++++++++++++++++++++++-- tests/oracle/test_p2_raw_xy.py | 36 +++++++++++++ 4 files changed, 141 insertions(+), 5 deletions(-) create mode 100644 tests/oracle/test_p2_raw_xy.py diff --git a/scripts/cloud/trial.sh b/scripts/cloud/trial.sh index d8f7dfe..6df63a8 100755 --- a/scripts/cloud/trial.sh +++ b/scripts/cloud/trial.sh @@ -16,7 +16,7 @@ REGION=${REGION:-us-east-1} TYPE=${TYPE:-m7g.4xlarge} 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:-200} +DISK_GB=${DISK_GB:-250} BUDGET_USD=${BUDGET_USD:-10} TAG=simurg-trial ACCOUNT=$(aws sts get-caller-identity --query Account --output text) diff --git a/scripts/cloud/userdata.sh b/scripts/cloud/userdata.sh index 5f51696..11a9848 100755 --- a/scripts/cloud/userdata.sh +++ b/scripts/cloud/userdata.sh @@ -33,7 +33,7 @@ echo "=== start $(date -u +%FT%TZ) on $(uname -m), $(nproc) cpus ===" export DEBIAN_FRONTEND=noninteractive apt-get update -q -apt-get install -y -q pigz git python3-venv curl unzip +apt-get install -y -q pigz git python3-venv curl unzip time snap install aws-cli --classic S3="/snap/bin/aws s3" ( while sleep 300; do upload_log; done ) & @@ -60,4 +60,14 @@ echo "=== validation $(date -u +%T) ===" --egt $D/HumanOmni2.5-4v1-Multi_B.egt --csv $D/HumanOmni2.5-4v1_B.csv --out /work/summary.json \ --workers "$(nproc)" || echo "validation reported differences (see summary.json)" $S3 cp /work/summary.json "s3://$BUCKET/$RUN/summary.json" --only-show-errors + +echo "=== simurg ingest at scale $(date -u +%T) ===" +/usr/bin/time -v .venv/bin/simurg ingest --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm --idats $D/data \ + --out /work/ingest --threads "$(nproc)" --quiet 2> /work/ingest_time.txt || echo "ingest exit code $?" +tail -25 /work/ingest_time.txt +for f in ingest_summary.tsv run_manifest.json; do + $S3 cp "/work/ingest/$f" "s3://$BUCKET/$RUN/ingest_$f" --only-show-errors || true +done +$S3 cp /work/ingest_time.txt "s3://$BUCKET/$RUN/ingest_time.txt" --only-show-errors || true +echo "ingest output size: $(du -sh /work/ingest | cut -f1)" echo "=== done $(date -u +%T) ===" diff --git a/scripts/validate_p1_large.py b/scripts/validate_p1_large.py index 31bbe6a..8cfb044 100644 --- a/scripts/validate_p1_large.py +++ b/scripts/validate_p1_large.py @@ -103,6 +103,59 @@ def check_idat_pair(args: tuple[str, str]) -> dict[str, Any]: return res +_P2: dict[str, Any] = {} + + +def _p2_init(manifest: Any) -> None: + _P2.clear() + _P2["m"] = manifest + + +def check_raw_xy(args: tuple[str, str, str]) -> dict[str, Any]: + """P2: raw X/Y + controls assembled from IDATs vs the arrays stored in the GTC (TOC 1000/1001/500/501).""" + from simurg.formats.gtc import read_gtc + from simurg.formats.idat import read_idat + from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy + + gtc, grn, red = args + res: dict[str, Any] = {"gtc": gtc, "ok": True} + try: + m = _P2["m"] + cache = _P2.setdefault("cache", LocusMapCache(m)) + ref = read_gtc(gtc, check_write_complete=False) + raw = assemble_raw_xy(read_idat(grn), read_idat(red), m, cache=cache) + dx = raw.raw_x != ref.raw_x + dy = raw.raw_y != ref.raw_y + diff = dx | dy + res.update( + n_diff=int(diff.sum()), + n_missing=raw.n_missing, + vendor_zero=int(((ref.raw_x == 0) & (ref.raw_y == 0)).sum()), + # plan ยง20 U1: would swapping channels explain the differences? + swapped=int(((raw.raw_x == ref.raw_y) & (raw.raw_y == ref.raw_x) & diff).sum()), + ours_missing_vendor_nonzero=int((raw.missing & ((ref.raw_x != 0) | (ref.raw_y != 0))).sum()), + vendor_zero_ours_present=int((~raw.missing & (ref.raw_x == 0) & (ref.raw_y == 0)).sum()), + controls_equal=bool( + np.array_equal(raw.controls_x, ref.controls_x) and np.array_equal(raw.controls_y, ref.controls_y) + ), + ) + if res["n_diff"]: + idx = np.flatnonzero(diff)[:5].tolist() + res["examples"] = [ + { + "locus": m.names[i], + "assay": int(m.assay_types[i]), + "ours": [int(raw.raw_x[i]), int(raw.raw_y[i])], + "vendor": [int(ref.raw_x[i]), int(ref.raw_y[i])], + } + for i in idx + ] + res["ok"] = res["n_diff"] == 0 and res["controls_equal"] + except Exception as e: + res.update(ok=False, error=f"{type(e).__name__}: {e}") + return res + + def main() -> int: ap = argparse.ArgumentParser() ap.add_argument("--root", required=True, type=Path) @@ -112,6 +165,7 @@ def main() -> int: ap.add_argument("--out", required=True, type=Path) ap.add_argument("--workers", type=int, default=os.cpu_count() or 1) ap.add_argument("--limit", type=int, default=0) + ap.add_argument("--no-p2", action="store_true", help="skip the phase-P2 raw X/Y check") a = ap.parse_args() from simurg.formats.bpm import read_bpm @@ -164,6 +218,28 @@ def main() -> int: "normalization_lookups": m.normalization_lookups().tolist() == b.normalization_lookups, "n_bins": len(np.unique(m.norm_bin_ids())), } + if not a.no_p2: + import multiprocessing as mp + import re + + key_re = re.compile(r"(\d+_r\d{2}c\d{2})", re.IGNORECASE) + triples = [] + for g in gtcs: + mk = key_re.search(Path(g).name) + k = mk.group(1).lower() if mk else Path(g).stem.lower() + if k in grn and k in red: + triples.append((g, str(grn[k]), str(red[k]))) + summary["p2_matched_gtc_idat"] = len(triples) + with ProcessPoolExecutor( + a.workers, mp_context=mp.get_context("spawn"), initializer=_p2_init, initargs=(m.lite(),) + ) as ex: + p2 = list(ex.map(check_raw_xy, triples, chunksize=4)) + summary["p2_ok"] = sum(r["ok"] for r in p2) + summary["p2_total_diff_loci"] = sum(r.get("n_diff", 0) for r in p2) + summary["p2_swapped_loci"] = sum(r.get("swapped", 0) for r in p2) + summary["p2_ours_missing_vendor_nonzero"] = sum(r.get("ours_missing_vendor_nonzero", 0) for r in p2) + summary["p2_vendor_zero_ours_present"] = sum(r.get("vendor_zero_ours_present", 0) for r in p2) + summary["p2_failures"] = [r for r in p2 if not r["ok"]][:30] e = read_egt(a.egt) with open(a.egt, "rb") as fh: be = module.ClusterFile.read_cluster_file(fh) @@ -200,7 +276,10 @@ def main() -> int: summary["gtc_transforms_vs_bins"] = [len(g0.normalization_transforms), summary["bpm"]["n_bins"]] summary["seconds"] = round(time.time() - t0, 1) - passed = ( + p2_passed = a.no_p2 or ( + summary.get("p2_ok") == summary.get("p2_matched_gtc_idat", -1) and summary.get("p2_matched_gtc_idat", 0) > 0 + ) + passed = p2_passed and ( summary["gtc_ok"] == len(gtcs) and summary["idat_ok"] == len(pairs) and all(v for k, v in summary["bpm"].items() if isinstance(v, bool)) @@ -213,8 +292,19 @@ def main() -> int: print( json.dumps( { - k: summary[k] - for k in ("n_gtc", "gtc_ok", "gtc_roundtrip_identical", "n_idat_pairs", "idat_ok", "passed", "seconds") + k: summary.get(k) + for k in ( + "n_gtc", + "gtc_ok", + "gtc_roundtrip_identical", + "n_idat_pairs", + "idat_ok", + "p2_matched_gtc_idat", + "p2_ok", + "p2_total_diff_loci", + "passed", + "seconds", + ) } ) ) diff --git a/tests/oracle/test_p2_raw_xy.py b/tests/oracle/test_p2_raw_xy.py new file mode 100644 index 0000000..2642bb4 --- /dev/null +++ b/tests/oracle/test_p2_raw_xy.py @@ -0,0 +1,36 @@ +"""Phase P2 acceptance on P0 data: raw X/Y and control intensities assembled from IDATs equal the arrays stored +in oracle GTCs (TOC 1000/1001, 500/501). The same check on vendor GTCs runs in scripts/validate_p1_large.py --p2. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from simurg.formats.bpm import read_bpm +from simurg.formats.gtc import read_gtc +from simurg.formats.idat import read_idat +from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy, check_coverage +from tests.conftest import GSA_BPM, GSA_IDATS, ORACLE_GTCS, need + +pytestmark = pytest.mark.data + + +def test_raw_xy_equals_oracle_gtcs() -> None: + need(GSA_BPM, GSA_IDATS, ORACLE_GTCS) + m = read_bpm(GSA_BPM) + cache = LocusMapCache(m) + gtcs = sorted(ORACLE_GTCS.glob("*.gtc")) + assert gtcs + for g in gtcs: + sid = g.stem + chip = sid.split("_")[0] + grn = read_idat(GSA_IDATS / chip / f"{sid}_Grn.idat") + red = read_idat(GSA_IDATS / chip / f"{sid}_Red.idat") + assert check_coverage(grn, m)[0] == 1.0 + raw = assemble_raw_xy(grn, red, m, cache=cache) + ref = read_gtc(g) + assert np.array_equal(raw.raw_x, ref.raw_x), sid + assert np.array_equal(raw.raw_y, ref.raw_y), sid + assert np.array_equal(raw.controls_x, ref.controls_x), sid + assert np.array_equal(raw.controls_y, ref.controls_y), sid