"""Analysis 4: hard steps versus geobiological pacing (Mills et al. 2025) as generative models. Data used (central values and intervals in config.py / data/events.csv), times measured from t0: transitions: life (only bounded: before LUCA), eukaryotes (LECA), animals (crown Metazoa), Homo sapiens (now) gates: habitability (t0), Great Oxidation (GOE), Neoproterozoic oxygenation (NOE), Paleozoic oxygenation (POE) window: T = t0 + future habitable time Three generative models for the K transition times, every one conditioned on all K transitions finishing before T (we are here): H ungated hard steps: every rate << 1/T -> times are uniform order statistics, density K!/T^K on the ordered simplex. Gates play no role. GH gated hard steps (the Watson 2025 picture): hard steps, but step j cannot be attempted before gate e_j opens -> uniform density on {t_1<..e_j}, i.e. 1/V(e,T). P pacing: transition j follows gate e_j after a lag ~ Exp(mean mu); mu << T means "happens as soon as permitted". mu gets a log-uniform prior (shared, or one per transition). All three are evaluated on the same data {t_1 <= t_L', t_2, t_3, t_4}, conditional on the gate times. Part (b) then asks what happens when the gate times are NOT handed to the models for free. """ from __future__ import annotations import os import numpy as np import pandas as pd import config as cfg MU_GRID = 10.0 ** np.linspace(-3, 1, 801) # Gyr # --------------------------------------------------------------------------------- def region_volume(gates, T): """Volume of {t_1<...e_j} for sorted gates e (M,K), T (M,). g_0 = 1, g_j(t) = int_{e_j}^{t} g_{j-1}(s) ds, V = g_K(T); polynomial coefficients per draw.""" gates = np.atleast_2d(gates) M, K = gates.shape c = np.zeros((M, K + 1)) c[:, 0] = 1.0 for j in range(K): new = np.zeros_like(c) new[:, 1:] = c[:, :-1] / np.arange(1, K + 1)[None, :] e = gates[:, j] const = sum(new[:, k] * e ** k for k in range(1, K + 1)) new[:, 0] = -const c = new T = np.broadcast_to(np.asarray(T, float), (M,)) return sum(c[:, k] * T ** k for k in range(K + 1)) def draw_dates(rng, M, fut, central=False, F_fixed=None, met_age=None, noe_age=None): if central: t0 = np.full(M, cfg.T0_ZIRCON) F = np.full(M, F_fixed) luca = np.full(M, cfg.LUCA_PAPER[0]) leca = np.full(M, cfg.LECA[0]) met = np.full(M, cfg.METAZOA[0] if met_age is None else met_age) goe, poe = np.full(M, cfg.GOE[0]), np.full(M, cfg.POE[0]) noe = np.full(M, cfg.NOE[0] if noe_age is None else noe_age) else: t0 = cfg.draw_t0(rng, M) F = cfg.draw_future(rng, M, fut) luca = np.minimum(cfg.split_normal(rng, *cfg.LUCA_PAPER, M), t0 - 0.005) leca, met = cfg.draw_uniform(rng, cfg.LECA, M), cfg.draw_uniform(rng, cfg.METAZOA, M) goe, noe, poe = (cfg.draw_uniform(rng, cfg.GOE, M), cfg.draw_uniform(rng, cfg.NOE, M), cfg.draw_uniform(rng, cfg.POE, M)) T = t0 + F return dict(t0=t0, T=T, tLp=t0 - luca, trans=np.stack([t0 - leca, t0 - met, t0], axis=1), # steps 2..4 gates=np.stack([np.zeros(M), t0 - goe, t0 - noe, t0 - poe], axis=1)) def pacing_like_mu(d, mu, use_life=True): """(M, n_mu) likelihood of the transitions under P for each mu on the grid (shared mu).""" lags = d["trans"] - d["gates"][:, 1:] # (M,3) ok = (lags >= 0).all(axis=1) S = lags.sum(axis=1) L = np.exp(-S[:, None] / mu[None, :]) / mu[None, :] ** 3 if use_life: L = L * (-np.expm1(-d["tLp"][:, None] / mu[None, :])) # selection: the last transition must happen before T L = L / (-np.expm1(-(d["T"] - d["gates"][:, 3])[:, None] / mu[None, :])) return L * ok[:, None], ok def pacing_like_indep(d, mu, w, use_life=True): """Independent mu_j per transition, each with prior weights w on the grid: (M,) likelihood.""" lags = d["trans"] - d["gates"][:, 1:] ok = (lags >= 0).all(axis=1) out = np.ones(len(ok)) for j in range(3): f = np.exp(-np.clip(lags[:, j], 0, None)[:, None] / mu[None, :]) / mu[None, :] if j == 2: f = f / (-np.expm1(-(d["T"] - d["gates"][:, 3])[:, None] / mu[None, :])) out = out * (f @ w) if use_life: out = out * ((-np.expm1(-d["tLp"][:, None] / mu[None, :])) @ w) return out * ok def mu_prior(lo, hi): w = ((MU_GRID >= lo) & (MU_GRID <= hi)).astype(float) return w / w.sum() def conditional_on_gates(rng, M=40_000): rows = [] cases = [("central dates, future 0.9 Gyr", dict(central=True, F_fixed=0.9, fut=None)), ("central dates, future 1.73 Gyr", dict(central=True, F_fixed=1.73, fut=None)), ("dates marginalised, future 'short'", dict(fut="short")), ("dates marginalised, future 'long'", dict(fut="long")), ("dates marginalised, future 'mix'", dict(fut="mix"))] for lab, kw in cases: m = 4 if kw.get("central") else M d = draw_dates(rng, m, **kw) for use_life in (True, False): K = 4 if use_life else 3 gates = d["gates"] if use_life else d["gates"][:, 1:] V = region_volume(gates, d["T"]) first = d["tLp"] if use_life else 1.0 LH = (first * np.array([1, 1, 2, 6, 24])[K] / d["T"] ** K).mean() lags = d["trans"] - d["gates"][:, 1:] gated_ok = (lags >= 0).all(axis=1) # GH also forbids a step before its gate LGH = (first / V * gated_ok).mean() p_feas = (V / (d["T"] ** K / np.array([1, 1, 2, 6, 24])[K])).mean() Lmu, ok = pacing_like_mu(d, MU_GRID, use_life) Lmu_mean = Lmu.mean(axis=0) row = dict(case=lab, transitions="life + eukaryotes + animals + humans" if use_life else "eukaryotes + animals + humans", K=K, frac_date_draws_with_every_transition_after_its_gate=float(ok.mean()), p_feasible_under_H=float(p_feas), BF_GH_over_H=float(LGH / LH), mu_best_Gyr=float(MU_GRID[np.argmax(Lmu_mean)]), BF_P_over_H_best_mu=float(Lmu_mean.max() / LH)) for plab, (lo, hi) in {"mu~LU[0.01,3]": (0.01, 3.0), "mu~LU[0.001,10]": (0.001, 10.0), "mu~LU[0.03,1]": (0.03, 1.0)}.items(): w = mu_prior(lo, hi) LP = float(Lmu_mean @ w) row[f"BF_P_over_H | {plab}"] = LP / LH row[f"BF_P_over_GH | {plab}"] = LP / LGH LPi = pacing_like_indep(d, MU_GRID, mu_prior(0.01, 3.0), use_life).mean() row["BF_P_over_H | independent mu_j~LU[0.01,3]"] = float(LPi / LH) row["BF_P_over_GH | independent mu_j~LU[0.01,3]"] = float(LPi / LGH) rows.append(row) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a4_bayes_factors_given_gates.csv"), index=False) return df # --------------------------------------------------------------------------------- def sample_gh(rng, gates, T, N): """Rejection-sample N draws from the uniform distribution on {t_1<..e_j}.""" K = len(gates) out = [] got = 0 while got < N: u = np.sort(rng.random((200_000, K)), axis=1) * T keep = u[(u > np.asarray(gates)[None, :]).all(axis=1)] out.append(keep) got += len(keep) return np.concatenate(out)[:N] def lag_test_under_gh(rng): """Frequentist check of GH: are the observed lags after the gates surprisingly short?""" out = {} for F in (0.9, 1.73): d = draw_dates(rng, 1, None, central=True, F_fixed=F) gates, T = d["gates"][0], d["T"][0] sim = sample_gh(rng, gates, T, 300_000) lag = sim - gates[None, :] obs_lag = d["trans"][0] - gates[1:] s_obs = obs_lag.sum() s_sim = lag[:, 1:].sum(axis=1) out[f"future {F}"] = dict( observed_lags_Gyr=[float(x) for x in obs_lag], observed_sum_Gyr=float(s_obs), GH_mean_sum_Gyr=float(s_sim.mean()), p_sum_of_3_lags_this_small=float((s_sim <= s_obs).mean()), p_also_life_before_LUCA=float(((s_sim <= s_obs) & (sim[:, 0] <= d["tLp"][0])).mean()), p_life_before_LUCA_alone=float((sim[:, 0] <= d["tLp"][0]).mean())) return out def biospheres_needed(rng, nsim=1500, mu_true=0.3): """How many independent Earth-like biospheres (same gates, all four transition times and the gate times known exactly) are needed before P and GH can be told apart at |ln BF| > ln 100?""" w = mu_prior(0.01, 3.0) out = {} for F in (0.9, 1.73): d = draw_dates(rng, 1, None, central=True, F_fixed=F) gates, T = d["gates"][0], d["T"][0] V = float(region_volume(gates[None, :], T)[0]) res = {} Nmax = 8 for truth in ("P", "GH"): if truth == "P": # exponential lags truncated so that every transition is complete before T lag = np.empty((nsim * Nmax, 4)) for j in range(4): room = T - gates[j] u = rng.random(nsim * Nmax) lag[:, j] = -mu_true * np.log1p(-u * (-np.expm1(-room / mu_true))) else: lag = sample_gh(rng, gates, T, nsim * Nmax) - gates[None, :] S = lag.sum(axis=1).reshape(nsim, Nmax) room4 = T - gates[3] med, frac = [], [] for N in range(1, Nmax + 1): Stot = S[:, :N].sum(axis=1) # log marginal likelihood of P (shared mu): prod over 4N lags, each planet truncated logL = (-Stot[:, None] / MU_GRID[None, :] - 4 * N * np.log(MU_GRID)[None, :] - N * np.log(-np.expm1(-room4 / MU_GRID))[None, :]) mx = logL.max(axis=1) logP = mx + np.log(np.exp(logL - mx[:, None]) @ w) lnbf = logP - N * (-np.log(V)) med.append(float(np.median(lnbf))) frac.append(float(np.mean(lnbf > np.log(100)) if truth == "P" else np.mean(lnbf < -np.log(100)))) res[truth] = dict(median_lnBF_P_over_GH_by_N=med, frac_decisive_by_N=frac, N_for_median_decisive=int(next((i + 1 for i, m in enumerate(med) if (m > np.log(100) if truth == "P" else m < -np.log(100))), -1))) out[f"future {F}"] = res return out # --------------------------------------------------------------------------------- def carter_coincidence(): """(b) When gate times are not given for free. Statistic: x = t_obs / T. H (K hard steps): p(t_obs | T) = K x^(K-1) / T P, driver time log-uniform on [tau_lo, tau_hi], conditioned on < T: p(t_obs | T) = 1 / (t_obs ln(T / tau_lo)) P, driver time uniform on (0, T): 1 / T. Bayes factor H:P = ratio of these densities.""" rows = [] t0 = cfg.T0_ZIRCON for flab, F in cfg.FUTURE.items(): T = t0 + F x = t0 / T for K in range(1, 7): row = dict(future=flab, x=x, K=K, BF_H_over_P_uniform_gate_time=K * x ** (K - 1)) for tau_lo in (0.1, 0.01, 0.001): row[f"BF_H_over_P_loguniform_gate_time_from_{tau_lo}_Gyr"] = K * x ** K * np.log(T / tau_lo) rows.append(row) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a4_carter_coincidence.csv"), index=False) return df def gate_prior_bf(rng, M=400_000): """(b') Same question for the whole set of three oxygenation gates. If the driving events are the hard steps themselves (order statistics, density 3!/T^3) versus geophysical gates whose successive waiting intervals are independent and log-uniform on [tau_lo, tau_hi] (conditioned on all three finishing before T), which predicts Earth's gate times better?""" out = {} t0 = cfg.T0_ZIRCON e = t0 - np.array([cfg.GOE[0], cfg.NOE[0], cfg.POE[0]]) d = np.diff(np.concatenate([[0.0], e])) for F in (0.9, 1.73): T = t0 + F for lo, hi in [(1e-3, 1e2), (1e-2, 1e1), (1e-1, 1e1)]: R = np.log(hi / lo) iv = np.exp(np.log(lo) + R * rng.random((M, 3))) psucc = float((iv.sum(axis=1) < T).mean()) dens_P = np.prod(1.0 / (d * R)) / psucc dens_H = 6.0 / T ** 3 out[f"future {F} | intervals log-uniform [{lo},{hi}] Gyr"] = dict( P_all_gates_open_before_T=psucc, BF_orderstats_over_loguniform=float(dens_H / dens_P)) return out def eukaryote_date_sensitivity(): """Strict gating requires every transition to follow its gate. Check alternative datings.""" t0 = cfg.T0_ZIRCON rows = [] for lab, age, gate in [("LECA 1.735 Ga vs GOE 2.33 Ga", 1.735, 2.33), ("mitochondrial FECA 2.25 Ga (Kay+ 2025) vs GOE 2.33 Ga", 2.25, 2.33), ("nuclear FECA 2.92 Ga (Kay+ 2025) vs GOE 2.33 Ga", 2.92, 2.33), ("crown Metazoa 0.60 Ga vs NOE 0.635 Ga", 0.60, 0.635), ("crown Metazoa 0.80 Ga (old clock) vs NOE 0.635 Ga", 0.80, 0.635), ("Metazoa 0.75 Ga (Lole Durbin+ 2026: 0.80-0.70) vs NOE 0.635 Ga", 0.75, 0.635), ("Metazoa 0.75 Ga (Lole Durbin+ 2026) vs NOE onset 0.80 Ga", 0.75, 0.80), ("Homo sapiens 0 Ga vs POE 0.41 Ga", 0.0, 0.41)]: rows.append(dict(pair=lab, lag_Gyr=gate - age, consistent_with_strict_gating=bool(gate - age >= 0))) return rows def gate_definition_sensitivity(): """Central-date Bayes factors (eukaryotes, animals, humans) as a function of the animal-origin date and of which point of the protracted Neoproterozoic oxygenation counts as 'the gate'. The baseline central values (animals 0.60 Ga, gate 0.635 Ga) differ by only 35 Myr, and Lole Durbin et al. (2026, Sci Adv 12:eaeg6289) now put the origin of Metazoa at ~0.80-0.70 Ga, i.e. BEFORE a 0.635 Ga gate, which both gated models forbid (referee report).""" rows = [] w = mu_prior(0.01, 3.0) for met, noe, lab in [(0.60, 0.635, "baseline: animals 0.60, gate 0.635 (mid-oxygenation)"), (0.574, 0.635, "animals at the fossil minimum 0.574, gate 0.635"), (0.75, 0.635, "animals 0.75 (Lole Durbin+ 2026), gate 0.635"), (0.60, 0.80, "animals 0.60, gate = onset 0.80"), (0.70, 0.80, "animals 0.70 (Lole Durbin+ 2026), gate = onset 0.80"), (0.75, 0.80, "animals 0.75 (Lole Durbin+ 2026), gate = onset 0.80"), (0.60, 0.539, "animals 0.60, gate = end 0.539")]: for F in (0.9, 1.73): d = draw_dates(None, 2, None, central=True, F_fixed=F, met_age=met, noe_age=noe) V = region_volume(d["gates"][:, 1:], d["T"]) LH = float((6.0 / d["T"] ** 3).mean()) lags = d["trans"] - d["gates"][:, 1:] ok = (lags >= 0).all(axis=1) LGH = float((ok / V).mean()) Lmu, _ = pacing_like_mu(d, MU_GRID, use_life=False) LP = float(Lmu.mean(axis=0) @ w) rows.append(dict(case=lab, animals_Ga=met, neoproterozoic_gate_Ga=noe, future_Gyr=F, animal_lag_after_gate_Gyr=float(lags[0, 1]), every_transition_after_its_gate=bool(ok.all()), BF_GH_over_H=LGH / LH, BF_P_over_H=LP / LH, BF_P_over_GH=(LP / LGH) if LGH > 0 else float("nan"))) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a4_gate_definition_sensitivity.csv"), index=False) return df def run(rng): out = {} df = conditional_on_gates(rng) out["given_gates"] = df.to_dict(orient="records") out["lag_test_under_GH"] = lag_test_under_gh(rng) out["biospheres_needed"] = biospheres_needed(rng) out["carter_coincidence"] = carter_coincidence().to_dict(orient="records") out["gate_time_prior"] = gate_prior_bf(rng) out["dating_sensitivity"] = eukaryote_date_sensitivity() out["gate_definition_sensitivity"] = gate_definition_sensitivity().to_dict(orient="records") return out if __name__ == "__main__": import json r = run(np.random.default_rng(cfg.SEED)) pd.set_option("display.width", 250) print(pd.DataFrame(r["given_gates"]).drop(columns=["transitions"]).round(3).T.to_string()) print(json.dumps({k: r[k] for k in ("lag_test_under_GH", "biospheres_needed", "gate_time_prior", "dating_sensitivity")}, indent=1, default=float)) print(pd.DataFrame(r["carter_coincidence"]).round(2).to_string())