"""Part 1: hard steps versus environmental gating on the classic dated transitions. Models (see chain.py): H all steps hard and ungated; Hg all steps hard but each waits for its gate; G every gated step is easy (gate + exponential lag, mean lag tau_i drawn from a prior); M_S the subset S of the candidate steps is hard (gated-hard by default) and the rest are easy. A marginal likelihood is the average, over Monte Carlo draws of the dates (published posteriors and ranges), gate dates, window length and lag means, of the chain density. Bayes factors are ratios of those averages. The gating likelihood is heavy-tailed (most of it comes from the few draws in which a transition falls just after its gate draw), so the honest Monte Carlo error is the batch-to-batch scatter, which is reported next to the bootstrap error. Gate definitions. 'conventional' = the three oxygenation events of the textbook narrative (Great Oxidation 2.22-2.46 Ga for eukaryotes, a Neoproterozoic window 0.539-0.80 Ga for animals, Palaeozoic oxygenation 0.40-0.42 Ga for humans). The first pass called these 'the gates named by proponents'; that was wrong for animals. Mills et al. 2025 call the animal oxygen threshold 'controversial and essentially unknown' and commit only to eukaryotes and animals being excluded before ~2.2 Ga and humans before ~0.40 Ga. The animal gate is therefore varied explicitly (`animal_gate`), including the Mills et al. reading (the Great Oxidation). Steps with no independently dated gate (oxygenic photosynthesis; land plants) are treated in two ways: 'flat' : under gating the step's time is flat over the whole window (an undated gate with a flat prior). This is NOT neutral between the models. Under Carter's model an extra early step is an ordered uniform draw, which raises the density of the observed sequence by n/L = 4/5.2 per Gyr, whereas the flat gating step contributes only 1/L; adding photosynthesis therefore lowers gating : hard steps by a factor of about 3.5 (8 -> 2.4) and gating : gated-hard from 3 to 2. Inside the subset posterior of Part 3 the step has the same kernel whether hard or easy, so there it carries no information about its own status. (The first pass called this treatment 'agnostic'.) 'strict' : the step must follow its predecessor after an exponential lag with the same prior on the mean. A third route conditions both models on the photosynthesis date (`t0_from='oxyphoto'`): the chain of three gated steps starts at the dated origin of oxygenic photosynthesis instead of at LUCA. """ import itertools import json import os import numpy as np import pandas as pd import chain import classic from config import DATA, SEED from curves import load_gcso from math import comb # The project's model puts 0.30/0.30/0.22/0.12/0.06 on 0..4 CANDIDATE steps. A candidate is hard only if its # success probability in the window is below 0.3, so the model's implied prior on the number of HARD steps is # different: its prior table gives P(no hard steps) = 0.349 and mean 1.15. We reproduce that by binomial # thinning of the candidate prior with q = 1.15 / 1.34 (which gives P(0) = 0.347, mean 1.15). MODEL_PRIOR_K = np.array([0.30, 0.30, 0.22, 0.12, 0.06]) Q_HARD = 1.15 / float((MODEL_PRIOR_K * np.arange(5)).sum()) IMPLIED_PRIOR_K = np.array([sum(MODEL_PRIOR_K[n] * comb(n, k) * Q_HARD ** k * (1 - Q_HARD) ** (n - k) for n in range(k, 5)) for k in range(5)]) TAU_PRIORS = { "logU_10_500": ("loguniform", 0.010, 0.500), "logU_1_1000": ("loguniform", 0.001, 1.000), "logU_10_2000": ("loguniform", 0.010, 2.000), "logU_10_3000": ("loguniform", 0.010, 3.000), # the timing study's lag prior (reconciliation row only) "logU_10_5000": ("loguniform", 0.010, 5.000), # 'easy' stretched to mean lags as long as the window "U_0_500": ("uniform", 0.0, 0.500), } SETS = { "C3": ["eukaryotes", "animals", "humans"], "C4": ["oxyphoto", "eukaryotes", "animals", "humans"], "C5": ["life", "oxyphoto", "eukaryotes", "animals", "humans"], "C6": ["life", "oxyphoto", "eukaryotes", "animals", "landplants", "humans"], } UNGATED_STEPS = {"oxyphoto", "landplants"} # no independently dated gate # the animal gate is not agreed: an equal-weight mixture of three published readings (weights are a judgement) ANIMAL_GATE_MIX = {"mix": ["noe_window", "tonian_onset", "goe_permanent"]} def draw_tau(rng, prior, shape): kind, lo, hi = TAU_PRIORS[prior] if kind == "loguniform": return np.exp(rng.uniform(np.log(lo), np.log(hi), size=shape)) return rng.uniform(max(lo, 1e-4), hi, size=shape) def make_draws(rng, D, steps, gate_version="conventional", tau_prior="logU_10_500", future=(1.0,), euk_def="eukaryotes", fix=None, physio_kw=None, gcso=None, start=(4.404, 4.51), animal_gate="noe_window", t0_from="life"): """Return forward-time arrays for the chain model. fix : dict transition -> estimate label, or -> (lo, hi) for a uniform range in Ga. animal_gate : ENV key or (lo, hi) in Ga for the animal gate of the 'conventional' version. t0_from : 'life' (chain starts at LUCA, or at the start of the window when life is itself a step) or 'oxyphoto' (chain starts at the dated origin of oxygenic photosynthesis).""" fix = fix or {} A0 = rng.uniform(start[0], start[1], size=D) F = rng.choice(np.asarray(future, float), size=D) T = A0 + F ages, labels = {}, {} for s in steps: key = euk_def if s == "eukaryotes" else s ages[s], labels[s] = classic.draw_transition(rng, key, D, which=fix.get(s)) luca, frac_luca_redrawn = classic.draw_luca_before(rng, A0) if "life" in steps: ages["life"] = luca t0 = np.zeros(D) elif t0_from == "oxyphoto": oxy, _ = classic.draw_transition(rng, "oxyphoto", D) t0 = A0 - oxy else: t0 = A0 - luca gate_age = {s: np.full(D, np.inf) for s in steps} # inf age = gate open from the start aux = {} if gate_version == "conventional": if "eukaryotes" in steps: gate_age["eukaryotes"] = rng.uniform(2.22, 2.46, size=D) if "animals" in steps: gate_age["animals"] = classic.draw_env(rng, animal_gate, D) if "humans" in steps: gate_age["humans"] = classic.draw_env(rng, "poe", D) elif gate_version == "physio": pg = classic.physio_gates(rng, D, gcso, **(physio_kw or {})) for s in ("eukaryotes", "animals", "humans"): if s in steps: gate_age[s] = pg[s] aux = pg elif gate_version == "none": pass else: raise ValueError(gate_version) if "life" in steps: gate_age["life"] = A0.copy() # habitable from the start of the window t = np.column_stack([A0 - ages[s] for s in steps]) g = np.column_stack([A0 - gate_age[s] for s in steps]) # -inf where no gate g = np.where(np.isfinite(g), g, -1.0) tau = draw_tau(rng, tau_prior, (D, len(steps))) return dict(steps=list(steps), t=t, g=g, lam=1.0 / tau, t0=t0, T=T, A0=A0, ages=ages, gate_age=gate_age, labels=labels, aux=aux, prior=TAU_PRIORS[tau_prior], frac_luca_redrawn=frac_luca_redrawn) def kinds_for(steps, hard_set, hard_kind="G", ungated="flat"): out = [] for s in steps: if s in UNGATED_STEPS: if ungated == "flat": out.append("H") # flat kernel whether hard or gated by an undated gate else: out.append("H" if s in hard_set else "E") # strict: easy means prompt after the predecessor else: if s in hard_set: out.append(hard_kind) else: out.append("E") return out def lik(dr, kinds): return chain.likelihood(kinds, dr["t"], dr["g"], dr["lam"], dr["t0"], dr["T"], prior=dr["prior"]) def boot_ratio(a, b, rng, nboot=400): """Ratio of means with a bootstrap standard error of its natural log.""" n = len(a) r = a.mean() / b.mean() if b.mean() > 0 else np.nan idx = rng.integers(0, n, size=(nboot, n)) am = a[idx].mean(axis=1) bm = b[idx].mean(axis=1) ok = (am > 0) & (bm > 0) se = float(np.std(np.log(am[ok] / bm[ok]))) if ok.sum() > 10 else np.nan return float(r), se def pure_models(dr, ungated="flat"): steps = dr["steps"] L_H = lik(dr, ["H"] * len(steps)) L_Hg = lik(dr, kinds_for(steps, set(steps), "G", ungated)) L_G = lik(dr, kinds_for(steps, set(), "G", ungated)) return L_H, L_Hg, L_G def bayes_factor_table(gcso, D=12000, seed=SEED, batches=3): """Each configuration is evaluated on `batches` independent batches of D draws, pooled. Two Monte Carlo errors are given for ln BF: a bootstrap over the pooled draws (lnse_*) and the scatter between batches divided by sqrt(batches) (lnse_batch_*). The second is the honest one for this heavy-tailed integrand.""" rng = np.random.default_rng(seed) rows = [] def add(label, group, **kw): ung = kw.pop("ungated", "flat") set_name = kw.pop("set_name") nb = kw.pop("batches", batches) parts, redrawn = [], [] for b in range(nb): dr = make_draws(rng, D, SETS[set_name], gcso=gcso, **kw) parts.append(pure_models(dr, ung)) redrawn.append(dr["frac_luca_redrawn"]) L_H, L_Hg, L_G = (np.concatenate([p[i] for p in parts]) for i in range(3)) r1, s1 = boot_ratio(L_G, L_H, rng) r2, s2 = boot_ratio(L_G, L_Hg, rng) r3, s3 = boot_ratio(L_Hg, L_H, rng) def batch_se(i, j): v = np.array([p[i].mean() / p[j].mean() if p[j].mean() > 0 else np.nan for p in parts]) v = v[np.isfinite(v) & (v > 0)] return float(np.std(np.log(v), ddof=1) / np.sqrt(len(v))) if len(v) > 1 else np.nan rows.append(dict(config=label, group=group, set=set_name, gate_version=kw.get("gate_version", "conventional"), animal_gate=str(kw.get("animal_gate", "noe_window")), t0_from=kw.get("t0_from", "life"), tau_prior=kw.get("tau_prior", "logU_10_500"), ungated_steps=ung, euk_def=kw.get("euk_def", "eukaryotes"), future=str(kw.get("future", (1.0,))), start=str(kw.get("start", (4.404, 4.51))), fix=json.dumps(kw.get("fix") or {}), physio=json.dumps(kw.get("physio_kw") or {}), BF_G_H=r1, lnse_G_H=s1, lnse_batch_G_H=batch_se(2, 0), BF_G_Hg=r2, lnse_G_Hg=s2, lnse_batch_G_Hg=batch_se(2, 1), BF_Hg_H=r3, lnse_Hg_H=s3, frac_draws_G_zero=float(np.mean(L_G == 0)), frac_luca_redrawn=float(np.mean(redrawn)), n_draws=D * nb)) cg = "conventional gates" # --- which steps are counted -------------------------------------------------------------------------- for sname in ("C3", "C4", "C5"): add(cg, "step set", set_name=sname) add(cg + ", three steps counted from oxygenic photosynthesis (both models conditioned on its date)", "step set", set_name="C3", t0_from="oxyphoto") add(cg + ", C4 strict (photosynthesis prompt after life)", "step set", set_name="C4", ungated="strict") add(cg + ", C5 strict", "step set", set_name="C5", ungated="strict") add(cg + ", C5, start of window 4.25-4.40 Ga (after last ocean-vaporising impact)", "step set", set_name="C5", start=(4.25, 4.40)) # --- which clock study -------------------------------------------------------------------------------- for lab in ("Parfrey2011", "Betts2018", "Strassert2021", "Kay2026"): add(cg + ", LECA " + lab, "clock study", set_name="C3", fix={"eukaryotes": lab}) for lab in ("dosReis2015", "LoleDurbin2026", "Carlisle2024"): add(cg + ", Metazoa " + lab, "clock study", set_name="C3", fix={"animals": lab}) add(cg + ", youngest dates (Kay LECA, Carlisle Metazoa)", "clock study", set_name="C3", fix={"eukaryotes": "Kay2026", "animals": "Carlisle2024"}) add(cg + ", eukaryotes = mitochondrial endosymbiosis", "clock study", set_name="C3", euk_def="mitochondria") add(cg + ", timing-study inputs (Kay LECA; animals uniform 0.574-0.80 Ga)", "clock study", set_name="C3", fix={"eukaryotes": "Kay2026", "animals": (0.574, 0.80)}) add(cg + ", timing-study inputs and its lag prior (log-uniform 10-3000 Myr)", "clock study", set_name="C3", fix={"eukaryotes": "Kay2026", "animals": (0.574, 0.80)}, tau_prior="logU_10_3000") # --- which event is the animal gate (eukaryote and human gates unchanged) ------------------------------ for lab, key in (("Tonian onset 0.75-0.85 Ga", "tonian_onset"), ("0.80 Ga exactly", "t800"), ("Marinoan termination 0.635 Ga", "marinoan_end"), ("Ediacaran ocean oxygenation 0.55-0.56 Ga (Dahl et al. 2010)", "ediacaran_ox"), ("Great Oxidation (Mills et al. 2025: animals excluded only before ~2.2 Ga)", "goe_permanent")): add("animal gate = " + lab, "animal gate", set_name="C3", animal_gate=key) for tp, txt in (("logU_10_2000", "10-2000 Myr"), ("logU_10_5000", "10-5000 Myr")): add("animal gate = Great Oxidation, lag prior log-uniform " + txt, "animal gate", set_name="C3", animal_gate="goe_permanent", tau_prior=tp) add("animal gate uncertain: equal mix of Neoproterozoic window, Tonian onset and Great Oxidation", "animal gate", set_name="C3", animal_gate=ANIMAL_GATE_MIX) # --- lag prior and remaining habitable time ----------------------------------------------------------- for tp in ("logU_1_1000", "logU_10_2000", "U_0_500"): add(cg + ", lag prior " + tp, "lag prior / window", set_name="C3", tau_prior=tp) for fut in ((0.9,), (1.73,), (0.9, 1.08, 1.73)): add(cg + ", future %s" % (fut,), "lag prior / window", set_name="C3", future=fut, batches=2 * batches) # --- physiological thresholds fixed in advance -------------------------------------------------------- for mp in ("all", "low", "high"): for an in ("animal_mills25", "animal_sponge", "animal_catling", "animal_catling_circ"): add("physiological gates, mid-Proterozoic %s, %s" % (mp, an), "threshold gates", set_name="C3", gate_version="physio", physio_kw=dict(midprot=mp, animal=an, human="human")) add("physiological gates, all, sponge, fire threshold for humans", "threshold gates", set_name="C3", gate_version="physio", physio_kw=dict(midprot="all", animal="animal_sponge", human="fire")) add("physiological gates, all, sponge, first crossing for eukaryotes", "threshold gates", set_name="C3", gate_version="physio", physio_kw=dict(midprot="all", animal="animal_sponge", human="human", euk_first=True)) add("physiological gates, all, sponge, mitochondria", "threshold gates", set_name="C3", gate_version="physio", euk_def="mitochondria", physio_kw=dict(midprot="all", animal="animal_sponge", human="human")) add("physiological gates, all, sponge, C4", "threshold gates", set_name="C4", gate_version="physio", physio_kw=dict(midprot="all", animal="animal_sponge", human="human")) df = pd.DataFrame(rows) df.to_csv(os.path.join(DATA, "part1_bayes_factors.csv"), index=False) return df def headline(gcso, D=12000, seed=SEED + 5, batches=8): """The reference configuration (three gated transitions, conventional oxygenation-event gates, primary lag prior, equal-weight mixture of clock studies) on several independent batches: pooled Bayes factors, their batch-to-batch scatter, and diagnostics of WHERE the gating likelihood comes from (which clock study, which draws). The referee's point: the average is dominated by one LECA study and by near-zero-lag draws.""" rng = np.random.default_rng(seed) per, parts = [], [] lab_e, lab_a, minlag = [], [], [] for b in range(batches): dr = make_draws(rng, D, SETS["C3"], gate_version="conventional", gcso=gcso) L_H, L_Hg, L_G = pure_models(dr) parts.append((L_H, L_Hg, L_G)) per.append([L_G.mean() / L_H.mean(), L_G.mean() / L_Hg.mean(), L_Hg.mean() / L_H.mean()]) lab_e.append(dr["labels"]["eukaryotes"]); lab_a.append(dr["labels"]["animals"]) t, g = dr["t"], dr["g"] lag_e = t[:, 0] - np.maximum(dr["t0"], g[:, 0]) lag_a = t[:, 1] - np.maximum(t[:, 0], g[:, 1]) minlag.append(np.minimum(lag_e, lag_a)) per = np.array(per) L_H, L_Hg, L_G = (np.concatenate([p[i] for p in parts]) for i in range(3)) lab_e, lab_a, minlag = np.concatenate(lab_e), np.concatenate(lab_a), np.concatenate(minlag) pooled = [L_G.mean() / L_H.mean(), L_G.mean() / L_Hg.mean(), L_Hg.mean() / L_H.mean()] tot = L_G.sum() e_names = [e[0] for e in classic.TRANSITIONS["eukaryotes"]["estimates"]] a_names = [e[0] for e in classic.TRANSITIONS["animals"]["estimates"]] near = (L_G > 0) & (minlag < 0.02) top = np.sort(L_G)[-max(1, len(L_G) // 100):] out = dict(n_batches=batches, draws_per_batch=D, BF_G_H=float(pooled[0]), BF_G_Hg=float(pooled[1]), BF_Hg_H=float(pooled[2]), batch_rel_sd_G_H=float(per[:, 0].std(ddof=1) / per[:, 0].mean()), batch_rel_sd_G_Hg=float(per[:, 1].std(ddof=1) / per[:, 1].mean()), batch_rel_sd_Hg_H=float(per[:, 2].std(ddof=1) / per[:, 2].mean()), pooled_rel_se_G_H=float(per[:, 0].std(ddof=1) / per[:, 0].mean() / np.sqrt(batches)), pooled_rel_se_G_Hg=float(per[:, 1].std(ddof=1) / per[:, 1].mean() / np.sqrt(batches)), batch_values_G_H=[float(x) for x in per[:, 0]], frac_draws_G_zero=float(np.mean(L_G == 0)), share_of_gating_likelihood_by_LECA_study={n: float(L_G[lab_e == k].sum() / tot) for k, n in enumerate(e_names)}, share_of_gating_likelihood_by_Metazoa_study={n: float(L_G[lab_a == k].sum() / tot) for k, n in enumerate(a_names)}, share_of_gating_likelihood_top_1pct_draws=float(top.sum() / tot), share_of_gating_likelihood_lag_under_20Myr=float(L_G[near].sum() / tot), frac_draws_lag_under_20Myr=float(near.mean()), BF_G_H_by_LECA_study={n: float(L_G[lab_e == k].mean() / L_H[lab_e == k].mean()) for k, n in enumerate(e_names)}, BF_G_Hg_by_LECA_study={n: float(L_G[lab_e == k].mean() / L_Hg[lab_e == k].mean()) for k, n in enumerate(e_names)}) with open(os.path.join(DATA, "part1_headline.json"), "w") as f: json.dump(out, f, indent=1) return out def per_step_factors(gcso, D=20000, seed=SEED + 1): """Single-step view: for each gated transition, the lag distribution and the ratio of the gating density (lag marginalised over the prior on its mean) to the gated-hard density 1 / (available interval), where the available interval runs from the start (gate or previous step) to the next transition (or the end of the window for the last step). A transparent decomposition; the full chain result is in the BF table.""" rng = np.random.default_rng(seed) rows = [] versions = [("conventional", None), ("physio", dict(midprot="low", animal="animal_sponge")), ("physio", dict(midprot="high", animal="animal_sponge")), ("physio", dict(midprot="all", animal="animal_sponge")), ("physio", dict(midprot="all", animal="animal_mills25")), ("physio", dict(midprot="all", animal="animal_catling"))] versions = [v + ("noe_window",) for v in versions] + [("conventional", None, "tonian_onset"), ("conventional", None, "goe_permanent")] for gv, kw, ag in versions: dr = make_draws(rng, D, SETS["C3"], gate_version=gv, physio_kw=kw, gcso=gcso, animal_gate=ag) for i, s in enumerate(dr["steps"]): lag_age = dr["gate_age"][s] - dr["ages"][s] # Gyr, positive = after the gate prev = dr["t0"] if i == 0 else dr["t"][:, i - 1] start = np.maximum(prev, dr["g"][:, i]) lag = dr["t"][:, i] - start nxt = dr["T"] if i == len(dr["steps"]) - 1 else dr["t"][:, i + 1] avail = nxt - start ok = (lag >= 0) & (avail > 0) lg = np.where(ok, np.maximum(lag, 1e-6), 1.0) dens_g = np.where(ok, chain.lag_density_marginal(lg, dr["prior"]), 0.0) dens_h = np.where(ok, 1.0 / np.where(avail > 0, avail, 1.0), 0.0) q = np.quantile(lag_age, [0.16, 0.5, 0.84]) rows.append(dict(gate_version=gv, physio=json.dumps(kw or {}), animal_gate=ag if gv == "conventional" else "", transition=s, gate_age_median_Ga=float(np.median(dr["gate_age"][s])), date_median_Ga=float(np.median(dr["ages"][s])), lag_p16_Gyr=q[0], lag_median_Gyr=q[1], lag_p84_Gyr=q[2], P_transition_before_gate=float(np.mean(lag_age < 0)), lag_over_available_median=float(np.median((lag / avail)[ok])) if ok.any() else np.nan, step_factor_gating_over_gatedhard=float(dens_g.mean() / dens_h.mean()) if dens_h.mean() > 0 else np.nan)) df = pd.DataFrame(rows) df.to_csv(os.path.join(DATA, "part1_lags.csv"), index=False) return df def subset_posterior(gcso, D=8000, seed=SEED + 2): """Posterior over which, and how many, of the four candidate steps are hard. Two priors on the number of hard steps are applied to the same likelihoods: 'candidate' : 0.30/0.30/0.22/0.12/0.06, the model's prior on the number of CANDIDATE steps, read as if every candidate were hard (the first pass used only this one); 'implied' : the model's prior on the number of HARD steps, P(0) = 0.35, mean 1.15 (IMPLIED_PRIOR_K). The transferable quantity is the Bayes factor 'no hard steps : at least one', which is given for both.""" rng = np.random.default_rng(seed) cands = SETS["C4"] rows, summ = [], [] cg = "conventional gates" configs = [ (cg + ", hard steps gated, photosynthesis flat", dict(), "G", "flat"), (cg + ", hard steps ungated, photosynthesis flat", dict(), "H", "flat"), (cg + ", hard steps gated, photosynthesis strict", dict(), "G", "strict"), ("threshold gates (mid-Proterozoic 0.1-10% PAL, sponge), hard gated, flat", dict(gate_version="physio", physio_kw=dict(midprot="all", animal="animal_sponge")), "G", "flat"), ("threshold gates (mid-Proterozoic 0.1-1% PAL, sponge), hard gated, flat", dict(gate_version="physio", physio_kw=dict(midprot="low", animal="animal_sponge")), "G", "flat"), ("threshold gates (mid-Proterozoic 1-10% PAL, sponge), hard gated, flat", dict(gate_version="physio", physio_kw=dict(midprot="high", animal="animal_sponge")), "G", "flat"), ("threshold gates (0.1-10% PAL, Mills 2025 animal threshold), hard gated, flat", dict(gate_version="physio", physio_kw=dict(midprot="all", animal="animal_mills25")), "G", "flat"), ("threshold gates (1-10% PAL, Mills 2025 animal threshold), hard gated, flat", dict(gate_version="physio", physio_kw=dict(midprot="high", animal="animal_mills25")), "G", "flat"), ("animal gate = Great Oxidation (Mills et al. 2025 reading), hard gated, flat", dict(animal_gate="goe_permanent"), "G", "flat"), ("animal gate = Great Oxidation, lag prior logU 10-2000 Myr, hard gated, flat", dict(animal_gate="goe_permanent", tau_prior="logU_10_2000"), "G", "flat"), ("animal gate = Great Oxidation, lag prior logU 10-5000 Myr, hard gated, flat", dict(animal_gate="goe_permanent", tau_prior="logU_10_5000"), "G", "flat"), ("animal gate = Tonian onset 0.75-0.85 Ga, hard gated, flat", dict(animal_gate="tonian_onset"), "G", "flat"), ("animal gate uncertain (equal mix of three readings), hard gated, flat", dict(animal_gate=ANIMAL_GATE_MIX), "G", "flat"), (cg + ", lag prior logU 1-1000 Myr, hard gated, flat", dict(tau_prior="logU_1_1000"), "G", "flat"), (cg + ", lag prior logU 10-2000 Myr, hard gated, flat", dict(tau_prior="logU_10_2000"), "G", "flat"), (cg + ", LECA Strassert only, hard gated, flat", dict(fix={"eukaryotes": "Strassert2021"}), "G", "flat"), (cg + ", youngest dates (Kay LECA, Carlisle Metazoa), hard gated, flat", dict(fix={"eukaryotes": "Kay2026", "animals": "Carlisle2024"}), "G", "flat"), (cg + ", mitochondria as the eukaryote event, hard gated, flat", dict(euk_def="mitochondria"), "G", "flat"), ] priors = {"": MODEL_PRIOR_K, "imp_": IMPLIED_PRIOR_K} def bf0(Lk, pri): return float(Lk[0] / ((pri[1:] * Lk[1:]).sum() / pri[1:].sum())) for label, kw, hk, ung in configs: dr = make_draws(rng, D, cands, gcso=gcso, **kw) Ls = {} for r in range(len(cands) + 1): for S in itertools.combinations(cands, r): Ls[S] = lik(dr, kinds_for(cands, set(S), hk, ung)) mean = {S: float(v.mean()) for S, v in Ls.items()} Lk = np.array([np.mean([mean[S] for S in mean if len(S) == k]) for k in range(len(cands) + 1)]) keys = list(Ls) mat = np.column_stack([Ls[S] for S in keys]) sizes = np.array([len(S) for S in keys]) idx = rng.integers(0, D, size=(300, D)) lk_bs = [] for b in range(300): mb = mat[idx[b]].mean(axis=0) lk_bs.append(np.array([mb[sizes == k].mean() for k in range(len(cands) + 1)])) lk_bs = np.array(lk_bs) for S in keys: rows.append(dict(config=label, hard_subset="+".join(S) if S else "(none)", n_hard=len(S), marginal_likelihood=mean[S], relative_to_all_gated=mean[S] / mean[()] if mean[()] > 0 else np.nan)) rec = dict(config=label, photosynthesis=ung, **{"L_k%d_rel" % k: Lk[k] / Lk[0] for k in range(5)}) for tag, pri in priors.items(): post = pri * Lk post = post / post.sum() bs = pri[None, :] * lk_bs bs = bs / bs.sum(axis=1, keepdims=True) rec.update(**{tag + "post_k%d" % k: post[k] for k in range(5)}, **{tag + "post_k%d_mcse" % k: float(bs[:, k].std()) for k in range(5)}, **{tag + "prior_k0": float(pri[0]), tag + "prior_mean_k": float((pri * np.arange(5)).sum()), tag + "post_mean_k": float((post * np.arange(5)).sum()), tag + "BF_k0_vs_kge1": bf0(Lk, pri), tag + "BF_k0_vs_kge1_lnse": float(np.std([np.log(bf0(x, pri)) for x in lk_bs]))}) # per-step marginal: posterior probability that each candidate is hard, under the candidate prior on k # with subsets of equal size equally likely a priori pri = {S: MODEL_PRIOR_K[len(S)] / sum(1 for X in keys if len(X) == len(S)) for S in keys} w = np.array([pri[S] * mean[S] for S in keys]) w = w / w.sum() step_post = {c: float(sum(w[i] for i, S in enumerate(keys) if c in S)) for c in cands} step_prior = {c: float(sum(pri[S] for S in keys if c in S)) for c in cands} rec.update(**{"P_hard_" + c: step_post[c] for c in cands}, **{"prior_hard_" + c: step_prior[c] for c in cands}) summ.append(rec) pd.DataFrame(rows).to_csv(os.path.join(DATA, "part1_subset_likelihoods.csv"), index=False) sm = pd.DataFrame(summ) sm.to_csv(os.path.join(DATA, "part3_posterior_hard_steps.csv"), index=False) return sm def _pool_density(t, t0, A0, gates_age, prior): """Lag density of every step under every candidate gate, with the chain rule of the main model: a step starts at max(previous transition, gate). t: (D, n) forward times; gates_age: (D, M) ages in Ga. Returns (n, D, M) in 1/Gyr; zero where the transition precedes its start. (The chain normaliser, which is within about 2 per cent of 1 for easy steps, is omitted here.)""" gate_t = A0[:, None] - gates_age out = [] for i in range(t.shape[1]): prev = t0 if i == 0 else t[:, i - 1] lag = t[:, i][:, None] - np.maximum(prev[:, None], gate_t) out.append(np.where(lag > 0, chain.lag_density_marginal(np.where(lag > 0, lag, 1.0), prior), 0.0)) return np.array(out) def _assignment_table(dens): """dens (3, D, M) -> (M, M, M) marginal likelihood of every assignment of one pool event to each step.""" D, M = dens.shape[1], dens.shape[2] pair = (dens[1][:, :, None] * dens[2][:, None, :]).reshape(D, M * M) return (dens[0].T @ pair).reshape(M, M, M) / D def _hard_step_fakes(rng, sub, blurred): """One synthetic history in which the first two transitions are hard steps: uniform order statistics between the start of the chain and today, the third (humans) kept at today. blurred=False: the fake dates are known exactly (first-pass version). blurred=True : each fake date carries the same dating-error distribution as the real one (the real draws re-centred on the fake date), so the fake record is dated as loosely as the real record.""" u = np.sort(rng.uniform(size=2)) span = sub["A0"] - sub["t0"] tt = sub["t"].copy() for i in range(2): tt[:, i] = sub["t0"] + u[i] * span if blurred: tt[:, i] += sub["t"][:, i] - np.median(sub["t"][:, i]) return tt def flexibility(gcso, D=12000, seed=SEED + 3, nnull=300, observed=None): """How much of the gating fit comes from the freedom to choose and date the gates. `observed` is the pooled reference Bayes factor (gating : hard steps) against which the nulls are compared.""" rng = np.random.default_rng(seed) out = {} steps = SETS["C3"] dr = make_draws(rng, D, steps, gate_version="conventional", gcso=gcso) L_H, L_Hg, L_G = pure_models(dr) out["this_batch_BF_G_H"] = float(L_G.mean() / L_H.mean()) obs = float(observed) if observed is not None else out["this_batch_BF_G_H"] out["observed_BF_G_H"] = obs # (a) shuffle the gates among the transitions perm_rows = [] for perm in itertools.permutations(range(len(steps))): d2 = dict(dr) d2["g"] = dr["g"][:, list(perm)] LG = lik(d2, ["E"] * len(steps)) perm_rows.append(dict(assignment=" ".join("%s<-gate of %s" % (steps[i], steps[p]) for i, p in enumerate(perm)), identity=perm == tuple(range(len(steps))), BF_G_H=float(LG.mean() / L_H.mean()), frac_zero=float(np.mean(LG == 0)))) out["shuffle"] = perm_rows # (b) random gate dates: three dates uniform over the window before today, sorted, assigned in order Dn = 2000 sub = {k: (v[:Dn] if isinstance(v, np.ndarray) else v) for k, v in dr.items()} LH_sub = L_H[:Dn].mean() null_bf = [] for r in range(nnull): ga = np.sort(rng.uniform(0.0, 4.4, size=3))[::-1] # ages, oldest first d2 = dict(sub) d2["g"] = sub["A0"][:, None] - ga[None, :] LG = lik(d2, ["E"] * 3) null_bf.append(float(LG.mean() / LH_sub)) null_bf = np.array(null_bf) out["random_gates"] = dict(n=nnull, frac_ge_observed=float(np.mean(null_bf >= obs)), frac_positive=float(np.mean(null_bf > 0)), median=float(np.median(null_bf)), p90=float(np.quantile(null_bf, 0.9)), p99=float(np.quantile(null_bf, 0.99)), max=float(null_bf.max())) # (c) transitions drawn from the hard-steps model (uniform order statistics before today), conventional gates for blurred in (False, True): null2 = [] for r in range(nnull): d2 = dict(sub) d2["t"] = _hard_step_fakes(rng, sub, blurred) LG = lik(d2, ["E"] * 3) LH = lik(d2, ["H"] * 3) null2.append(float(LG.mean() / LH.mean()) if LH.mean() > 0 else 0.0) null2 = np.array(null2) out["hard_step_null_" + ("blurred_dates" if blurred else "exact_dates")] = dict( n=nnull, frac_ge_observed=float(np.mean(null2 >= obs)), frac_positive=float(np.mean(null2 > 0)), median=float(np.median(null2)), p90=float(np.quantile(null2, 0.9)), p99=float(np.quantile(null2, 0.99))) # (d) candidate pool: each transition's gate is one of M listed step changes conv_age = np.column_stack([dr["gate_age"][s] for s in steps]) conv = np.prod([_pool_density(dr["t"], dr["t0"], dr["A0"], conv_age[:, [i]], dr["prior"])[i, :, 0] for i in range(3)], axis=0) out["conventional_gates_same_unnormalised_likelihood_BF"] = float(conv.mean() / L_H.mean()) pool_rows = [] for pool_name, pool in (("oxygen events only", classic.POOL_OXYGEN), ("all listed events", classic.POOL_ALL)): M = len(pool) gates_age = np.column_stack([classic.draw_env(rng, k, D) for k in pool]) # (D, M) dens = _pool_density(dr["t"], dr["t0"], dr["A0"], gates_age, dr["prior"]) # (3, D, M) LH = L_H.mean() E = _assignment_table(dens) best = np.unravel_index(np.argmax(E), E.shape) # events that precede the transition in at least half of the date draws ('possible' gates) feas = [np.mean(dens[i] > 0, axis=0) >= 0.5 for i in range(3)] feas_mix = np.prod([dens[i][:, feas[i]].mean(axis=1) for i in range(3)], axis=0) row = dict(pool=pool_name, M=M, ceiling_M_cubed=M ** 3, BF_equal_prior_pool_vs_H=float(E.mean() / LH), n_possible_gates_per_step={s: int(feas[i].sum()) for i, s in enumerate(steps)}, BF_equal_prior_over_possible_gates_vs_H=float(feas_mix.mean() / LH), BF_best_single_assignment_vs_H=float(E.max() / LH), best_single_assignment={s: pool[best[i]] for i, s in enumerate(steps)}, flexibility_factor_best_assignment=float(E.max() / E.mean()), BF_per_draw_maximum_vs_H=float(np.prod(dens.max(axis=2), axis=0).mean() / LH), mean_density_per_Gyr={s: {pool[j]: float(dens[i, :, j].mean()) for j in range(M)} for i, s in enumerate(steps)}, n_assignments_with_zero_likelihood=int((E == 0).sum()), n_assignments_ge_conventional=int((E / LH >= out["conventional_gates_same_unnormalised_likelihood_BF"]).sum())) # calibration: the same post hoc procedure (best single assignment) applied to hard-step histories ga_sub = gates_age[:Dn] for blurred in (False, True): cal = [] for r in range(nnull): tt = _hard_step_fakes(rng, sub, blurred) ordered = (tt[:, 0] > sub["t0"]) & (tt[:, 1] > tt[:, 0]) & (tt[:, 2] > tt[:, 1]) LHf = (6.0 / (sub["T"] - sub["t0"]) ** 3 * ordered).mean() Ef = _assignment_table(_pool_density(tt, sub["t0"], sub["A0"], ga_sub, dr["prior"])) cal.append(float(Ef.max() / LHf) if LHf > 0 else 0.0) cal = np.array(cal) tag = "null_blurred_" if blurred else "null_exact_" row.update({tag + "median": float(np.median(cal)), tag + "p90": float(np.quantile(cal, 0.9)), tag + "frac_ge_observed": float(np.mean(cal >= row["BF_best_single_assignment_vs_H"]))}) pool_rows.append(row) out["pool"] = pool_rows with open(os.path.join(DATA, "part1_flexibility.json"), "w") as f: json.dump(out, f, indent=1) return out def spacing_test(D=20000, seed=SEED + 4): """Carter's corollary: under hard steps the n+1 intervals (including the time remaining) are exchangeable and Dirichlet(1,...,1). Statistics: coefficient of variation of the intervals, and the remaining time.""" rng = np.random.default_rng(seed) rows = [] for set_name in ("C4", "C5", "C6"): for fut in (0.9, 1.0, 1.73): dr = make_draws(rng, D, SETS[set_name], gate_version="none", future=(fut,)) t = dr["t"] ok = np.all(np.diff(np.column_stack([dr["t0"], t, dr["T"]]), axis=1) > 0, axis=1) edges = np.column_stack([dr["t0"], t, dr["T"]])[ok] gaps = np.diff(edges, axis=1) frac = gaps / gaps.sum(axis=1, keepdims=True) n = t.shape[1] cv_obs = frac.std(axis=1) / frac.mean(axis=1) sim = rng.dirichlet(np.ones(n + 1), size=200000) cv_null = sim.std(axis=1) / sim.mean(axis=1) p_low = np.mean(cv_null[None, :2000] <= cv_obs[:2000, None]) # P(null CV <= observed), averaged rem_obs = frac[:, -1] p_rem = np.mean(sim[:2000, -1][None, :] <= rem_obs[:2000, None]) rows.append(dict(set=set_name, n_steps=n, future_Gyr=fut, frac_ordered_draws=float(ok.mean()), mean_gap_Gyr=float(gaps.mean()), sd_gap_Gyr=float(gaps.std(axis=1).mean()), cv_observed_median=float(np.median(cv_obs)), cv_null_median=float(np.median(cv_null)), P_null_cv_le_observed=float(p_low), remaining_fraction=float(np.median(rem_obs)), expected_remaining_fraction=1.0 / (n + 1), P_null_remaining_le_observed=float(p_rem))) df = pd.DataFrame(rows) df.to_csv(os.path.join(DATA, "part1_spacing_test.csv"), index=False) return df def run_all(): _, gcso = load_gcso() classic.write_tables() res = {} res["selftest"] = chain.selftest() res["headline"] = headline(gcso) res["bf"] = bayes_factor_table(gcso) res["lags"] = per_step_factors(gcso) res["post"] = subset_posterior(gcso) res["flex"] = flexibility(gcso, observed=res["headline"]["BF_G_H"]) res["spacing"] = spacing_test() return res if __name__ == "__main__": import time t0 = time.time() pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70) r = run_all() print(json.dumps(r["headline"], indent=1)) print(r["bf"][["config", "set", "BF_G_H", "lnse_G_H", "lnse_batch_G_H", "BF_G_Hg", "lnse_G_Hg", "BF_Hg_H", "frac_draws_G_zero", "frac_luca_redrawn"]].round(3).to_string()) print(r["lags"].round(3).to_string()) print(r["post"].round(3).T.to_string()) print(json.dumps(r["flex"], indent=1)) print(r["spacing"].round(3).to_string()) print("seconds", time.time() - t0)