"""Recompute every number of F0005 from data/ and assert it matches the recorded values (finding.yaml in the findings repo, expected_numbers.json in a downloaded copy). python check.py # assert python check.py --print # show computed values """ import csv import json import statistics as st import sys from pathlib import Path HERE = Path(__file__).parent def pct(k, n): return round(100 * k / n, 1) if n else 0.0 def table(name): return list(csv.DictReader(open(HERE / "data" / name, newline=""))) def compute(): screen = table("screen.csv") rows = table("search_outcomes.csv") kept = [r for r in screen if r["kept"] == "1"] left = [r for r in screen if r["status"] == "match" and r["kept"] == "0"] done = [r for r in rows if r["status"] == "done"] lower = [r for r in done if r["outcome"] == "lower"] kind = lambda k: [float(r["lower_by_kcal"]) for r in lower if r["best_kind"] == k] gap = [float(r["lower_by_kcal"]) for r in lower] med = lambda xs: round(st.median(xs), 2) if xs else None atoms = [int(r["n_atoms"]) for r in kept] status = lambda s: sum(r["status"] == s for r in screen) drop = lambda r: max(float(r["start_end_drop_kcal"]), float(r["end_end_drop_kcal"])) dropped = {r["reaction_id"] for r in rows if r["status"] == "no_route" and r["routes_dropped"] not in ("", "0")} # why no route: the reactions where no route was found at all (not those whose routes were all dropped) nr = [r for r in table("search_diagnostics.csv") if r["status"] == "no_route" and r["reaction_id"] not in dropped] q = lambda r: r["irc_pair_queued"] == "1" dn = lambda r: r["irc_pair_path_done"] == "1" sp = lambda r: r["irc_pair_path_split"] == "1" n_nr = sum(r["status"] == "no_route" for r in rows) return { # the screen: which Transition1x reactions make the benchmark set "screened": len(screen), "set": len(kept), "set_pct": pct(len(kept), len(screen)), "left_out_irc": len(left), "screen_one_end": status("start_only") + status("end_only"), "screen_mismatch": status("mismatch"), "screen_degenerate_irc": status("degenerate_irc"), "screen_end_not_minimized": status("irc_end_not_minimized") + status("irc_failed"), "screen_not_first_order": status("ts_wrong_order"), "screen_not_converged": status("ts_not_converged"), "screen_stereo_undetermined": status("stereo_undetermined"), "screen_error": status("error"), "left_out_drop_over_20_kcal": sum(drop(r) > 20 for r in left), "left_out_soft_saddle": sum(float(r["ts_imaginary_cm"]) > -200 for r in left), "set_soft_saddle": sum(float(r["ts_imaginary_cm"]) > -200 for r in kept), "set_median_end_drop_kcal": med([drop(r) for r in kept]), "set_median_barrier_kcal": med([float(r["barrier_kcal"]) for r in kept]), "n_atoms_min": min(atoms), "n_atoms_max": max(atoms), # the IRC retry (fmax 0.01) where a first IRC end changed species on minimization "irc_retried": sum(r["irc_retried"] != "0" for r in screen), "irc_retry_failed": sum(r["irc_retried"] == "failed" for r in screen), "irc_retry_kept": sum(r["irc_retried"] == "1" and r["kept"] == "1" for r in screen), "irc_retry_kept_was_match": sum(r["irc_retried"] == "1" and r["kept"] == "1" and r["loose_irc_status"] == "match" for r in screen), # the search "searched": len(rows), "with_route": len(done), "lower": len(lower), "lower_pct_of_searched": pct(len(lower), len(rows)), "lower_single_step": len(kind("direct")), "lower_multistep": len(kind("multistep")), "matched": sum(r["outcome"] == "matched" for r in done), "higher": sum(r["outcome"] == "higher" for r in done), "no_route": n_nr, "no_route_pct_of_searched": pct(n_nr, len(rows)), "time_limit": sum(r["status"] == "timeout" for r in rows), "crashed": sum(r["status"] not in ("done", "no_route", "timeout") for r in rows), "lower_median_kcal": med(gap), "lower_median_single_step_kcal": med(kind("direct")), "lower_median_multistep_kcal": med(kind("multistep")), "lower_max_kcal": round(max(gap), 2) if gap else None, "search_median_minutes": round(st.median(float(r["seconds"]) for r in rows if r["seconds"]) / 60), "top_n": 10, # route steps: a step counts only if its IRC ends keep their species (retried at fmax 0.01 once) "rescore_reactions_step_retried": sum(r["step_ircs_retried"] not in ("", "0") for r in rows), "rescore_step_ircs_retried": sum(int(r["step_ircs_retried"] or 0) for r in rows), "rescore_routes_dropped_reactions": sum(r["routes_dropped"] not in ("", "0") for r in rows), "rescore_no_route_all_dropped": sum(r["status"] == "no_route" and r["routes_dropped"] not in ("", "0") for r in rows), # why no route "no_route_none_found": len(nr), "no_route_irc_pair_not_queued": sum(not q(r) for r in nr), "no_route_irc_ends_not_in_pools": sum(r["irc_ends_in_pools"] != "1" for r in nr), "no_route_irc_pair_not_reached": sum(q(r) and not dn(r) for r in nr), "no_route_irc_pair_split": sum(dn(r) and sp(r) for r in nr), "no_route_irc_pair_single_step": sum(dn(r) and not sp(r) for r in nr), "no_route_budget_used_up": sum(r["budget_reached"] == "1" for r in nr), "no_route_found_other_ts": sum(r["n_offtarget"] not in ("", "0") for r in nr), } def expected(): if (HERE / "finding.yaml").exists(): import yaml return yaml.safe_load((HERE / "finding.yaml").read_text())["evidence"]["numbers"] return json.loads((HERE / "expected_numbers.json").read_text()) if __name__ == "__main__": got = compute() if "--print" in sys.argv: print(json.dumps(got, indent=1)) sys.exit(0) want = expected() bad = [f"{k}: recorded {want.get(k)!r}, data gives {v!r}" for k, v in got.items() if want.get(k) != v] bad += [f"{k}: recorded but not computed" for k in want if k not in got] if bad: print("\n".join(bad)) sys.exit(1) print(f"all {len(got)} numbers reproduce from data/")