"""Score the native-proposal replay against verified union-retrieval results.""" import argparse import collections import hashlib import itertools import json import random import sys from pathlib import Path B = Path(__file__).resolve().parent E = B.parent CV = E / "pangenome-reference-cv-20261014" sys.path.insert(0, str(CV / "asian-seven-comparison")) from endpoints import canonical_slots, named_pair, unmatched_cores def load(path): return json.loads(path.read_text()) def sha(path): return hashlib.sha256(path.read_bytes()).hexdigest() def uncertainty(pairs, group_of): clusters = collections.defaultdict(list) ungrouped = sorted({r["donor"] for r in pairs if r["donor"] not in group_of}) for r in pairs: clusters[group_of.get(r["donor"], (r["donor"],))].append(r["difference"]) values = [(sum(v), len(v)) for v in clusters.values()] if not values: return {"family_clusters": 0, "difference_pp": None, "bootstrap95_pp": None, "sign_flip_p": None, "singleton_fallback_donors": ungrouped} rng = random.Random(20261020) boots = [] for _ in range(4000): draw = [rng.choice(values) for _ in values] boots.append(100 * sum(x[0] for x in draw) / sum(x[1] for x in draw)) boots.sort() differences = [v[0] for v in values if v[0]] observed = abs(sum(differences)) if len(differences) <= 18: outcomes = [sum(s * d for s, d in zip(signs, differences)) for signs in itertools.product((-1, 1), repeat=len(differences))] pvalue = sum(abs(x) >= observed for x in outcomes) / len(outcomes) pmethod = "exact_family_sign_flip" else: hits = sum(abs(sum(rng.choice((-1, 1)) * d for d in differences)) >= observed for _ in range(20000)) pvalue = (hits + 1) / 20001 pmethod = "monte_carlo_family_sign_flip_20000" return {"family_clusters": len(values), "singleton_fallback_donors": ungrouped, "difference_pp": 100 * sum(x[0] for x in values) / sum(x[1] for x in values), "bootstrap95_pp": [boots[100], boots[3899]], "sign_flip_p": pvalue, "sign_flip_method": pmethod} def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--runs", type=Path, nargs="+", default=[E / "native-proposal-ablation-20261020/collected/runs"]) parser.add_argument("--graph-prefix", choices=["DogoHLA", "DogoHLA_genomic"], default="DogoHLA") parser.add_argument("--nested-native", action="store_true") parser.add_argument("--output-name", default="GRAPH_RETRIEVAL_ABLATION.json") a = parser.parse_args() original = load(B / "GRAPH_RESULTS.json") index = {(r["donor"], r["gene"], r["method"]): r for r in original["rows"]} folds = load(CV / "FOLDS.json") group_of = {d: tuple(sorted(g)) for g in folds["groups"] for d in g} truth = collections.defaultdict(list) for r in load(CV / "SOURCE_RECORDS.json")["records"]: truth[r["donor"], r["gene"]].append(r) rows, receipts, fragments = [], {}, [] pattern = "*/*/native_proposals_only/COMPLETE.json" if a.nested_native else "*/*/COMPLETE.json" seen = set() for receipt_path in sorted(p for root in a.runs for p in root.glob(pattern)): receipt = load(receipt_path) assert receipt["status"] == "native_proposal_ablation_complete" assert receipt["source_complete_sha256"] in original["sources_sha256"].values() arm_folder = receipt_path.parent.parent if a.nested_native else receipt_path.parent donor, panel = arm_folder.parent.name, arm_folder.name assert (donor, panel) not in seen, "Duplicate donor-panel replay" seen.add((donor, panel)) assert Path(receipt["source"]).parts[-2:] == (donor, panel) receipts[str(receipt_path)] = sha(receipt_path) fragments.append({"donor": donor, "panel": panel, **receipt["retrieval_counts"]}) for mode in ("E1", "E2", "E3"): folder = receipt_path.parent / mode payload = {} for name in ("calls.json", "sequence_calls.json", "decisions.json"): assert sha(folder / name) == receipt["output_sha256"][mode + "/" + name] payload[name] = load(folder / name) for gene in ("A", "B", "C", "DPA1", "DPB1", "DQA1", "DQB1", "DRB1"): union = index[donor, gene, f"{a.graph_prefix}/{panel}/{mode}"] slots = canonical_slots(payload["calls.json"].get(gene, {})) if len(slots) != 2: slots = [{"labels": [], "label_unknown": True}] * 2 for slot, core in zip(slots, payload["sequence_calls.json"][gene]["slots"]): slot.update(core) named = named_pair(slots, truth[donor, gene]) assert named["eligible"] == union["four_field"]["eligible"] rows.append({ "donor": donor, "panel": panel, "mode": mode, "gene": gene, "eligible": named["eligible"], "native_only": named, "union": union["four_field"], "native_only_unmatched": unmatched_cores(slots, truth[donor, gene]), "union_unmatched": union["unmatched"], "difference": int(union["four_field"]["correct"]) - int(named["correct"]), "native_only_reason": payload["decisions.json"]["decisions"].get(gene, {}).get("reason"), "union_reason": union["reason"], }) contrasts = [] for panel in ("hprc", "hprc_asian"): for mode in ("E1", "E2", "E3"): subset = [r for r in rows if r["panel"] == panel and r["mode"] == mode] eligible = [r for r in subset if r["eligible"]] contrasts.append({ "panel": panel, "mode": mode, "donors": len({r["donor"] for r in subset}), "eligible": len(eligible), "union_correct": sum(r["union"]["correct"] for r in eligible), "native_only_correct": sum(r["native_only"]["correct"] for r in eligible), "union_only_correct": sum(r["difference"] == 1 for r in eligible), "native_only_correct_union_wrong": sum(r["difference"] == -1 for r in eligible), "unnamed": sum(r["union_unmatched"]["eligible_copies"] for r in subset), "union_unnamed_recovered": sum(r["union_unmatched"]["exact_copies_recovered"] for r in subset), "native_only_unnamed_recovered": sum(r["native_only_unmatched"]["exact_copies_recovered"] for r in subset), **uncertainty(eligible, group_of), }) out = { "status": "scored_available_replays", "graph_prefix": a.graph_prefix, "completed_donor_panel_arms": len(receipts), "contrasts": contrasts, "rows": rows, "retrieval_counts": fragments, "receipt_sha256": receipts, "graph_results_sha256": sha(B / "GRAPH_RESULTS.json"), "folds_sha256": sha(CV / "FOLDS.json"), "script_sha256": sha(Path(__file__)), "limitations": [ "Graph contribution means additional graph-proposed placements within this fixed hybrid pipeline.", ("Full-genomic/flanked catalogue." if a.graph_prefix == "DogoHLA_genomic" else "Original.fa/flanked catalogue; no claim of superiority to full-genomic IPD."), "Family bootstrap/sign-flip conditional on frozen trained folds; overlapping training sets are not independently resampled.", "Six descriptive mode/panel contrasts; no multiplicity-adjusted confirmatory claim.", "Partial completion must not be reported as the planned full cohort.", ], } (B / a.output_name).write_text(json.dumps(out, indent=2) + "\n") print(json.dumps(contrasts, indent=2)) if __name__ == "__main__": main()