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

@@ -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)

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) ==="

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",
+                )
             }
         )
     )

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