Add self-destructing AWS trial and large-scale P1 validation script

scripts/cloud/trial.sh launches one spot instance that validates the
parsers on the 1000 Genomes Omni2.5 set (2141 samples, vendor GTCs)
streamed from the public S3 bucket, uploads the summary and terminates
itself. Cost guards: spot price cap, hard time limit (shutdown =
terminate), EXIT-trap shutdown, CloudWatch idle alarm that terminates
from outside, disk deleted with the instance, monthly budget alert, and
an `audit` command listing anything still billable.

scripts/validate_p1_large.py compares every GTC with BeadArrayFiles,
checks byte-identical round-trips, IDAT pairs and BPM/EGT/CSV.
This commit is contained in:
kemal committed 2026-09-27 19:30:27 +01:00
1 parent 5a20e7131c
commit 3bfb7b5757
4 files changed
+412

No files matched your search

+2
View File
@@ -16,3 +16,5 @@ __pycache__/
.ruff_cache/ .ruff_cache/
# rust # rust
target/ target/
.trial-last
results/
+123
View File
@@ -0,0 +1,123 @@
#!/usr/bin/env bash
# Self-destructing, budget-guarded AWS validation run for SimurgArray.
#
# bash scripts/cloud/trial.sh setup <alert-email> # 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 <alert-email>}
echo "-> budget alert: \$$BUDGET_USD/month, e-mail $EMAIL at 50%, 100% actual and 100% forecast"
cat > /tmp/simurg-budget.json <<EOF
{"BudgetName":"$TAG-monthly","BudgetLimit":{"Amount":"$BUDGET_USD","Unit":"USD"},"TimeUnit":"MONTHLY","BudgetType":"COST"}
EOF
notif() { echo "{\"Notification\":{\"NotificationType\":\"$1\",\"ComparisonOperator\":\"GREATER_THAN\",\"Threshold\":$2,\"ThresholdType\":\"PERCENTAGE\"},\"Subscribers\":[{\"SubscriptionType\":\"EMAIL\",\"Address\":\"$EMAIL\"}]}"; }
echo "[$(notif ACTUAL 50),$(notif ACTUAL 100),$(notif FORECASTED 100)]" > /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
+62
View File
@@ -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) ==="
+225
View File
@@ -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())