diff --git a/.gitignore b/.gitignore index 3f56cfc..b91a2b1 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,5 @@ __pycache__/ .ruff_cache/ # rust target/ +.trial-last +results/ diff --git a/scripts/cloud/trial.sh b/scripts/cloud/trial.sh new file mode 100755 index 0000000..d8f7dfe --- /dev/null +++ b/scripts/cloud/trial.sh @@ -0,0 +1,123 @@ +#!/usr/bin/env bash +# Self-destructing, budget-guarded AWS validation run for SimurgArray. +# +# bash scripts/cloud/trial.sh setup # one-time: budget alert, results bucket, instance role +# bash scripts/cloud/trial.sh launch # start one spot run (terminates itself when done) +# bash scripts/cloud/trial.sh status # instance state + last lines of the remote log +# bash scripts/cloud/trial.sh results # download log/summary of the latest run +# bash scripts/cloud/trial.sh kill # emergency: terminate every trial instance now +# bash scripts/cloud/trial.sh audit # list anything still costing money +# +# Cost guards: spot price cap, hard time cap (shutdown = terminate), EXIT trap shutdown, CloudWatch idle alarm +# that terminates from outside, disk deleted with the instance, monthly budget e-mail alert. +set -euo pipefail +cd "$(dirname "$0")/../.." +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} +BUDGET_USD=${BUDGET_USD:-10} +TAG=simurg-trial +ACCOUNT=$(aws sts get-caller-identity --query Account --output text) +BUCKET=simurg-trial-$ACCOUNT +ROLE=simurg-trial-role +AWS="aws --region $REGION" + +cmd=${1:-help} +case "$cmd" in +setup) + EMAIL=${2:?usage: trial.sh setup } + echo "-> budget alert: \$$BUDGET_USD/month, e-mail $EMAIL at 50%, 100% actual and 100% forecast" + cat > /tmp/simurg-budget.json < /tmp/simurg-notif.json + aws budgets create-budget --account-id "$ACCOUNT" --budget file:///tmp/simurg-budget.json \ + --notifications-with-subscribers file:///tmp/simurg-notif.json 2>/dev/null || echo " (budget already exists)" + echo "-> results bucket s3://$BUCKET (objects expire after 30 days)" + $AWS s3api create-bucket --bucket "$BUCKET" >/dev/null 2>&1 || true + $AWS s3api put-public-access-block --bucket "$BUCKET" --public-access-block-configuration \ + BlockPublicAcls=true,IgnorePublicAcls=true,BlockPublicPolicy=true,RestrictPublicBuckets=true + $AWS s3api put-bucket-lifecycle-configuration --bucket "$BUCKET" --lifecycle-configuration \ + '{"Rules":[{"ID":"expire","Status":"Enabled","Filter":{},"Expiration":{"Days":30}}]}' + echo "-> instance role $ROLE (read/write only this bucket)" + aws iam create-role --role-name $ROLE --assume-role-policy-document \ + '{"Version":"2012-10-17","Statement":[{"Effect":"Allow","Principal":{"Service":"ec2.amazonaws.com"},"Action":"sts:AssumeRole"}]}' \ + >/dev/null 2>&1 || true + aws iam put-role-policy --role-name $ROLE --policy-name bucket-only --policy-document \ + "{\"Version\":\"2012-10-17\",\"Statement\":[{\"Effect\":\"Allow\",\"Action\":[\"s3:GetObject\",\"s3:PutObject\"],\"Resource\":\"arn:aws:s3:::$BUCKET/*\"}]}" + aws iam create-instance-profile --instance-profile-name $ROLE >/dev/null 2>&1 || true + aws iam add-role-to-instance-profile --instance-profile-name $ROLE --role-name $ROLE 2>/dev/null || true + echo "setup done (budget alerts need you to confirm nothing; AWS e-mails $EMAIL when thresholds are crossed)" + ;; +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; } + 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 + $AWS s3 cp /tmp/simurg-repo.tar.gz "s3://$BUCKET/$RUN/repo.tar.gz" --only-show-errors + sed -e "s/__BUCKET__/$BUCKET/" -e "s/__RUN__/$RUN/" -e "s/__MAX_MINUTES__/$MAX_MINUTES/" \ + 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 \ + --query Parameter.Value --output text) + echo "-> launching $TYPE spot (cap \$$MAX_PRICE/h, hard limit ${MAX_MINUTES} min, ${DISK_GB} GB disk deleted on exit)" + ID=$($AWS ec2 run-instances --image-id "$AMI" --instance-type "$TYPE" --count 1 \ + --iam-instance-profile Name=$ROLE \ + --instance-initiated-shutdown-behavior terminate \ + --instance-market-options "MarketType=spot,SpotOptions={MaxPrice=$MAX_PRICE,SpotInstanceType=one-time,InstanceInterruptionBehavior=terminate}" \ + --block-device-mappings "DeviceName=/dev/sda1,Ebs={VolumeSize=$DISK_GB,VolumeType=gp3,DeleteOnTermination=true}" \ + --metadata-options HttpTokens=required \ + --user-data file:///tmp/simurg-userdata.sh \ + --tag-specifications "ResourceType=instance,Tags=[{Key=Project,Value=$TAG},{Key=Run,Value=$RUN}]" \ + "ResourceType=volume,Tags=[{Key=Project,Value=$TAG},{Key=Run,Value=$RUN}]" \ + --query 'Instances[0].InstanceId' --output text) + echo "$ID $RUN" > .trial-last + echo "-> idle alarm: terminate if average CPU < 2% for 30 minutes" + $AWS cloudwatch put-metric-alarm --alarm-name "$TAG-idle-$ID" --namespace AWS/EC2 --metric-name CPUUtilization \ + --dimensions Name=InstanceId,Value="$ID" --statistic Average --period 300 --evaluation-periods 6 \ + --threshold 2 --comparison-operator LessThanThreshold --treat-missing-data notBreaching \ + --alarm-actions "arn:aws:automate:$REGION:ec2:terminate" + echo "launched $ID ($RUN). Check with: bash scripts/cloud/trial.sh status" + ;; +status) + read -r ID RUN < .trial-last + $AWS ec2 describe-instances --instance-ids "$ID" \ + --query 'Reservations[0].Instances[0].[InstanceId,State.Name,InstanceType,LaunchTime]' --output text + $AWS s3 cp "s3://$BUCKET/$RUN/log.txt" - 2>/dev/null | tail -n "${LINES_:-15}" || echo "(no log uploaded yet)" + ;; +results) + read -r ID RUN < .trial-last + mkdir -p "results/$RUN" + $AWS s3 cp "s3://$BUCKET/$RUN/" "results/$RUN/" --recursive --exclude repo.tar.gz --only-show-errors + ls -la "results/$RUN" + ;; +kill) + ids=$($AWS ec2 describe-instances --filters Name=tag:Project,Values=$TAG \ + Name=instance-state-name,Values=pending,running,stopping,stopped --query 'Reservations[].Instances[].InstanceId' --output text) + [[ -n "$ids" ]] && $AWS ec2 terminate-instances --instance-ids $ids --output text || echo "no trial instances" + ;& +audit) + echo "-- instances (any project) --" + $AWS ec2 describe-instances --filters Name=instance-state-name,Values=pending,running,stopping,stopped \ + --query 'Reservations[].Instances[].[InstanceId,State.Name,InstanceType]' --output text + echo "-- volumes --" + $AWS ec2 describe-volumes --query 'Volumes[].[VolumeId,State,Size]' --output text + echo "-- snapshots owned --" + $AWS ec2 describe-snapshots --owner-ids self --query 'Snapshots[].[SnapshotId,VolumeSize]' --output text + echo "-- elastic IPs --" + $AWS ec2 describe-addresses --query 'Addresses[].PublicIp' --output text + echo "-- trial alarms (deleting alarms of terminated instances) --" + for a in $($AWS cloudwatch describe-alarms --alarm-name-prefix "$TAG-idle-" --query 'MetricAlarms[].AlarmName' --output text); do + iid=${a#"$TAG-idle-"} + st=$($AWS ec2 describe-instances --instance-ids "$iid" --query 'Reservations[0].Instances[0].State.Name' --output text 2>/dev/null || echo gone) + if [[ "$st" == "terminated" || "$st" == "gone" ]]; then $AWS cloudwatch delete-alarms --alarm-names "$a"; echo "deleted $a"; else echo "$a ($st)"; fi + done + ;; +*) + sed -n '2,15p' "$0" + ;; +esac diff --git a/scripts/cloud/userdata.sh b/scripts/cloud/userdata.sh new file mode 100755 index 0000000..a700ded --- /dev/null +++ b/scripts/cloud/userdata.sh @@ -0,0 +1,62 @@ +#!/bin/bash +# EC2 user-data for a self-destructing SimurgArray validation run (Ubuntu 24.04, arm64). +# Placeholders (__X__) are filled in by scripts/cloud/trial.sh. +# +# Safety layers inside the instance: +# 1. `shutdown -h +MAX_MINUTES` is scheduled first -> hard time cap +# 2. an EXIT trap uploads the log and shuts down -> ends as soon as the job ends, success or failure +# Shutdown = TERMINATE (instance-initiated-shutdown-behavior) and the only disk has DeleteOnTermination=true, +# so nothing keeps costing money afterwards. A CloudWatch alarm (layer 3) terminates an idle instance from +# outside, in case the OS hangs. +set -uo pipefail +BUCKET="__BUCKET__" +RUN="__RUN__" +MAX_MINUTES="__MAX_MINUTES__" +shutdown -h +"$MAX_MINUTES" "SimurgArray trial hard time limit" + +LOG=/var/log/simurg-trial.log +exec > >(tee -a "$LOG") 2>&1 +S3="" +upload_log() { [[ -n "$S3" ]] && $S3 cp "$LOG" "s3://$BUCKET/$RUN/log.txt" --only-show-errors || true; } +finish() { + rc=$? + echo "=== finished rc=$rc at $(date -u +%FT%TZ) ===" + echo "$rc" > /tmp/rc && [[ -n "$S3" ]] && $S3 cp /tmp/rc "s3://$BUCKET/$RUN/exit_code.txt" --only-show-errors + upload_log + shutdown -h now +} +trap finish EXIT +set -e +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 +snap install aws-cli --classic +S3="/snap/bin/aws s3" +( while sleep 300; do upload_log; done ) & + +mkdir -p /work && cd /work +$S3 cp "s3://$BUCKET/$RUN/repo.tar.gz" . --only-show-errors +mkdir repo && tar -xzf repo.tar.gz -C repo && cd repo +curl -LsSf https://astral.sh/uv/install.sh | sh +export PATH="$HOME/.local/bin:$PATH" +uv venv -q --python 3.12 .venv && uv pip install -q -e ".[dev]" +mkdir -p refs && git clone -q --depth 1 https://github.com/Illumina/BeadArrayFiles.git refs/BeadArrayFiles + +SRC=s3://1000genomes/release/20130502/supporting/hd_genotype_chip/broad_intensities +D=/work/omni && mkdir -p $D/data +for f in HumanOmni2.5-4v1-Multi_B.bpm HumanOmni2.5-4v1-Multi_B.egt HumanOmni2.5-4v1_B.csv; do + $S3 cp --no-sign-request "$SRC/$f" $D/ --only-show-errors +done +echo "=== streaming 107 GB IDAT+GTC archive (no copy of the .tgz is kept) $(date -u +%T) ===" +$S3 cp --no-sign-request "$SRC/Omni25_idats_gtcs_2141_samples.tgz" - | pigz -dc | tar -x -C $D/data +echo "=== extracted $(du -sh $D/data | cut -f1) at $(date -u +%T) ===" + +echo "=== validation $(date -u +%T) ===" +.venv/bin/python -m pytest tests/unit -q -p no:cacheprovider +.venv/bin/python scripts/validate_p1_large.py --root $D/data --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm \ + --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 "=== done $(date -u +%T) ===" diff --git a/scripts/validate_p1_large.py b/scripts/validate_p1_large.py new file mode 100644 index 0000000..3e96f2a --- /dev/null +++ b/scripts/validate_p1_large.py @@ -0,0 +1,225 @@ +#!/usr/bin/env python3 +"""Phase P1 validation on a large data set (e.g. 1000 Genomes Omni2.5, 2141 samples with vendor GTCs). + +For every GTC: parse, compare every field with Illumina BeadArrayFiles (BSD-2) and test the byte-identical +round-trip. For every IDAT pair: parse, pair-check and compare the address set with the manifest. BPM/EGT/CSV are +compared with BeadArrayFiles once. Writes a JSON summary; exit code 0 only if everything matched. + +Usage: validate_p1_large.py --root DIR --bpm F --egt F [--csv F] --out summary.json [--workers N] [--limit N] +""" + +from __future__ import annotations + +import argparse +import json +import os +import sys +import time +import traceback +from concurrent.futures import ProcessPoolExecutor +from pathlib import Path +from typing import Any + +import numpy as np + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT / "refs" / "BeadArrayFiles")) + + +def check_gtc(path: str) -> dict[str, Any]: + from module.GenotypeCalls import GenotypeCalls # type: ignore[import-not-found] + + from simurg.formats.gtc import encode_gtc, read_gtc + + res: dict[str, Any] = {"file": path, "ok": True, "diffs": []} + try: + g = read_gtc(path, check_write_complete=False) + raw = Path(path).read_bytes() + res["version"] = g.version + res["autocall_version"] = g.autocall_version + res["roundtrip_identical"] = encode_gtc(g) == raw + if not res["roundtrip_identical"]: + res["roundtrip_canonical_identical"] = encode_gtc(g, canonical=True) == raw + res["size_ours_vs_file"] = [len(encode_gtc(g)), len(raw)] + b = GenotypeCalls(path, check_write_complete=False) + pairs = { + "raw_x": (g.raw_x, b.get_raw_x_intensities()), + "raw_y": (g.raw_y, b.get_raw_y_intensities()), + "genotypes": (g.genotypes, np.array(b.get_genotypes(), dtype=np.uint8)), + "scores": (g.genotype_scores, b.get_genotype_scores()), + } + if g.version >= 4: + pairs["baf"] = (g.b_allele_freqs, b.get_ballele_freqs()) + pairs["lrr"] = (g.logr_ratios, b.get_logr_ratios()) + for k, (a, c) in pairs.items(): + if not np.array_equal(np.asarray(a), np.asarray(c), equal_nan=a.dtype.kind == "f"): + res["diffs"].append(k) + scalars = { + "call_rate": (g.call_rate, b.get_call_rate()), + "gender": (g.gender, b.get_gender()), + "logr_dev": (g.logr_dev, b.get_logr_dev()), + "gc10": (g.gc10, b.get_gc10()), + "gc50": (g.gc50, b.get_gc50()), + "num_calls": (g.num_calls, b.get_num_calls()), + "sample_name": (g.sample_name, b.get_sample_name()), + "base_calls": (g.base_call_strings(), [x.decode() for x in b.get_base_calls()]), + } + for k, (a, c) in scalars.items(): + same = (a == c) or (isinstance(a, float) and np.isnan(a) and np.isnan(c)) + if not same: + res["diffs"].append(k) + xf_b = [ + (t.offset_x, t.offset_y, t.scale_x, t.scale_y, t.shear, t.theta) for t in b.get_normalization_transforms() + ] + xf_g = [ + tuple(float(x[k]) for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta")) + for x in g.normalization_transforms + ] + if [tuple(float(v) for v in t) for t in xf_b] != xf_g: + res["diffs"].append("normalization_transforms") + res["ok"] = not res["diffs"] + except Exception as e: # noqa: BLE001 - report everything + res.update(ok=False, error=f"{type(e).__name__}: {e}", trace=traceback.format_exc(limit=3)) + return res + + +def check_idat_pair(args: tuple[str, str]) -> dict[str, Any]: + from simurg.formats.idat import check_pair, read_idat + + grn, red = args + res: dict[str, Any] = {"grn": grn, "ok": True} + try: + g, r = read_idat(grn), read_idat(red) + check_pair(g, r) + res.update( + n=g.n_addresses, + chip=g.chip_type, + midblock_itemsize=g.mid_block_itemsize, + scan=g.scan_time_iso, + name_ok=(g.barcode in Path(grn).name), + ) + except Exception as e: # noqa: BLE001 + 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) + ap.add_argument("--bpm", required=True, type=Path) + ap.add_argument("--egt", required=True, type=Path) + ap.add_argument("--csv", type=Path) + 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) + a = ap.parse_args() + + from simurg.formats.bpm import read_bpm + from simurg.formats.discovery import walk_files + from simurg.formats.egt import read_egt + + t0 = time.time() + summary: dict[str, Any] = {"root": str(a.root), "started": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime())} + + gtcs = [str(p) for p in walk_files(a.root, (".gtc",))] + idats = walk_files(a.root, (".idat", ".idat.gz")) + if a.limit: + gtcs = gtcs[: a.limit] + grn = {p.name.lower().replace("_grn.idat", ""): p for p in idats if "_grn.idat" in p.name.lower()} + red = {p.name.lower().replace("_red.idat", ""): p for p in idats if "_red.idat" in p.name.lower()} + pairs = [(str(grn[k]), str(red[k])) for k in sorted(set(grn) & set(red))] + if a.limit: + pairs = pairs[: a.limit] + summary["n_gtc"] = len(gtcs) + summary["n_idat_pairs"] = len(pairs) + summary["unpaired_idats"] = len(set(grn) ^ set(red)) + + with ProcessPoolExecutor(a.workers) as ex: + gres = list(ex.map(check_gtc, gtcs, chunksize=4)) + ires = list(ex.map(check_idat_pair, pairs, chunksize=8)) + summary["gtc_ok"] = sum(r["ok"] for r in gres) + summary["gtc_roundtrip_identical"] = sum(bool(r.get("roundtrip_identical")) for r in gres) + summary["gtc_versions"] = sorted({r.get("version") for r in gres if "version" in r}) + summary["gtc_autocall_versions"] = sorted({r.get("autocall_version") for r in gres if "autocall_version" in r}) + summary["gtc_failures"] = [r for r in gres if not r["ok"]][:50] + summary["gtc_roundtrip_failures"] = [r for r in gres if r.get("roundtrip_identical") is False][:20] + summary["idat_ok"] = sum(r["ok"] for r in ires) + summary["idat_failures"] = [r for r in ires if not r["ok"]][:50] + summary["idat_midblock_itemsizes"] = sorted({r.get("midblock_itemsize") for r in ires if r["ok"]}) + summary["idat_chip_types"] = sorted({r.get("chip") for r in ires if r["ok"]}) + + # manifest / cluster file vs BeadArrayFiles + import module # type: ignore[import-not-found] + + m = read_bpm(a.bpm) + b = module.BeadPoolManifest(str(a.bpm)) + summary["bpm"] = { + "num_loci": m.num_loci, + "names": m.names == b.names, + "snps": m.snps == b.snps, + "chroms": m.chroms == b.chroms, + "map_infos": m.map_infos.tolist() == b.map_infos, + "address_a": m.address_a.tolist() == b.addresses, + "assay_types": m.assay_types.tolist() == b.assay_types, + "normalization_lookups": m.normalization_lookups().tolist() == b.normalization_lookups, + "n_bins": int(len(np.unique(m.norm_bin_ids()))), + } + e = read_egt(a.egt) + with open(a.egt, "rb") as fh: + be = module.ClusterFile.read_cluster_file(fh) + bad = 0 + + def same(a: float, b: float) -> bool: + return bool(np.float32(a) == np.float32(b) or (np.isnan(a) and np.isnan(b))) + + for i, name in enumerate(e.names): + rec = be.get_record(name) + st = (rec.aa_cluster_stats, rec.ab_cluster_stats, rec.bb_cluster_stats) + if not all( + same(e.theta_mean[i, k], s.theta_mean) + and same(e.r_mean[i, k], s.r_mean) + and same(e.theta_dev[i, k], s.theta_dev) + and same(e.r_dev[i, k], s.r_dev) + and e.n[i, k] == s.N + for k, s in enumerate(st) + ): + bad += 1 + summary["egt"] = { + "num_records": e.num_records, + "dbv": e.data_block_version, + "extra_floats": e.extra_per_record, + "record_mismatches": bad, + } + if a.csv: + from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest + + summary["bpm_csv_problems"] = compare_manifests(m, read_csv_manifest(a.csv, load_sequences=False)) + if gtcs: + g0 = __import__("simurg.formats.gtc", fromlist=["read_gtc"]).read_gtc(gtcs[0], check_write_complete=False) + summary["gtc_vs_bpm_num_snps"] = [g0.num_snps, m.num_loci] + summary["gtc_transforms_vs_bins"] = [len(g0.normalization_transforms), summary["bpm"]["n_bins"]] + + summary["seconds"] = round(time.time() - t0, 1) + passed = ( + 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)) + and summary["egt"]["record_mismatches"] == 0 + and not summary.get("bpm_csv_problems") + ) + summary["passed"] = passed + a.out.parent.mkdir(parents=True, exist_ok=True) + a.out.write_text(json.dumps(summary, indent=2, default=str) + "\n") + print( + json.dumps( + { + k: summary[k] + for k in ("n_gtc", "gtc_ok", "gtc_roundtrip_identical", "n_idat_pairs", "idat_ok", "passed", "seconds") + } + ) + ) + return 0 if passed else 1 + + +if __name__ == "__main__": + sys.exit(main())