"""Analysis 1: fast vs slow abiogenesis. (a) reproduce Kipping (2020, PNAS 117:11995) and Kipping (2025, Astrobiology 25:323); (b) the two likelihood plateaus and their ratio B as a function of what follows abiogenesis: the all-hard-steps value B(m), and how far B moves when the later steps are not all hard (it is neither a floor nor a ceiling in general); (c) sensitivity to the prior family and bounds; (d) sensitivity to the dates, using the actual spread of LUCA ages in the Moody et al. (2024) repository, and a proper marginalisation over date uncertainty. """ from __future__ import annotations import os import numpy as np import pandas as pd import config as cfg from timing_lib import (bmax_hardsteps, hypoexp_logpdf, kipping_bound, kipping_like, log_hypoexp_cdf, log_hypoexp_cdf2, plateau_ratio_general, plateau_ratio_one_later_step, weighted_quantile) BIG, SMALL = 1e7, 1e-9 # --------------------------------------------------------------------------------- def reproduce_2020(rng): T = 5.304 out = {} x, w = np.polynomial.legendre.leggauss(600) for name, tLp in [("optimistic_zircon_4.1Ga", 0.304), ("conservative_microfossil_3.465Ga", 0.939)]: tIp = 4.404 - tLp r = {} z10 = float(kipping_like(BIG, SMALL, tLp, tIp, T)) z00 = float(kipping_like(SMALL, SMALL, tLp, tIp, T)) r["Z10/Z00_numeric"] = z10 / z00 r["Z10/Z00_formula"] = kipping_bound(tLp, tIp, T) for prior in ["jeffreys", "uniform"]: if prior == "jeffreys": # f ~ Beta(1/2,1/2): f = sin^2(theta), theta ~ U(0, pi/2) th = (x + 1) * np.pi / 4 ww = w / 2 lam = -2 * np.log(np.cos(th)) / T else: # f ~ U(0,1) f = (x + 1) / 2 ww = w / 2 lam = -np.log1p(-f) / T l_fast = kipping_like(BIG, lam, tLp, tIp, T) @ ww l_slow = kipping_like(SMALL, lam, tLp, tIp, T) @ ww i_slow = kipping_like(lam, SMALL, tLp, tIp, T) @ ww ens = ww @ kipping_like(lam[:, None], lam[None, :], tLp, tIp, T) @ ww r[f"ZLfast/ZLslow_{prior}"] = float(l_fast / l_slow) r[f"ZIslow/Zens_{prior}"] = float(i_slow / ens) # independent Monte Carlo cross-check of the Jeffreys ensemble evidence M = 2_000_000 fL = rng.beta(0.5, 0.5, M) fI = rng.beta(0.5, 0.5, M) # evaluated in chunks of 100,000 draws: the quadrature branch of kipping_like builds an # (M x 200) array, which needed > 8 GB for M = 2e6 in one call (referee report). Same # random stream, identical numbers, peak memory < 2 GB. CH = 100_000 like = np.concatenate([np.nan_to_num(kipping_like(-np.log1p(-fL[i:i + CH]) / T, -np.log1p(-fI[i:i + CH]) / T, tLp, tIp, T)) for i in range(0, M, CH)]) ens_mc = like.mean() th = (x + 1) * np.pi / 4 lam = -2 * np.log(np.cos(th)) / T i_slow = kipping_like(lam, SMALL, tLp, tIp, T) @ (w / 2) # Diagnostic only (revision): the two published Jeffreys odds (1.638, 1.572) that we cannot # reproduce are matched if the ensemble integral leaves out f_I < 10^-2.5 = 0.0032, i.e. the # integrable spike of the Jeffreys prior at the origin (3.6% of its mass). One free number # fitted to two values: a coincidence is possible and no error in the paper is implied. eps = 10.0 ** -2.5 a = np.arcsin(np.sqrt(eps)) thI = a + (np.pi / 2 - a) * (x + 1) / 2 wI = (w / 2) * (np.pi / 2 - a) / (np.pi / 2) lamI = -2 * np.log(np.cos(thI)) / T ens_trunc = 0.0 for j in range(0, len(lam), 100): # rows in blocks to bound memory ens_trunc += (w[j:j + 100] / 2) @ kipping_like(lam[j:j + 100, None], lamI[None, :], tLp, tIp, T) @ wI r["ZIslow/Zens_jeffreys_if_ensemble_omits_fI_below_0.0032"] = float(i_slow / ens_trunc) r["ZIslow/Zens_jeffreys_MC"] = float(i_slow / ens_mc) r["ZIslow/Zens_jeffreys_MC_se"] = float(i_slow / ens_mc * like.std() / np.sqrt(M) / ens_mc) out[name] = r out["published"] = { "optimistic_zircon_4.1Ga": {"Z10/Z00": 8.73, "ZLfast/ZLslow_jeffreys": 9.538, "ZLfast/ZLslow_uniform": 9.648, "ZIslow/Zens_jeffreys": 1.638, "ZIslow/Zens_uniform": 1.474}, "conservative_microfossil_3.465Ga": {"Z10/Z00": 2.83, "ZLfast/ZLslow_jeffreys": 3.110, "ZLfast/ZLslow_uniform": 3.137, "ZIslow/Zens_jeffreys": 1.572, "ZIslow/Zens_uniform": 1.416}} return out def reproduce_2025(rng): T = 5.304 out = {"published": {"Isua_3.7Ga": 3.8, "zircon_4.1Ga": 8.7, "LUCA_4.2Ga": 13.0, "LUCA_68pct": [9.9, 18.4], "LUCA_future_1.8Gyr": 15.2, "min_over_T": 11.3}} out["Isua_3.7Ga"] = kipping_bound(0.7, 4.404 - 0.7, T) out["zircon_4.1Ga"] = kipping_bound(0.304, 4.1, T) out["LUCA_4.2Ga"] = kipping_bound(0.204, 4.2, T) out["LUCA_future_1.8Gyr"] = kipping_bound(0.204, 4.2, 4.404 + 1.8) # minimum over T (Kipping 2025 appendix) tLp, tIp = 0.204, 4.2 out["min_over_T"] = (tIp + 2 * tLp + np.sqrt(tIp ** 2 + 4 * tIp * tLp + 2 * tLp ** 2)) / (4 * tLp) out["T_crit"] = float(np.sqrt(tIp ** 2 + 4 * tIp * tLp + 2 * tLp ** 2)) # propagate the LUCA 95% interval (4.09-4.33) as a two-piece normal luca = cfg.split_normal(rng, *cfg.LUCA_PAPER, 400_000) tl = 4.404 - luca ok = tl > 0 bf = T / (2 * tl[ok]) out["LUCA_BF_quantiles_16_50_84"] = [float(q) for q in np.percentile(bf, [15.87, 50, 84.13])] out["LUCA_frac_draws_older_than_4.404"] = float(1 - ok.mean()) return out # --------------------------------------------------------------------------------- def brute_force_bmax(m, tLp, tobs, T, later_rate=1e-4): """Numerical L(lambda_L -> inf)/L(lambda_L -> 0) for abiogenesis + m slow steps, straight from hypoexponential pdfs/cdfs (no closed form), to validate bmax_hardsteps.""" lam_later = np.full((1, m), later_rate) * (1 + 0.01 * np.arange(m)) # distinct rates # fast limit: abiogenesis at t=0 -> density of observers at tobs / P(done by T) l_inf = np.exp(hypoexp_logpdf(lam_later, tobs) - log_hypoexp_cdf(lam_later, T))[0] # slow limit: lambda_L = eps eps = 1e-9 x, w = np.polynomial.legendre.leggauss(64) t = 0.5 * tLp * (x + 1) wt = 0.5 * tLp * w dens = np.array([np.exp(hypoexp_logpdf(lam_later, tobs - ti))[0] for ti in t]) num = np.sum(wt * eps * np.exp(-eps * t) * dens) lam_all = np.concatenate([[eps], lam_later[0]])[None, :] den = np.exp(log_hypoexp_cdf(lam_all, T))[0] return float(l_inf / (num / den)) def bmax_tables(): rows = [] t0 = cfg.T0_ZIRCON for life_lab, life_age in cfg.LIFE_DATES.items(): tLp = t0 - life_age for fut_lab, F in cfg.FUTURE.items(): T = t0 + F tobs = t0 # absolute-time pacing: observers appear at a planet-set time; abiogenesis only has to # precede a deadline c. Two deadlines: the GOE (life must exist well before it) and t_obs. for lab, c in [("pacing, deadline = GOE (2.33 Ga)", t0 - cfg.GOE[0]), ("pacing, deadline = today", tobs)]: rows.append(dict(structure=lab, m_later_hard_steps=0, delay_Gyr=0.0, life_evidence=life_lab, future=fut_lab, tLp=tLp, T=T, B_max=c / tLp)) for m in range(1, 7): rows.append(dict(structure=f"{m} later hard step(s)", m_later_hard_steps=m, delay_Gyr=0.0, life_evidence=life_lab, future=fut_lab, tLp=tLp, T=T, B_max=float(bmax_hardsteps(m, tLp, tobs, T)))) for m in range(1, 4): rows.append(dict(structure=f"{m} later hard step(s) + 1.5 Gyr Earth-system delay (Watson 2025)", m_later_hard_steps=m, delay_Gyr=1.5, life_evidence=life_lab, future=fut_lab, tLp=tLp, T=T, B_max=float(bmax_hardsteps(m, tLp, tobs, T, delay=1.5)))) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a1_bmax_table.csv"), index=False) # validation of the closed form val = [] for m in (1, 2, 3, 4): val.append(dict(m=m, closed_form=float(bmax_hardsteps(m, 0.204, 4.404, 5.304)), brute_force=brute_force_bmax(m, 0.204, 4.404, 5.304))) return df, val def theorem_curves(tLp=0.204, tobs=4.404, T=5.304): """L(lambda_L)/L(lambda_L -> 0) for several 'what comes after' structures (later rates -> 0).""" lam = 10.0 ** np.linspace(-4, 3.5, 301) x, w = np.polynomial.legendre.leggauss(3000) curves = {} def hard(m, D=0.0): a, Tm = tobs - D, T - D num = np.empty_like(lam) den = np.empty_like(lam) for i, l in enumerate(lam): # u = 1 - exp(-l t): int_0^{tLp} l e^{-lt} (a - t)^{m-1} dt = int_0^{u1} (a - t(u))^{m-1} du u1 = -np.expm1(-l * tLp) u = 0.5 * u1 * (x + 1) num[i] = np.sum(0.5 * u1 * w * (a + np.log1p(-u) / l) ** (m - 1)) u2 = -np.expm1(-l * Tm) u = 0.5 * u2 * (x + 1) den[i] = np.sum(0.5 * u2 * w * np.clip(Tm + np.log1p(-u) / l, 0, None) ** m) L = num / den L0 = (m + 1) * (a ** m - (a - tLp) ** m) / (m * Tm ** (m + 1)) # analytic lambda -> 0 limit return L / L0 def finite(r): """one later step that is NOT hard (rate r): g(s) = r e^(-r s), F(s) = 1 - e^(-r s)""" num = np.empty_like(lam) den = np.empty_like(lam) for i, l in enumerate(lam): u1 = -np.expm1(-l * tLp) u = 0.5 * u1 * (x + 1) num[i] = np.sum(0.5 * u1 * w * r * np.exp(-r * (tobs + np.log1p(-u) / l))) u2 = -np.expm1(-l * T) u = 0.5 * u2 * (x + 1) den[i] = np.sum(0.5 * u2 * w * (-np.expm1(-r * np.clip(T + np.log1p(-u) / l, 0, None)))) L0 = (np.exp(-r * (tobs - tLp)) - np.exp(-r * tobs)) / (T + np.expm1(-r * T) / r) return (num / den) / L0 for m in (1, 2, 4): curves[f"{m} later hard step(s)"] = hard(m) curves["1 later step at 1.4 per Gyr (not hard)"] = finite(1.4) curves["1 later hard step + 1.5 Gyr delay"] = hard(1, 1.5) for lab, c in [("pacing, deadline = GOE", cfg.T0_ZIRCON - cfg.GOE[0]), ("pacing, deadline = today", tobs)]: curves[lab] = (-np.expm1(-lam * tLp)) / (-np.expm1(-lam * c)) / (tLp / c) return lam, curves # --------------------------------------------------------------------------------- def later_process_sensitivity(): """How far the fast:slow plateau ratio moves when the steps after abiogenesis are NOT all hard. B(m) is the value in the corner where every later rate -> 0; it is not a ceiling.""" out = {} scen = {"LUCA 4.20, start 4.404, end +0.9": (0.204, 4.404, 5.304), "LUCA 4.20, start 4.404, end +1.73": (0.204, 4.404, 6.134), "LUCA 4.20, start 4.51, end +0.9": (0.310, 4.510, 5.410), "fossils 3.48, start 4.404, end +0.9": (0.924, 4.404, 5.304)} r = 10.0 ** np.linspace(-5, 3, 3201) rows = [] for lab, (tLp, tobs, T) in scen.items(): b_exact, b_int = plateau_ratio_one_later_step(r, tLp, T) i = int(np.argmax(b_exact)) e = {"all_hard_m1": float(bmax_hardsteps(1, tLp, tobs, T)), "all_hard_m5": float(bmax_hardsteps(5, tLp, tobs, T)), "one_later_step_exact_time_max": float(b_exact[i]), "rate_at_max_per_Gyr": float(r[i]), "one_later_step_exact_time_at_rate_10": float(np.interp(1.0, np.log10(r), b_exact)), "kipping_interval_limit_fast_intelligence_T_over_tLp": float(T / tLp), "pacing_deadline_today_tobs_over_tLp": float(tobs / tLp), "pacing_deadline_GOE": float((tobs - cfg.GOE[0]) / tLp)} # two later steps: one hard, one with a finite rate for rr in (0.1, 1.0, 10.0): e[f"two_later_steps_one_hard_one_rate_{rr:g}"] = float( plateau_ratio_general([1e-7, rr], tLp, tobs, T)) # two later steps, both rates free: coarse maximum g = 10.0 ** np.linspace(-2, 1.5, 36) best = max(((plateau_ratio_general([a, b], tLp, tobs, T), a, b) for k, a in enumerate(g) for b in g[k:]), key=lambda t: t[0]) e["two_later_steps_max"] = float(best[0]) e["two_later_steps_rates_at_max"] = [float(best[1]), float(best[2])] vals = [e[k] for k in ("all_hard_m1", "all_hard_m5", "one_later_step_exact_time_max", "kipping_interval_limit_fast_intelligence_T_over_tLp", "pacing_deadline_today_tobs_over_tLp", "two_later_steps_max")] e["range_over_structures"] = [float(min(vals)), float(max(vals))] out[lab] = e rows.append(dict(scenario=lab, **{k: v for k, v in e.items() if not isinstance(v, list)})) pd.DataFrame(rows).to_csv(os.path.join(cfg.DATA, "a1_later_process_sensitivity.csv"), index=False) # validation of the closed forms and of the general routine (baseline scenario) tLp, tobs, T = scen["LUCA 4.20, start 4.404, end +0.9"] val = [] for rr in (1e-3, 0.3, 1.4, 5.0): cf_exact, cf_int = plateau_ratio_one_later_step(rr, tLp, T) num_int = float(kipping_like(BIG, rr, tLp, tobs - tLp, T) / kipping_like(SMALL, rr, tLp, tobs - tLp, T)) num_exact = float(like2_density(BIG, rr, tLp, tobs, T) / like2_density(SMALL, rr, tLp, tobs, T)) val.append(dict(rate=rr, exact_closed=float(cf_exact), exact_numeric=num_exact, exact_general_routine=float(plateau_ratio_general([rr], tLp, tobs, T)), interval_closed=float(cf_int), interval_numeric=num_int)) out["validation"] = val out["curve_baseline"] = {"rate": [float(v) for v in r[::40]], "exact": [float(v) for v in plateau_ratio_one_later_step(r, tLp, T)[0][::40]], "interval": [float(v) for v in plateau_ratio_one_later_step(r, tLp, T)[1][::40]]} return out # --------------------------------------------------------------------------------- def like2_density(lL, lI, tLp, tobs, T): """Density-form 2-step likelihood: abiogenesis in [0, tLp], observers exactly at tobs.""" from timing_lib import log_exprel lL, lI = np.broadcast_arrays(lL, lI) num = (np.log(lL) + np.log(lI) - lI * tobs + np.log(tLp) + log_exprel((lI - lL) * tLp)) return np.exp(num - log_hypoexp_cdf2(lL, lI, T)) def prior_sensitivity(): from scipy import stats t0, F = cfg.T0_ZIRCON, 0.9 T, tobs = t0 + F, t0 gL = np.arange(-100, 3 + 1e-9, 0.02) # log10 lambda_L grid gI = np.arange(-6, 2 + 1e-9, 0.05) # log10 lambda_I grid (log-uniform prior, marginalised) lL = 10.0 ** gL lI = 10.0 ** gI rows = [] for life_lab in ["LUCA 4.2 Ga (Moody+ 2024)", "Dresser 3.48 Ga (accepted)"]: tLp = t0 - cfg.LIFE_DATES[life_lab] like = like2_density(lL[:, None], lI[None, :], tLp, tobs, T) # (nL, nI) like_margI = like.mean(axis=1) # log-uniform prior on lambda_I like_hardI = like2_density(lL, 1e-9, tLp, tobs, T) # lambda_I -> 0 limit priors = {} for floor in (-3, -10, -30, -100): p = np.where(gL >= floor - 1e-9, 1.0, 0.0) priors[f"log-uniform lambda_L in [1e{floor}, 1e3] /Gyr"] = (p / p.sum(), 0.0, 0.0) # log-normal in log10(lambda_L T), sigma = 50 dex (Sandberg-Drexler-Ord style); tails beyond the # grid sit on the two likelihood plateaus and are added analytically mu, sig = -np.log10(T), 50.0 edges = np.concatenate([[gL[0] - 0.01], gL + 0.01]) cdf = stats.norm.cdf(edges, mu, sig) priors["log-normal, median 1 event per window, sigma = 50 dex"] = (np.diff(cdf), cdf[0], 1 - cdf[-1]) lam_pdf = T / (np.pi * np.sqrt(np.expm1(np.minimum(lL * T, 700)))) # Jeffreys in f p = lam_pdf * lL priors["Jeffreys in f (Kipping 2020)"] = (p / p.sum(), 0.0, 0.0) p = T * np.exp(-lL * T) * lL priors["uniform in f"] = (p / p.sum(), 0.0, 0.0) fast = gL > 0.0 # lambda_L > 1 per Gyr for plab, (p, lo_tail, hi_tail) in priors.items(): for ilab, lk in [("lambda_I marginalised (log-uniform 1e-6..1e2)", like_margI), ("lambda_I -> 0 (hard)", like_hardI)]: norm = p.sum() + lo_tail + hi_tail prior_fast = (p[fast].sum() + hi_tail) / norm z = (p * lk).sum() + lo_tail * lk[0] + hi_tail * lk[-1] post_fast = ((p * lk)[fast].sum() + hi_tail * lk[-1]) / z f_life = -np.expm1(-lL * T) post_f50 = ((p * lk)[f_life > 0.5].sum() + hi_tail * lk[-1]) / z rows.append(dict(life_evidence=life_lab, prior=plab, intelligence_rate=ilab, prior_P_fast=prior_fast, post_P_fast=post_fast, implied_BF=(post_fast / (1 - post_fast)) / (prior_fast / (1 - prior_fast)), post_P_flife_gt_0p5=post_f50)) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a1_prior_sensitivity.csv"), index=False) return df # --------------------------------------------------------------------------------- def moody_bf(rng): d = pd.read_csv(os.path.join(cfg.DATA, "moody_node_ages.csv")) luca = d[d.node == "LUCA"].copy() rows = [] for _, r in luca.iterrows(): row = dict(analysis=r.analysis, clock=r.clock, stage=r.stage, copy=r["copy"], LUCA_mean_Ga=r.mean_Ma / 1e3, LUCA_q025_Ga=r.q025_Ma / 1e3, LUCA_q975_Ga=r.q975_Ma / 1e3) for lab, t0 in [("t0=4.404", cfg.T0_ZIRCON), ("t0=4.51", cfg.T0_MOON)]: tl = t0 - r.mean_Ma / 1e3 row[f"tLp_{lab}"] = tl row[f"BF_kipping_{lab}"] = (t0 + 0.9) / (2 * tl) if tl > 0 else np.inf rows.append(row) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(cfg.DATA, "a1_moody_luca_bf.csv"), index=False) post = tab[tab.stage == "posterior"] b44, b451 = post["BF_kipping_t0=4.404"], post["BF_kipping_t0=4.51"] worst = post.loc[post["tLp_t0=4.404"].idxmin()] summ = { "n_posterior_estimates": int(len(post)), "LUCA_mean_range_Ga": [float(post.LUCA_mean_Ga.min()), float(post.LUCA_mean_Ga.max())], # Point values T/(2 tLp) diverge as a LUCA estimate approaches the assumed start date, so the # raw range across variants is NOT a meaningful summary (referee report). We keep the # median, the quartiles and the share above 10; the date-marginalised factors below are # the numbers to quote. "BF_t0_4.404_median": float(b44.median()), "BF_t0_4.404_quartiles": [float(b44.quantile(0.25)), float(b44.quantile(0.75))], "BF_t0_4.404_min": float(b44.min()), "BF_t0_4.51_median": float(b451.median()), "BF_t0_4.51_quartiles": [float(b451.quantile(0.25)), float(b451.quantile(0.75))], "BF_t0_4.51_min": float(b451.min()), "BF_t0_4.51_max": float(b451.max()), "frac_BF_gt_10_t0_4.404": float((b44 > 10).mean()), "frac_BF_gt_10_t0_4.51": float((b451 > 10).mean()), "n_variants_whose_95pct_interval_contains_4.404": int((post.LUCA_q975_Ga >= cfg.T0_ZIRCON).sum()), "divergent_point_value_t0_4.404": { "analysis": str(worst.analysis), "clock": str(worst.clock), "copy": int(worst["copy"]), "LUCA_mean_Ga": float(worst.LUCA_mean_Ga), "LUCA_95pct_Ga": [float(worst.LUCA_q025_Ga), float(worst.LUCA_q975_Ga)], "gap_to_start_Myr": float(1e3 * worst["tLp_t0=4.404"]), "point_BF": float(worst["BF_kipping_t0=4.404"]), "note": "division by a 2.8 Myr gap; an artefact, not a result"}, } # calibration-only (no sequence data) priors on the LUCA node, by tree configuration, and the root pri = tab[tab.stage == "prior"] def _grp(g): return dict(n=int(len(g)), mean_Ga=float(g.LUCA_mean_Ga.mean()), q025_Ga=float(g.LUCA_q025_Ga.mean()), q975_Ga=float(g.LUCA_q975_Ga.mean())) is_nb = pri.analysis.str.contains("not cross-braced") is_fb = pri.analysis.str.contains("fossil-node cross-braced") summ["calibration_only_prior_LUCA"] = {"cross-braced (A)": _grp(pri[~is_nb & ~is_fb]), "fossil-node cross-braced (B)": _grp(pri[is_fb]), "not cross-braced": _grp(pri[is_nb])} root = d[d.node == "ROOT"] rm = root[root.analysis.str.contains("paper main")] summ["root_paper_main"] = [dict(clock=r.clock, stage=r.stage, mean_Ga=r.mean_Ma / 1e3, q025_Ga=r.q025_Ma / 1e3, q975_Ga=r.q975_Ma / 1e3) for r in rm.itertuples()] summ["root_hard_maximum_Ga"] = 4.520 main = tab[tab.analysis.str.contains("paper main")] summ["paper_main"] = main[["clock", "stage", "LUCA_mean_Ga", "LUCA_q025_Ga", "LUCA_q975_Ga", "BF_kipping_t0=4.404", "BF_kipping_t0=4.51"]].to_dict(orient="records") # the two independently dated copies of the LUCA node when cross-bracing is switched off nb = post[post.analysis.str.contains("not cross-braced")] gaps = [] for (an, ck), g in nb.groupby(["analysis", "clock"]): if len(g) == 2: gaps.append(dict(analysis=an, clock=ck, copy1_Ga=float(g.LUCA_mean_Ga.max()), copy2_Ga=float(g.LUCA_mean_Ga.min()), gap_Gyr=float(g.LUCA_mean_Ga.max() - g.LUCA_mean_Ga.min()))) summ["non_cross_braced_copy_gaps"] = gaps # ---- proper marginalisation over date uncertainty -------------------------------------- M = 400_000 res = {} def luca_draws(kind): if kind == "paper main (ILN+GBM)": mm = main[main.stage == "posterior"] else: # every posterior estimate in the archive, equally weighted (clock + method systematics) mm = post idx = rng.integers(0, len(mm), M) c = mm.LUCA_mean_Ga.values[idx] lo = mm.LUCA_q025_Ga.values[idx] hi = mm.LUCA_q975_Ga.values[idx] z = rng.standard_normal(M) return c + np.where(z < 0, (c - lo) / 1.959964, (hi - c) / 1.959964) * z for kind in ["paper main (ILN+GBM)", "all 38 posterior estimates"]: for t0kind in ["t0 = 4.404 fixed", "t0 ~ U(4.404, 4.50)", "t0 = 4.51 fixed"]: for fut in ["0.9 fixed", "short", "long", "mix"]: luca_age = luca_draws(kind) if t0kind == "t0 = 4.404 fixed": t0 = np.full(M, cfg.T0_ZIRCON) elif t0kind == "t0 = 4.51 fixed": t0 = np.full(M, cfg.T0_MOON) else: t0 = cfg.draw_t0(rng, M) F = np.full(M, 0.9) if fut == "0.9 fixed" else cfg.draw_future(rng, M, fut) ok = luca_age < t0 - 0.005 # LUCA cannot predate habitability tl = (t0 - luca_age)[ok] tobs = t0[ok] T = (t0 + F)[ok] entry = {"frac_draws_LUCA_older_than_t0": float(1 - ok.mean()), "median_tLp": float(np.median(tl)), "mean_tLp": float(tl.mean())} for m in (1, 2, 3, 4, 5): l_inf = m * tobs ** (m - 1) / T ** m l_0 = (m + 1) * (tobs ** m - (tobs - tl) ** m) / T ** (m + 1) entry[f"B_marg_m{m}"] = float(l_inf.mean() / l_0.mean()) if m == 1: b = l_inf / l_0 entry["B_pointwise_median_m1"] = float(np.median(b)) entry["B_pointwise_16_84_m1"] = [float(q) for q in np.percentile(b, [15.87, 84.13])] # absolute-time pacing with deadline at the GOE c = t0[ok] - cfg.GOE[0] entry["B_marg_pacing_GOE_deadline"] = float(1.0 / (tl / c).mean()) res[f"{kind} | {t0kind} | future {fut}"] = entry summ["marginalised"] = res return tab, summ def run(rng): out = {} out["kipping2020"] = reproduce_2020(rng) out["kipping2025"] = reproduce_2025(rng) bt, val = bmax_tables() out["bmax_validation"] = val base = bt[(bt.life_evidence.str.startswith("LUCA")) & (bt.future.str.startswith("0.9"))] out["bmax_LUCA_future0.9"] = {r.structure: round(r.B_max, 2) for r in base.itertuples()} base = bt[(bt.life_evidence.str.startswith("LUCA")) & (bt.future.str.startswith("1.73"))] out["bmax_LUCA_future1.73"] = {r.structure: round(r.B_max, 2) for r in base.itertuples()} base = bt[(bt.life_evidence.str.startswith("Dresser")) & (bt.future.str.startswith("0.9"))] out["bmax_Dresser_future0.9"] = {r.structure: round(r.B_max, 2) for r in base.itertuples()} ps = prior_sensitivity() out["prior_sensitivity"] = ps.to_dict(orient="records") out["later_process_sensitivity"] = later_process_sensitivity() tab, summ = moody_bf(rng) out["moody"] = summ # the intelligence-rate likelihood: monotone decreasing iff t_obs > T/2 out["intelligence_rate_mle_is_zero_if_tobs_gt_T_over_2"] = { "tobs": 4.404, "T_needed_for_interior_maximum": 2 * 4.404, "remaining_time_needed_Gyr": 4.404} return out if __name__ == "__main__": import json res = run(np.random.default_rng(cfg.SEED)) print(json.dumps(res, indent=1, default=float)[:6000])