"""Single entry point of the gate_lag study: regenerates every number and figure. python3 run_all.py (11 to 15 minutes on one core; raw inputs are cached after the first run) Steps 1. fetch_data : download raw inputs into fermi_raw/gate_lag (skipped when already present) 2. chain : self-test of the likelihood engine against closed forms 3. part1 : classic transitions, hard steps versus gating (Bayes factors, lags, flexibility, spacing) and the posterior over the number of hard steps (Part 3) 4. part2 : Phanerozoic innovations against oxygen, temperature and CO2 reconstructions 5. figures : five PNG files in ../figures 6. results.json in ../data with the headline numbers quoted in the report All random numbers come from numpy Generators seeded from config.SEED, so reruns reproduce the tables exactly. """ import json import os import sys import time import numpy as np import pandas as pd sys.dont_write_bytecode = True sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) # so that `python3 -I run_all.py` also works import chain import fetch_data import figures import part1 import part2 from config import DATA, SEED def _f(x, nd=3): return None if x is None or (isinstance(x, float) and np.isnan(x)) else round(float(x), nd) def main(): t0 = time.time() fetch_data.main() st = chain.selftest() print("likelihood engine self-test (max relative error):", {k: "%.1e" % v for k, v in st.items()}) assert max(st.values()) < 5e-3, "normaliser self-test failed" r1 = part1.run_all() print("part 1 done, %.0f s" % (time.time() - t0)) r2 = part2.run_all() print("part 2 done, %.0f s" % (time.time() - t0)) figs = figures.run_all() bf, post, flex, hl, sp = r1["bf"], r1["post"], r1["flex"], r1["headline"], r1["spacing"] agn = post[post["photosynthesis"] == "flat"] def bfrow(sub): r = bf[bf["config"].str.contains(sub, regex=False)].iloc[0] return dict(config=r["config"], set=r["set"], BF_G_H=_f(r["BF_G_H"]), BF_G_Hg=_f(r["BF_G_Hg"]), lnse_batch_G_H=_f(r["lnse_batch_G_H"]), lnse_batch_G_Hg=_f(r["lnse_batch_G_Hg"])) def spread(group): g = bf[bf["group"] == group] return dict(n=int(len(g)), BF_G_H=[_f(g["BF_G_H"].min()), _f(g["BF_G_H"].max())], BF_G_Hg=[_f(g["BF_G_Hg"].min()), _f(g["BF_G_Hg"].max())]) tests = r2["tests"] prim = tests[(tests["statistic"] == "rise_20Myr") & (tests["subset"] == "all") & tests["dates"].str.startswith("fossil")] sens = r2["sens"] gm = sens[(sens["family"] == "GEOCARBSULFOR MC") & (sens["statistic"] == "rise_20Myr")] res = dict( seed=SEED, part1=dict( headline=hl, BF_G_H_range_all_configs=[_f(bf["BF_G_H"].min()), _f(bf["BF_G_H"].max())], BF_G_Hg_range_all_configs=[_f(bf["BF_G_Hg"].min()), _f(bf["BF_G_Hg"].max())], n_configs=int(len(bf)), n_configs_G_H_above_3=int((bf["BF_G_H"] > 3).sum()), n_configs_G_H_above_10=int((bf["BF_G_H"] > 10).sum()), n_configs_G_H_below_1=int((bf["BF_G_H"] < 1).sum()), n_configs_G_Hg_above_3=int((bf["BF_G_Hg"] > 3).sum()), n_configs_G_Hg_below_1=int((bf["BF_G_Hg"] < 1).sum()), spread_by_group={g: spread(g) for g in bf["group"].unique()}, step_set=[bfrow(c) for c in ("conventional gates", "three steps counted from oxygenic")] + [dict(config=r["config"], set=r["set"], BF_G_H=_f(r["BF_G_H"]), BF_G_Hg=_f(r["BF_G_Hg"])) for _, r in bf[(bf["config"] == "conventional gates") & (bf["set"] != "C3")].iterrows()], by_LECA_study=[bfrow("LECA " + c) for c in ("Betts2018", "Kay2026", "Parfrey2011", "Strassert2021")], by_Metazoa_study=[bfrow("Metazoa " + c) for c in ("LoleDurbin2026", "dosReis2015", "Carlisle2024")], by_animal_gate=[dict(config=r["config"], BF_G_H=_f(r["BF_G_H"], 4), BF_G_Hg=_f(r["BF_G_Hg"], 4)) for _, r in bf[bf["group"] == "animal gate"].iterrows()], timing_study_reconciliation=[bfrow("timing-study inputs (Kay"), bfrow("timing-study inputs and its")], late_window_start=dict(**bfrow("start of window 4.25"), frac_luca_redrawn=_f(bf[bf["config"].str.contains("start of window")] ["frac_luca_redrawn"].iloc[0])), flexibility=flex, spacing=sp.to_dict(orient="records"), ), part3=dict( candidate_step_prior_k=[_f(x, 4) for x in part1.MODEL_PRIOR_K], implied_hard_step_prior_k=[_f(x, 4) for x in part1.IMPLIED_PRIOR_K], implied_prior_thinning_q=_f(part1.Q_HARD, 4), posterior=post[["config", "photosynthesis"] + ["post_k%d" % k for k in range(5)] + ["post_mean_k", "BF_k0_vs_kge1", "BF_k0_vs_kge1_lnse"] + ["imp_post_k%d" % k for k in range(5)] + ["imp_post_mean_k", "imp_BF_k0_vs_kge1"] + [c for c in post.columns if c.startswith("P_hard_")]].round(4).to_dict(orient="records"), n_variants_flat=int(len(agn)), post_k0_range_flat=[_f(agn["post_k0"].min()), _f(agn["post_k0"].max())], post_k0_range_all=[_f(post["post_k0"].min()), _f(post["post_k0"].max())], imp_post_k0_range_flat=[_f(agn["imp_post_k0"].min()), _f(agn["imp_post_k0"].max())], imp_post_k0_range_all=[_f(post["imp_post_k0"].min()), _f(post["imp_post_k0"].max())], BF_k0_range_flat=[_f(agn["BF_k0_vs_kge1"].min()), _f(agn["BF_k0_vs_kge1"].max())], BF_k0_range_all=[_f(post["BF_k0_vs_kge1"].min()), _f(post["BF_k0_vs_kge1"].max())], imp_BF_k0_range_flat=[_f(agn["imp_BF_k0_vs_kge1"].min()), _f(agn["imp_BF_k0_vs_kge1"].max())], post_mean_k_range_flat=[_f(agn["post_mean_k"].min()), _f(agn["post_mean_k"].max())], imp_post_mean_k_range_flat=[_f(agn["imp_post_mean_k"].min()), _f(agn["imp_post_mean_k"].max())], ), part2=dict( n_innovations=int(len(r2["tab"])), n_pbdb_dated=int((r2["tab"]["fad_basis"] == "PBDB").sum()), primary=prim[["curve", "n_innovations", "observed", "null_shift_mean", "null_shift_p05", "null_shift_p95", "p_shift_ge", "p_jitter_ge", "p_pbdb_ge", "z_shift"]] .round(3).to_dict(orient="records"), geocarbsulfor_mc=dict(n=int(len(gm)), z_median=_f(gm["z"].median()), z_min=_f(gm["z"].min()), z_max=_f(gm["z"].max()), frac_p_below_05=_f((gm["p_ge"] < 0.05).mean())), summary=r2["summary"], window=r2["window"].round(3).to_dict(orient="records"), power=r2["power"].round(3).to_dict(orient="records"), power_vs_n=r2["power_n"].round(3).to_dict(orient="records"), multiplicity=r2["multiplicity"], literature_dates_all=tests[(tests["statistic"] == "rise_20Myr") & (tests["subset"] == "all") & tests["dates"].str.startswith("literature")] [["curve", "n_innovations", "observed", "null_shift_mean", "p_shift_ge"]].round(3).to_dict(orient="records"), pbdb_younger_than_literature=r2["tab"][r2["tab"]["pbdb_younger_than_literature"]] [["id", "pbdb_fad_taxon", "pbdb_fad_old", "pbdb_fad_young", "lit_old", "lit_young"]].to_dict(orient="records"), curve_agreement=r2["agreement"].round(3).to_dict(orient="records"), ), figures=[os.path.basename(f) for f in figs], runtime_seconds=round(time.time() - t0, 1), ) with open(os.path.join(DATA, "results.json"), "w") as f: json.dump(res, f, indent=1, default=str) pd.set_option("display.width", 220); pd.set_option("display.max_colwidth", 80) print("\n=== Part 1: reference configuration (three gated transitions, conventional oxygenation-event gates, " "equal-weight mixture of clock studies, lag prior log-uniform 10-500 Myr) ===") print(json.dumps({k: v for k, v in hl.items() if k != "batch_values_G_H"}, indent=1)) print("\n=== Part 1: Bayes factors, all configurations ===") print(bf[["config", "set", "BF_G_H", "lnse_G_H", "lnse_batch_G_H", "BF_G_Hg", "lnse_G_Hg", "lnse_batch_G_Hg", "BF_Hg_H", "frac_draws_G_zero", "frac_luca_redrawn"]].round(3).to_string()) print("\n=== Part 1: lags ===") print(r1["lags"].round(3).to_string()) print("\n=== Part 1: flexibility ===") print(json.dumps(flex, indent=1)) print("\n=== Part 1: spacing test ===") print(sp.round(3).to_string()) print("\n=== Part 3: posterior over the number of hard steps ===") print("candidate-step prior", np.round(part1.MODEL_PRIOR_K, 3), " implied hard-step prior", np.round(part1.IMPLIED_PRIOR_K, 3)) print(post[["config"] + ["post_k%d" % k for k in range(5)] + ["post_mean_k", "BF_k0_vs_kge1"] + ["imp_post_k0", "imp_post_mean_k", "imp_BF_k0_vs_kge1"] + [c for c in post.columns if c.startswith("P_hard_")]].round(3).to_string()) print("\n=== Part 2: primary test (mean O2 change over the 20 Myr before origin; circular-shift null) ===") print(prim[["curve", "n_innovations", "observed", "null_shift_mean", "null_shift_p05", "null_shift_p95", "p_shift_ge", "p_jitter_ge", "p_pbdb_ge", "z_shift"]].round(3).to_string()) print("\n=== Part 2: all-literature dates, window sensitivity, multiplicity ===") print(tests[(tests["statistic"] == "rise_20Myr") & (tests["subset"] == "all")][ ["curve", "dates", "n_innovations", "observed", "null_shift_mean", "p_shift_ge"]].round(3).to_string()) print(r2["window"].round(3).to_string()) print(json.dumps({k: v for k, v in r2["multiplicity"].items() if not k.startswith("significant")}, indent=1)) print(json.dumps(r2["summary"], indent=1)) print("\n=== Part 2: power ===") print(r2["power"].round(2).to_string()) print(r2["power_n"].round(2).to_string()) print("\nfigures:", figs) print("total seconds: %.0f" % (time.time() - t0)) if __name__ == "__main__": main()