""" Step 4. Itemised search completeness for a second origin of life, and the likelihood of the null result to date. Input : ../data/completeness_factors.csv (STRUCTURED EXPERT JUDGEMENT, itemised: every row is my judgement with a low / most-likely / high value and a stated basis; no row is a measurement; change a row and rerun to see what it does) Output: ../data/completeness_by_habitat.csv, ../data/completeness_pathways.csv, ../data/null_likelihood_curves.csv, ../data/completeness_summary.json, ../data/completeness_draws.npz Definitions (lambda = expected independent origins of life per 5.6-Gyr habitable window of an Earth-like planet, the model's abiogenesis parameter): q_h probability that habitat h has been habitable long enough for an origin o_h its opportunity for an origin, in units of that window (duration x suitability relative to Earth) c_h search completeness: probability that the evidence in hand would already have revealed (and the community accepted) an independent origin in h, given that ONE happened there. Built from pathways: c_h = 1 - prod_groups [ 1 - s_g * d_g ], d_g = 1 - prod_pathways (1 - success * f * p) s = the life persists, or leaves preservable remains, where this group of searches looks f = share of the places it would be that has actually been sampled p = the instruments used would have recognised it and the result been accepted For Earth: c = persist * c_alive, c_alive = sum over four exclusive classes (sequenceable or not) x (abundant or rare). The likelihood of the null. How the chance of detection grows with the NUMBER of origins in a habitat is a modelling choice that the factor table does not settle, and it matters in the model's fast camp (lambda 1 to 1000), which is where the whole update comes from. Three structures are computed side by side (mu = o_h * lambda, expected origins in h): 'habitat' fates fully correlated within a habitat: every factor acts once per habitat, however many origins L_h = 1 - q c (1 - exp(-mu)) saturates at 1 - q c 'origin_persist' each origin persists / is preserved independently (Poisson thinning by s_g, or by `persist` on Earth); sampling and recognition (f, p; c_alive on Earth) still act once per habitat L_h = 1 - q + q * E_N[ prod_g (1 - d_g (1 - (1 - s_g)^N)) ] (exact, N ~ Poisson(mu)) 'origin_all' every origin is an independent chance of detection for every factor: the functional form of the model's own likelihood, exp(-sg_frac * lambda) L_h = 1 - q + q * exp(-mu c) All three agree for lambda << 1, where L ~ exp(-lambda F) with F = sum_h q_h c_h o_h (the small-lambda exposure). F is NOT comparable with the model's sg_frac = 0.3 in the fast camp, because the model's form never saturates; the like-for-like comparison is the detection chance D_fast = 1 - over the fast camp (log10 lambda uniform on 0..3). Usage: python3 completeness.py [n_draws=20000] """ import csv import json import os import sys import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "..", "data") EPOCHS = ["today", "funded", "unfunded"] # cumulative: each includes the ones before it EPOCH_LABEL = {"today": "evidence in hand, October 2026", "funded": "plus funded missions (results by about 2038)", "unfunded": "plus two missions not funded (Mars cores returned; Enceladus lander; 2040s-2050s)"} HABITATS = ["Earth", "Mars", "Europa", "Enceladus", "Titan", "Venus", "Small bodies"] LAM_GRID = 10 ** np.linspace(-8, 3, 1101) # 0.01 dex FAST = (LAM_GRID >= 1.0) # the model's fast camp: log10 lambda ~ U(0, 3) STRUCTURES = ["habitat", "origin_persist", "origin_all"] STRUCT_LABEL = {"habitat": "per habitat (fates fully correlated)", "origin_persist": "per origin for persistence, per habitat for sampling and recognition", "origin_all": "per origin for every factor (the model's own functional form)"} def draw(rng, lo, mode, hi, dist, n): if dist == "const" or hi <= lo: return np.full(n, mode) if dist == "tri": return rng.triangular(lo, mode, hi, n) if dist == "logtri": return 10 ** rng.triangular(np.log10(lo), np.log10(mode), np.log10(hi), n) raise ValueError(dist) def load_factors(): rows = list(csv.DictReader(open(os.path.join(DATA, "completeness_factors.csv")))) for r in rows: for k in ("low", "mode", "high"): r[k] = float(r[k]) r["year"] = int(r["year"]) return rows def sample_all(rows, n, seed=20261012, point=None): """One draw of every factor. point = 'low' | 'mode' | 'high' returns that column instead (n ignored -> length 1).""" rng = np.random.default_rng(seed) vals = {} for i, r in enumerate(rows): key = (r["epoch"], r["habitat"], r["group"], r["pathway"], r["factor"]) if point: vals[key] = np.array([r[point]]) else: vals[key] = draw(rng, r["low"], r["mode"], r["high"], r["dist"], n) return vals def habitat_terms(rows, vals, upto): """q, o, c per habitat for the cumulative epoch `upto`; also per-pathway detection probabilities.""" use = EPOCHS[:EPOCHS.index(upto) + 1] n = len(next(iter(vals.values()))) out, paths = {}, [] for h in HABITATS: g = lambda f: vals[("today", h, "-", "-", f)] q, o = g("q"), g("o") if h == "Earth": pna, pab = g("p_na"), g("p_abund") c_today = pna * pab * g("c_na_abund") + pna * (1 - pab) * g("c_na_rare") + (1 - pna) * pab * g("c_x_abund") + (1 - pna) * (1 - pab) * g("c_x_rare") c = g("persist") * c_today out[h] = dict(q=q, o=o, c=c, c_alive_today=c_today, persist=g("persist"), groups=[(g("persist"), c_today)]) continue groups = sorted({r["group"] for r in rows if r["habitat"] == h and r["group"] != "-"}) miss_h = np.ones(n) glist = [] for grp in groups: s = vals[("today", h, grp, "*", "s")] pw = sorted({(r["epoch"], r["pathway"], r["year"]) for r in rows if r["habitat"] == h and r["group"] == grp and r["pathway"] != "*" and r["epoch"] in use}, key=lambda t: (EPOCHS.index(t[0]), t[1])) miss_g = np.ones(n) for ep, pth, yr in pw: a = vals.get((ep, h, grp, pth, "success"), np.ones(n)) d = a * vals[(ep, h, grp, pth, "f")] * vals[(ep, h, grp, pth, "p")] miss_g *= 1 - d paths.append(dict(habitat=h, group=grp, pathway=pth, epoch=ep, year=yr, c_path=s * d, s=s, success=a, f=vals[(ep, h, grp, pth, "f")], p=vals[(ep, h, grp, pth, "p")])) miss_h *= 1 - s * (1 - miss_g) glist.append((s, 1 - miss_g)) out[h] = dict(q=q, o=o, c=1 - miss_h, groups=glist) return out, paths def null_likelihood(terms, lam=LAM_GRID, only=None, structure="habitat", ceiling=None, kc=1.0, ko=1.0): """L(lambda) for every draw: array (n_draws, n_lambda). structure: 'habitat' | 'origin_persist' | 'origin_all' (see the module docstring). ceiling: None the searches as itemised 'surviving' everything that survives or is preserved today is found (d_g = 1; limited by s_g / persist) 'ever' any life that ever arose is found (c = 1 everywhere, Earth included); the same in every structure kc, ko: common-mode multipliers on every detection factor (capped at 1) and on every opportunity.""" L = None for h, t in terms.items(): if only and h not in only: continue q = t["q"][:, None] mu = ko * t["o"][:, None] * lam[None, :] groups = [(s, np.ones_like(d) if ceiling == "surviving" else np.minimum(1.0, kc * d)) for s, d in t["groups"]] if ceiling == "ever": term = 1 - q * -np.expm1(-mu) elif structure == "habitat": c = 1 - np.prod([1 - s * d for s, d in groups], axis=0) if ceiling is None and kc != 1.0: c = np.minimum(1.0, kc * t["c"]) # the first-pass definition of the common-mode scaling term = 1 - (t["q"] * c)[:, None] * -np.expm1(-mu) elif structure == "origin_all": c = 1 - np.prod([1 - s * d for s, d in groups], axis=0) term = 1 - q * -np.expm1(-mu * c[:, None]) elif structure == "origin_persist": # exact expectation over N ~ Poisson(mu) of prod_g [1 - d_g + d_g (1 - s_g)^N], expanded over subsets of groups inner = np.zeros_like(mu) ng = len(groups) for mask in range(1 << ng): coef = np.ones_like(t["q"]); keep = np.ones_like(t["q"]) for i, (sg, dg) in enumerate(groups): if mask >> i & 1: coef = coef * dg; keep = keep * (1 - sg) else: coef = coef * (1 - dg) inner += coef[:, None] * np.exp(-mu * (1 - keep)[:, None]) term = 1 - q + q * inner else: raise ValueError(structure) L = term if L is None else L * term return L def qs(x, q=(0.05, 0.5, 0.95)): return [float(v) for v in np.quantile(x, q)] def raup_valentine(k=10): """Per-origin survival probability p (independent clades) below which 'only one of k survives, given any does' is more likely than not: the condition Raup & Valentine (1983) report for 10 origins.""" p = np.linspace(1e-4, 0.9, 9000) one = k * p * (1 - p) ** (k - 1) / (1 - (1 - p) ** k) return float(p[np.argmin(np.abs(one - 0.5))]) def main(): n = int(sys.argv[1]) if len(sys.argv) > 1 else 20000 rows = load_factors() vals = sample_all(rows, n) summary = dict(n_draws=n, n_factor_rows=len(rows), epochs=EPOCH_LABEL, structures=STRUCT_LABEL, raup_valentine_p_max_for_10_origins=raup_valentine(10), habitats={}, totals={}) hab_rows, path_rows, curves = [], [], {"lambda": LAM_GRID} save = {} for ep in EPOCHS: terms, paths = habitat_terms(rows, vals, ep) F = sum(t["q"] * t["c"] * t["o"] for t in terms.values()) # small-lambda exposure (NOT comparable with sg_frac in the fast camp) save[f"F_{ep}"] = F tot = dict(F_exposure_q=qs(F), F_exposure_mean=float(F.mean())) for st in STRUCTURES: L = null_likelihood(terms, structure=st) Lbar = L.mean(0) # the likelihood with the factor uncertainty integrated out sfx = "" if st == "habitat" else "_" + st curves[f"L_mean_{ep}{sfx}"] = Lbar if st != "origin_all": curves[f"L_p05_{ep}{sfx}"], curves[f"L_p50_{ep}{sfx}"], curves[f"L_p95_{ep}{sfx}"] = np.quantile(L, [0.05, 0.5, 0.95], axis=0) D_fast = 1 - L[:, FAST].mean(1) # P(the evidence would have shown a second origin | life starts easily), per draw save[f"Dfast_{ep}{sfx}"] = D_fast tot[f"D_fast_q{sfx}"] = qs(D_fast); tot[f"D_fast_mean{sfx}"] = float(D_fast.mean()) if st != "habitat": continue D_sat = 1 - np.prod([1 - t["q"] * t["c"] for t in terms.values()], axis=0) # the same if every habitable habitat had an origin tot.update(D_fast_from_mean_likelihood=float(1 - Lbar[FAST].mean()), D_saturated_q=qs(D_sat), D_saturated_mean=float(D_sat.mean())) del L # the model's own form taken literally: exp(-F lambda), F drawn with the factors, and with F fixed at its median Lm = np.exp(-F[:, None] * LAM_GRID[None, :]) curves[f"L_mean_{ep}_model_form"] = Lm.mean(0) curves[f"L_{ep}_model_form_F_median"] = np.exp(-np.median(F) * LAM_GRID) tot["D_fast_model_form_q"] = qs(1 - Lm[:, FAST].mean(1)); tot["D_fast_model_form_mean"] = float(1 - Lm[:, FAST].mean()) tot["D_fast_model_form_F_median"] = float(1 - np.exp(-np.median(F) * LAM_GRID[FAST]).mean()) del Lm summary["totals"][ep] = tot for h, t in terms.items(): dfh = 1 - null_likelihood({h: t})[:, FAST].mean(1) dfo = 1 - null_likelihood({h: t}, structure="origin_persist")[:, FAST].mean(1) row = dict(epoch=ep, habitat=h, q_p05=qs(t["q"])[0], q_p50=qs(t["q"])[1], q_p95=qs(t["q"])[2], o_p05=qs(t["o"])[0], o_p50=qs(t["o"])[1], o_p95=qs(t["o"])[2], c_p05=qs(t["c"])[0], c_p50=qs(t["c"])[1], c_p95=qs(t["c"])[2], c_mean=float(t["c"].mean()), exposure_qco_p05=qs(t["q"] * t["c"] * t["o"])[0], exposure_qco_p50=qs(t["q"] * t["c"] * t["o"])[1], exposure_qco_p95=qs(t["q"] * t["c"] * t["o"])[2], exposure_qco_mean=float((t["q"] * t["c"] * t["o"]).mean()), D_fast_p05=qs(dfh)[0], D_fast_p50=qs(dfh)[1], D_fast_p95=qs(dfh)[2], D_fast_mean=float(dfh.mean()), D_fast_origin_persist_p05=qs(dfo)[0], D_fast_origin_persist_p50=qs(dfo)[1], D_fast_origin_persist_p95=qs(dfo)[2], D_fast_origin_persist_mean=float(dfo.mean())) if h == "Earth": row.update(c_alive_today_p05=qs(t["c_alive_today"])[0], c_alive_today_p50=qs(t["c_alive_today"])[1], c_alive_today_p95=qs(t["c_alive_today"])[2], c_alive_today_mean=float(t["c_alive_today"].mean())) hab_rows.append(row) save[f"c_{ep}_{h.replace(' ', '_')}"] = t["c"] save[f"q_{ep}_{h.replace(' ', '_')}"] = t["q"] save[f"o_{ep}_{h.replace(' ', '_')}"] = t["o"] if h == "Earth": save[f"calive_{ep}_Earth"] = t["c_alive_today"] save[f"persist_{ep}_Earth"] = t["persist"] if ep == "unfunded": # pathways once, from the widest epoch for p in paths: path_rows.append(dict(habitat=p["habitat"], group=p["group"], pathway=p["pathway"], first_epoch=p["epoch"], year=p["year"], c_path_p05=qs(p["c_path"])[0], c_path_p50=qs(p["c_path"])[1], c_path_p95=qs(p["c_path"])[2], c_path_mean=float(p["c_path"].mean()), s_p50=qs(p["s"])[1], success_p50=qs(p["success"])[1], f_p50=qs(p["f"])[1], p_p50=qs(p["p"])[1])) if ep == "today": # Two ceilings, each defined the same way for every body (the first pass mixed them: Earth capped by persistence, other bodies not): # 'ever' any life that ever arose on any body is found (c = 1 everywhere) # 'surviving' everything that still lives or is preserved where a search could look is found (limited by persist and s) Le = null_likelihood(terms, ceiling="ever") curves["L_mean_ceiling_ever"] = Le.mean(0) Fe = sum(t["q"] * t["o"] for t in terms.values()) summary["totals"]["ceiling_ever"] = dict(F_exposure_q=qs(Fe), F_exposure_mean=float(Fe.mean()), D_fast_q=qs(1 - Le[:, FAST].mean(1)), D_fast_mean=float(1 - Le[:, FAST].mean())) del Le Fsv = sum(t["q"] * (1 - np.prod([1 - s for s, d in t["groups"]], axis=0)) * t["o"] for t in terms.values()) summary["totals"]["ceiling_surviving"] = dict(F_exposure_q=qs(Fsv), F_exposure_mean=float(Fsv.mean())) for st in ("habitat", "origin_persist"): Lv = null_likelihood(terms, ceiling="surviving", structure=st) sfx = "" if st == "habitat" else "_" + st curves[f"L_mean_ceiling_surviving{sfx}"] = Lv.mean(0) summary["totals"]["ceiling_surviving"][f"D_fast_q{sfx}"] = qs(1 - Lv[:, FAST].mean(1)) summary["totals"]["ceiling_surviving"][f"D_fast_mean{sfx}"] = float(1 - Lv[:, FAST].mean()) del Lv # Solar System only (no Earth term), for comparison with the model's sg_frac, which is about other bodies Ls = null_likelihood(terms, only=[h for h in HABITATS if h != "Earth"]) curves["L_mean_today_no_earth"] = Ls.mean(0) Fs = sum(t["q"] * t["c"] * t["o"] for h, t in terms.items() if h != "Earth") summary["totals"]["today_solar_system_only"] = dict(F_exposure_q=qs(Fs), F_exposure_mean=float(Fs.mean()), D_fast_q=qs(1 - Ls[:, FAST].mean(1)), D_fast_mean=float(1 - Ls[:, FAST].mean())) del Ls Fo = sum(t["q"] * t["o"] for h, t in terms.items() if h != "Earth") summary["totals"]["opportunity_sum_q_o_solar_system"] = dict(q=qs(Fo), mean=float(Fo.mean())) # low / mode / high corner cases (all factors at one end together: a deliberate over-statement of the range) for pt in ("low", "mode", "high"): v1 = sample_all(rows, 1, point=pt) for ep in EPOCHS: terms, _ = habitat_terms(rows, v1, ep) L = null_likelihood(terms)[0] curves[f"L_all_{pt}_{ep}"] = L Lo = null_likelihood(terms, structure="origin_persist")[0] curves[f"L_all_{pt}_{ep}_origin_persist"] = Lo summary["totals"].setdefault("corner_" + pt, {})[ep] = dict(F_exposure=float(sum(t["q"][0] * t["c"][0] * t["o"][0] for t in terms.values())), D_fast=float(1 - L[FAST].mean()), D_fast_origin_persist=float(1 - Lo[FAST].mean()), c={h: float(t["c"][0]) for h, t in terms.items()}) curves["L_model_thorough"] = np.exp(-0.3 * LAM_GRID) # the model's own 'thorough search finds none' (sg_frac = 0.3), without the constant (1 - sg_fp) summary["totals"]["model_thorough"] = dict(sg_frac=0.3, D_fast=float(1 - np.exp(-0.3 * LAM_GRID[FAST]).mean())) with open(os.path.join(DATA, "completeness_by_habitat.csv"), "w", newline="") as f: keys = list(dict.fromkeys(k for r in hab_rows for k in r)) w = csv.DictWriter(f, fieldnames=keys); w.writeheader(); w.writerows(hab_rows) with open(os.path.join(DATA, "completeness_pathways.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(path_rows[0])); w.writeheader(); w.writerows(path_rows) with open(os.path.join(DATA, "null_likelihood_curves.csv"), "w", newline="") as f: keys = list(curves) w = csv.writer(f); w.writerow(keys) for i in range(len(LAM_GRID)): w.writerow([f"{curves[k][i]:.8g}" for k in keys]) np.savez_compressed(os.path.join(DATA, "completeness_draws.npz"), **{k: v.astype(np.float32) for k, v in save.items()}) json.dump(summary, open(os.path.join(DATA, "completeness_summary.json"), "w"), indent=1) print("Raup & Valentine check: per-origin survival below", round(summary["raup_valentine_p_max_for_10_origins"], 3), "makes a single surviving tree the likelier outcome of 10 origins") for ep in EPOCHS: t = summary["totals"][ep] print(f"\n== {ep}: {EPOCH_LABEL[ep]}") print(f" small-lambda exposure F = sum q c o: median {t['F_exposure_q'][1]:.4f} [{t['F_exposure_q'][0]:.4f}, {t['F_exposure_q'][2]:.4f}] mean {t['F_exposure_mean']:.4f}") print(" D_fast = chance the evidence shows a second origin if life starts easily, by likelihood structure (median [5%, 95%] mean over factor draws):") for st in STRUCTURES: sfx = "" if st == "habitat" else "_" + st q_ = t[f"D_fast_q{sfx}"] print(f" {st:15s} {q_[1]:.3f} [{q_[0]:.3f}, {q_[2]:.3f}] mean {t[f'D_fast_mean{sfx}']:.3f}") q_ = t["D_fast_model_form_q"] print(f" {'exp(-F lambda)':15s} {q_[1]:.3f} [{q_[0]:.3f}, {q_[2]:.3f}] mean {t['D_fast_model_form_mean']:.3f} (F fixed at its median: {t['D_fast_model_form_F_median']:.3f})") for r in hab_rows: if r["epoch"] == ep: extra = f" c(alive today) {r['c_alive_today_p50']:.3f} [{r['c_alive_today_p05']:.3f}, {r['c_alive_today_p95']:.3f}]" if r["habitat"] == "Earth" else "" print(f" {r['habitat']:13s} q {r['q_p50']:.2f} o {r['o_p50']:.4f} c {r['c_p50']:.4f} [{r['c_p05']:.4f}, {r['c_p95']:.4f}] q*c*o {r['exposure_qco_p50']:.5f} D_fast {r['D_fast_p50']:.4f} (per-origin persistence {r['D_fast_origin_persist_p50']:.4f}){extra}") for k in ("ceiling_ever", "ceiling_surviving", "today_solar_system_only", "model_thorough"): print(k, summary["totals"][k]) print("sum of q*o over the Solar System habitats (no Earth):", summary["totals"]["opportunity_sum_q_o_solar_system"]) print("corners:", {pt: {ep: (round(v["F_exposure"], 5), round(v["D_fast"], 4), round(v["D_fast_origin_persist"], 4)) for ep, v in summary["totals"]["corner_" + pt].items()} for pt in ("low", "mode", "high")}) if __name__ == "__main__": main()