""" Step 5. What the search record to date is worth: re-weight the model's posterior samples by the likelihood of the record. Inputs : ../data/model_samples.npz (+ meta), ../data/completeness_factors.csv (through completeness.py) Outputs: ../data/update_scenarios.csv, ../data/update_camps.csv, ../data/update_range_over_completeness.csv, ../data/ideal_search_map.csv, ../data/detection_forecast.csv, ../data/claim_value.csv, ../data/cheyava_sensitivity.csv, ../data/cheyava_breakeven.csv, ../data/independence.csv, ../data/transfer_literature.csv, ../data/update_summary.json The model's own likelihood for 'thorough search finds none' is (1 - sg_fp) * exp(-sg_frac * lambda), sg_frac = 0.3. The itemised likelihood used here is a product over habitats, and how it grows with the number of origins in a habitat is a modelling choice (see completeness.py). Every result is therefore given for three structures: habitat L = prod_h [1 - q c (1 - exp(-o lambda))] (fates fully correlated within a habitat) origin_persist persistence / preservation thins origins independently; sampling and recognition act once per habitat origin_all L = prod_h [1 - q + q exp(-o lambda c)], and its literal version exp(-F lambda), F = sum q c o (the model's own functional form with sg_frac = F) The constant (1 - sg_fp) cancels in every posterior. The record is not a pure null: the Cheyava Falls 'potential biosignature' (Hurowitz et al. 2025) is scored separately with a likelihood ratio LR that nobody can pin down today; the break-even LR at which it cancels the whole null update is reported. """ import csv import json import os import numpy as np import completeness as C HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "..", "data") SCHOOLS = ["site", "ssa_all", "ssa_now", "ssa_obs"] CLASSES = ["empty_universe", "rare", "en_route", "stuck", "unnoticed", "hidden", "excluded"] LG = np.log10(C.LAM_GRID) N_DRAW = 4000 ST_MAIN = ["habitat", "origin_persist"] # the two structures carried through every table ST_NAME = {"habitat": "per habitat", "origin_persist": "per origin (persistence)", "origin_all": "per origin (every factor)", "model_form": "model's own form, exp(-F lambda)"} EP_NAME = {"today": "Null to date", "funded": "Null after the funded missions (about 2038)", "unfunded": "Null after two unfunded missions as well (2050s)"} CHEYAVA_LR = (1.0, 1.1, 1.2, 1.5, 2.0, 3.0, 10.0) # likelihood ratio P(features | life in that mud) / P(features | no life there) CHEYAVA_SB_FIXED = (0.1, 0.2, 0.4, 0.8) # fixed values of sB shown beside the range in the factor table # Worth, Sigurdsson & House 2013, Astrobiology 13, 1155, published Tables 4 and 5 (PMC3870607): 3-m rocks ejected below # sterilising shock that reach the target in 3.5 Gyr. The arXiv v1 preprint has different outer-moon values and is not used. TRANSFER = [("Earth -> Mars", 3.6e5, 4e4), ("Mars -> Earth", 2.1e7, 5e5), ("Earth -> Venus", 2.6e7, 3e5), ("Mars -> Venus", 1.2e7, 4e5), ("Earth -> Europa", 6, 0.9), ("Mars -> Europa", 20, 4), ("Earth -> Titan", 1, 0.9), ("Earth -> Enceladus", 1, 0.7)] def load(): S = dict(np.load(os.path.join(DATA, "model_samples.npz"))) meta = json.load(open(os.path.join(DATA, "model_samples_meta.json"))) S["lam"] = 10.0 ** S["log10_lam"].astype(np.float64) return S, meta def interp_L(curve, log10_lam, left=1.0): """Likelihood at each sample; below the grid (lambda < 1e-8) the null says nothing: L = 1.""" return np.interp(log10_lam, LG, curve, left=left, right=curve[-1]) def stats(S, w): ws = w.sum() out = dict(p_fast=float(w[S["fast"]].sum() / ws), p_alone=float(w[S["alone"]].sum() / ws)) for i, k in enumerate(CLASSES): out[k] = float(w[S["cls"] == i].sum() / ws) out["nobody"] = out["empty_universe"] + out["rare"] out["spread"] = out["en_route"] + out["unnoticed"] + out["hidden"] + out["excluded"] out["ess"] = float(ws ** 2 / (w ** 2).sum()) return out def per_seed(S, w0, like, key): """Value of `key` and its change from the base, seed by seed (paired), for the Monte Carlo error.""" vals, diffs = [], [] for sd in np.unique(S["seed"]): m = S["seed"] == sd Sm = {k: v[m] for k, v in S.items()} a = stats(Sm, w0[m])[key] b = stats(Sm, w0[m] * like[m])[key] vals.append(b); diffs.append(b - a) return float(np.std(vals, ddof=1) / np.sqrt(len(vals))), float(np.mean(diffs)), float(np.std(diffs, ddof=1) / np.sqrt(len(diffs))) def bisect(fn, lo, hi, n=50): """Root of a decreasing function fn on [lo, hi]; returns nan if there is no sign change.""" if fn(lo) < 0 or fn(hi) > 0: return float("nan") for _ in range(n): mid = 0.5 * (lo + hi) if fn(mid) > 0: lo = mid else: hi = mid return 0.5 * (lo + hi) def main(): S, meta = load() rows = C.load_factors() vals = C.sample_all(rows, N_DRAW) terms = {ep: C.habitat_terms(rows, vals, ep)[0] for ep in C.EPOCHS} Ldraw = {(ep, st): C.null_likelihood(terms[ep], structure=st) for ep in C.EPOCHS for st in ST_MAIN} # (draws, grid) Lbar = {k: v.mean(0) for k, v in Ldraw.items()} Fdraw = {ep: sum(t["q"] * t["c"] * t["o"] for t in terms[ep].values()) for ep in C.EPOCHS} for ep in C.EPOCHS: Lbar[(ep, "origin_all")] = C.null_likelihood(terms[ep], structure="origin_all").mean(0) Lbar[(ep, "model_form")] = np.exp(-Fdraw[ep][:, None] * C.LAM_GRID[None, :]).mean(0) others = [h for h in C.HABITATS if h != "Earth"] part = {(nm, st): C.null_likelihood(terms["today"], only=only, structure=st).mean(0) for st in ST_MAIN for nm, only in (("Earth", ["Earth"]), ("Mars", ["Mars"]), ("others", others))} part[("Earth", "origin_all")] = C.null_likelihood(terms["today"], only=["Earth"], structure="origin_all").mean(0) ceil_ever = C.null_likelihood(terms["today"], ceiling="ever").mean(0) ceil_surv = {st: C.null_likelihood(terms["today"], ceiling="surviving", structure=st).mean(0) for st in ST_MAIN} corner = {} for pt in ("low", "mode", "high"): v1 = C.sample_all(rows, 1, point=pt) for ep in C.EPOCHS: for st in ST_MAIN: corner[(pt, ep, st)] = C.null_likelihood(C.habitat_terms(rows, v1, ep)[0], structure=st)[0] def scaled(ep, st, kc=1.0, ko=1.0): """Mean likelihood with every detection factor multiplied by kc (capped at 1) and every opportunity by ko: a common-mode error in my judgement.""" return C.null_likelihood(terms[ep], structure=st, kc=kc, ko=ko).mean(0) lam, ll = S["lam"], S["log10_lam"].astype(np.float64) fp, frac = meta["sg_fp"], meta["sg_frac"] F_med = float(np.median(Fdraw["today"])) # ------------------------------------------------------------------ scenarios ------------------------------------------------------------------ scen = {"No update (model as published)": ("-", "-", np.ones_like(lam))} for st in ST_MAIN + ["origin_all", "model_form"]: for ep in C.EPOCHS: scen[f"{EP_NAME[ep]}; structure: {ST_NAME[st]}"] = (ep, st, interp_L(Lbar[(ep, st)], ll)) scen[f"Null to date; model's own form with sg_frac fixed at the median exposure F = {F_med:.4f}"] = ("today", "model_form_F_median", np.exp(-F_med * lam)) for st in ST_MAIN: for nm, lab in (("Earth", "Earth's single tree only"), ("Mars", "Mars only"), ("others", "all Solar System bodies, no Earth term")): scen[f"Null to date, {lab}; structure: {ST_NAME[st]}"] = ("today", st, interp_L(part[(nm, st)], ll)) for pt, lab in (("low", "low end"), ("mode", "most likely value"), ("high", "high end")): scen[f"Null to date, every factor at its {lab}; structure: {ST_NAME[st]}"] = ("today", st, interp_L(corner[(pt, "today", st)], ll)) scen[f"Null after the funded missions, every factor at its high end; structure: {ST_NAME[st]}"] = ("funded", st, interp_L(corner[("high", "funded", st)], ll)) for kc, lab in ((3.0, "3 times higher"), (10.0, "10 times higher"), (1 / 3.0, "3 times lower")): scen[f"Null to date, every completeness value {lab}; structure: {ST_NAME[st]}"] = ("today", st, interp_L(scaled("today", st, kc), ll)) for ko, lab in ((3.0, "3 times higher"), (1 / 3.0, "3 times lower")): scen[f"Null to date, every opportunity {lab}; structure: {ST_NAME[st]}"] = ("today", st, interp_L(scaled("today", st, 1.0, ko), ll)) scen[f"Null after the funded missions, every completeness value 3 times higher; structure: {ST_NAME[st]}"] = ("funded", st, interp_L(scaled("funded", st, 3.0), ll)) scen[f"Ceiling: everything that survives or is preserved today is found; structure: {ST_NAME[st]}"] = ("ceiling", st, interp_L(ceil_surv[st], ll)) # hybrids (Earth's factors and the other bodies' factors are separate rows, drawn independently, so the mean likelihoods multiply) scen["Null to date; hybrid: Earth per origin (persistence), other bodies per habitat"] = ("today", "hybrid_1", interp_L(part[("Earth", "origin_persist")] * part[("others", "habitat")], ll)) scen["Null to date; hybrid: Earth per origin (every factor), other bodies per origin (persistence)"] = ("today", "hybrid_2", interp_L(part[("Earth", "origin_all")] * part[("others", "origin_persist")], ll)) scen["Ceiling: any life that ever arose on any body is found (c = 1 everywhere)"] = ("ceiling", "-", interp_L(ceil_ever, ll)) scen["Model's 'thorough search finds none' (sg_frac = 0.3)"] = ("model", "-", np.exp(-frac * lam)) scen["Model's 'independent origin confirmed' (sg_fp = 0.15)"] = ("model", "-", fp + (1 - fp) * -np.expm1(-frac * lam)) out_rows = [] summary = dict(model_md5=meta["model_md5"], n_samples=int(len(lam)), n_per_run=meta["n_per_run"], seeds=meta["seeds"], n_factor_draws=N_DRAW, F_exposure_median_today=F_med, scenarios={}) base_by_school = {sch: stats(S, S["w_" + sch].astype(np.float64)) for sch in SCHOOLS} for name, (ep, st, like) in scen.items(): for sch in SCHOOLS: w0 = S["w_" + sch].astype(np.float64) stt = stats(S, w0 * like) se_a, d_a, se_da = per_seed(S, w0, like, "p_alone") se_f, d_f, se_df = per_seed(S, w0, like, "p_fast") base = base_by_school[sch] out_rows.append(dict(scenario=name, epoch=ep, structure=st, convention=sch, p_fast=stt["p_fast"], p_alone=stt["p_alone"], nobody=stt["nobody"], stuck=stt["stuck"], spread=stt["spread"], d_p_alone=stt["p_alone"] - base["p_alone"], d_p_fast=stt["p_fast"] - base["p_fast"], mc_se_p_alone=se_a, mc_se_d_p_alone=se_da, mc_se_p_fast=se_f, mc_se_d_p_fast=se_df, odds_factor_on_fast=(stt["p_fast"] / (1 - stt["p_fast"])) / (base["p_fast"] / (1 - base["p_fast"])), ess=stt["ess"])) if sch == "site": summary["scenarios"][name] = out_rows[-1] with open(os.path.join(DATA, "update_scenarios.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(out_rows[0])); w.writeheader(); w.writerows(out_rows) w0 = S["w_site"].astype(np.float64) base = base_by_school["site"] # baseline and its Monte Carlo error (seed-to-seed scatter of the absolute value; the shifts are paired differences and far better known) seed_vals = [stats({k: v[S["seed"] == sd] for k, v in S.items()}, w0[S["seed"] == sd]) for sd in np.unique(S["seed"])] summary["baseline"] = dict(p_alone=base["p_alone"], p_fast=base["p_fast"], p_alone_mc_se=float(np.std([v["p_alone"] for v in seed_vals], ddof=1) / np.sqrt(len(seed_vals))), p_fast_mc_se=float(np.std([v["p_fast"] for v in seed_vals], ddof=1) / np.sqrt(len(seed_vals))), p_alone_per_seed=[v["p_alone"] for v in seed_vals]) def shift(like): return stats(S, w0 * like)["p_alone"] - base["p_alone"] # ---- the sg_frac that, put into the model's own likelihood exp(-sg_frac lambda), gives the same shift in P(alone) ---- eff = {} for st in ST_MAIN + ["origin_all"]: for ep in C.EPOCHS: target = shift(interp_L(Lbar[(ep, st)], ll)) x = bisect(lambda lx: target - shift(np.exp(-10 ** lx * lam)), -7.0, 1.0) eff[f"{ep}__{st}"] = dict(d_p_alone=target, equivalent_sg_frac=float(10 ** x), share_of_model_thorough_0p3=float(10 ** x / frac)) summary["equivalent_sg_frac"] = eff # ---- the same update under other prior choices (camps) ---- Cs = dict(np.load(os.path.join(DATA, "model_samples_camps.npz"))) camp_rows = [] for cname in meta["camps"]: cl = Cs[f"{cname}__log10_lam"].astype(np.float64); cw = Cs[f"{cname}__w_site"].astype(np.float64) cal, cf, cst = Cs[f"{cname}__alone"], Cs[f"{cname}__fast"], Cs[f"{cname}__stuck"] lst = [("No update", np.ones_like(cl))] for st in ST_MAIN: lst += [(f"{EP_NAME[ep]}; structure: {ST_NAME[st]}", interp_L(Lbar[(ep, st)], cl)) for ep in C.EPOCHS] lst += [(f"Null to date; structure: {ST_NAME['model_form']}", interp_L(Lbar[("today", "model_form")], cl)), ("Ceiling: any life that ever arose on any body is found", interp_L(ceil_ever, cl)), ("Model's thorough null (sg_frac = 0.3)", np.exp(-frac * 10.0 ** cl))] for label, like in lst: w = cw * like camp_rows.append(dict(camp=cname, scenario=label, p_fast=float(w[cf].sum() / w.sum()), p_alone=float(w[cal].sum() / w.sum()), stuck=float(w[cst].sum() / w.sum()))) with open(os.path.join(DATA, "update_camps.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(camp_rows[0])); w.writeheader(); w.writerows(camp_rows) summary["camps"] = camp_rows # ---- cross-check against the numbers computed inside the model itself (same seeds) ---- chk = {} for key in ("base", "none_thorough", "found"): e = np.array([meta["per_seed"][str(sd)][key]["site"]["evidence"] for sd in meta["seeds"]]) chk[key] = dict(p_alone=float(np.dot([meta["per_seed"][str(sd)][key]["site"]["alone"] for sd in meta["seeds"]], e) / e.sum()), p_fast=float(np.dot([meta["per_seed"][str(sd)][key]["site"]["fast"] for sd in meta["seeds"]], e) / e.sum())) summary["cross_check_model_internal"] = chk # ---- range over INDEPENDENT factor draws only: one posterior per factor draw (binned in lambda, 0.01 dex), per structure ---- idx = np.clip(np.round((ll - LG[0]) / (LG[1] - LG[0])).astype(int), -1, len(LG) - 1) idx[ll < LG[0] - 0.005] = -1 rng_rows = [] for sch in SCHOOLS: ws = S["w_" + sch].astype(np.float64) def binned(mask): b = np.bincount(idx[mask] + 1, weights=ws[mask], minlength=len(LG) + 1) return b[0], b[1:] u_all, W_all = binned(np.ones(len(lam), bool)); u_al, W_al = binned(S["alone"]); u_f, W_f = binned(S["fast"]); u_s, W_s = binned(S["cls"] == 3) for st in ST_MAIN: for ep in C.EPOCHS: den = u_all + Ldraw[(ep, st)] @ W_all pa = (u_al + Ldraw[(ep, st)] @ W_al) / den pf = (u_f + Ldraw[(ep, st)] @ W_f) / den ps = (u_s + Ldraw[(ep, st)] @ W_s) / den for nm, v in (("p_alone", pa), ("p_fast", pf), ("stuck", ps)): q = np.quantile(v, [0.05, 0.25, 0.5, 0.75, 0.95]) rng_rows.append(dict(convention=sch, structure=st, epoch=ep, quantity=nm, p05=q[0], p25=q[1], p50=q[2], p75=q[3], p95=q[4], mean=float(v.mean()), min=float(v.min()), max=float(v.max()))) with open(os.path.join(DATA, "update_range_over_completeness.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(rng_rows[0])); w.writeheader(); w.writerows(rng_rows) summary["range_over_independent_factor_draws_site"] = {f"{r['structure']}__{r['epoch']}__{r['quantity']}": [r["p05"], r["p50"], r["p95"]] for r in rng_rows if r["convention"] == "site"} # ---- the general map: a search that would find a second origin with probability Cs wherever one arose, over exposure F ---- map_rows = [] for F in (0.003, 0.01, 0.03, 0.1, 0.3, 1.0): for Cc in (0.0, 0.03, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 1.0): stt = stats(S, w0 * (1 - Cc * -np.expm1(-F * lam))) map_rows.append(dict(opportunity_F=F, detection_probability_C=Cc, p_fast=stt["p_fast"], p_alone=stt["p_alone"], stuck=stt["stuck"], d_p_alone=stt["p_alone"] - base["p_alone"])) with open(os.path.join(DATA, "ideal_search_map.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(map_rows[0])); w.writeheader(); w.writerows(map_rows) # slope: change in P(alone) per unit of D_fast (the chance the search would have succeeded if life starts easily) D = np.linspace(0, 1, 21) pa = np.array([stats(S, w0 * np.where(S["fast"], 1 - d, 1.0))["p_alone"] for d in D]) summary["p_alone_vs_Dfast"] = dict(D_fast=D.tolist(), p_alone=pa.tolist(), slope_at_zero=float((pa[1] - pa[0]) / (D[1] - D[0]))) # ---- forecasts: chance that the searches turn up a real second origin, under today's posterior; per structure ---- post_today = {st: w0 * interp_L(Lbar[("today", st)], ll) for st in ST_MAIN} fc_rows = [] for st in ST_MAIN: pt = post_today[st] for ep_to, label in (("funded", "funded missions, results by about 2038"), ("unfunded", "funded plus the two unfunded missions, 2050s")): ratio = interp_L(Lbar[(ep_to, st)], ll) / np.maximum(interp_L(Lbar[("today", st)], ll), 1e-300) # P(new searches find nothing | lambda, nothing so far) stt = stats(S, pt * (1 - ratio)) # posterior if a real independent origin is found fc_rows.append(dict(structure=st, searches=label, p_true_detection=float((pt * (1 - ratio)).sum() / pt.sum()), p_true_detection_if_life_starts_easily=float((pt * (1 - ratio))[S["fast"]].sum() / pt[S["fast"]].sum()), p_fast_after_true_detection=stt["p_fast"], p_alone_after_true_detection=stt["p_alone"], stuck_after_true_detection=stt["stuck"])) # what the evidence to date would have shown, a priori (before any search) like0 = interp_L(Lbar[("today", st)], ll) fc_rows.append(dict(structure=st, searches="everything done to date (retrospective: chance it would have found one)", p_true_detection=float((w0 * (1 - like0)).sum() / w0.sum()), p_true_detection_if_life_starts_easily=float((w0 * (1 - like0))[S["fast"]].sum() / w0[S["fast"]].sum()), p_fast_after_true_detection=float("nan"), p_alone_after_true_detection=float("nan"), stuck_after_true_detection=float("nan"))) with open(os.path.join(DATA, "detection_forecast.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(fc_rows[0])); w.writeheader(); w.writerows(fc_rows) summary["detection_forecast"] = fc_rows # ---- how much an accepted claim is worth, as a function of the false-alarm rate (funded missions); the base-rate argument of Foote et al. 2023 ---- cl_rows = [] for fpr in (0.0, 0.001, 0.003, 0.01, 0.03, 0.05, 0.15, 0.30): for st in ST_MAIN: pt = post_today[st] ratio = interp_L(Lbar[("funded", st)], ll) / np.maximum(interp_L(Lbar[("today", st)], ll), 1e-300) p_true = float((pt * (1 - ratio)).sum() / pt.sum()) stt = stats(S, pt * (fpr + (1 - fpr) * (1 - ratio))) p_claim = float((pt * (fpr + (1 - fpr) * (1 - ratio))).sum() / pt.sum()) cl_rows.append(dict(search=f"funded missions, itemised completeness, structure: {ST_NAME[st]}", structure=st, false_alarm_rate=fpr, p_claim=p_claim, share_of_claims_real=float((1 - fpr) * p_true / max(p_claim, 1e-300)), p_fast_after_claim=stt["p_fast"], p_alone_after_claim=stt["p_alone"], stuck_after_claim=stt["stuck"])) stt = stats(S, w0 * (fpr + (1 - fpr) * -np.expm1(-frac * lam))) p_claim = float((w0 * (fpr + (1 - fpr) * -np.expm1(-frac * lam))).sum() / w0.sum()) p_true = float((w0 * -np.expm1(-frac * lam)).sum() / w0.sum()) cl_rows.append(dict(search="model's thorough search (sg_frac = 0.3)", structure="model", false_alarm_rate=fpr, p_claim=p_claim, share_of_claims_real=float((1 - fpr) * p_true / max(p_claim, 1e-300)), p_fast_after_claim=stt["p_fast"], p_alone_after_claim=stt["p_alone"], stuck_after_claim=stt["stuck"])) with open(os.path.join(DATA, "claim_value.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(cl_rows[0])); w.writeheader(); w.writerows(cl_rows) summary["claim_value"] = cl_rows # ---- Part 3: is Martian life, if found, an independent origin? ---- transfer = [] for route, rocks, unc in TRANSFER: t = dict(route=route, rocks=rocks, unc=unc) for ps in (1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-3): t[f"T_if_seed_prob_{ps:g}"] = float(-np.expm1(-rocks * ps)) t["seed_prob_for_T_half"] = float(np.log(2) / rocks) t["source"] = "Worth, Sigurdsson & House 2013, Astrobiology 13, 1155, published Tables 4 and 5 (PMC3870607)" transfer.append(t) with open(os.path.join(DATA, "transfer_literature.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(transfer[0])); w.writeheader(); w.writerows(transfer) summary["transfer_ratio_mars_to_earth_over_earth_to_mars"] = TRANSFER[1][1] / TRANSFER[0][1] # Shared ancestry can arise by EITHER route: # T_em probability that Earth life seeded Mars and took hold there (whether or not Mars had a native origin) # T_me probability that Earth's own tree descends from Martian life, given that Mars had a native origin # (Mars-to-Earth transfer is 58 times more efficient in rocks delivered than Earth-to-Mars) # Native Martian origin: pM(lambda) = q (1 - exp(-o lambda)). Life is found on Mars if pM or seeding: pM + (1 - pM) T_em. # It is an independent origin if Mars had a native origin, Earth's tree is not descended from it (1 - T_me), and the find is the # native lineage rather than an Earth-seeded one (probability one half if both are present). # P(independent | found) = / # The result depends strongly on the prior for a native origin: it is given under the model's base posterior and within the fast camp only. pt = post_today["habitat"] qM, oM = terms["today"]["Mars"]["q"], terms["today"]["Mars"]["o"] pM_draw = qM[:, None] * -np.expm1(-oM[:, None] * C.LAM_GRID[None, :]) pM = np.interp(ll, LG, pM_draw.mean(0), left=0.0) ind_rows = [] T_grid = (0.0, 0.01, 0.03, 0.1, 0.3, 0.5, 0.9) for prior, m in (("model's base posterior (after the null to date)", np.ones(len(lam), bool)), ("fast camp only (life starts easily)", S["fast"])): p_nat = float((pt * pM)[m].sum() / pt[m].sum()) for T_em in T_grid: for T_me_lab, T_me in (("0", 0.0), ("same as T_em", T_em), ("0.5", 0.5), ("0.9", 0.9)): num = (pt * pM * (1 - T_em / 2) * (1 - T_me))[m].sum() den = (pt * (pM + (1 - pM) * T_em))[m].sum() ind_rows.append(dict(body="Mars", prior=prior, p_native_origin=p_nat, T_earth_to_target=T_em, T_target_to_earth_label=T_me_lab, T_target_to_earth=T_me, p_inhabited_ever=float(den / pt[m].sum()), p_independent_given_life_found=float(num / max(den, 1e-300)), note="")) # per-rock parametrisation: the same per-rock probability of taking hold in both directions (rock counts from Worth et al. 2013) for ps in (1e-9, 1e-8, 1e-7, 1e-6, 1e-5): T_em = float(-np.expm1(-TRANSFER[0][1] * ps)); T_me = float(-np.expm1(-TRANSFER[1][1] * ps)) for prior, m in (("model's base posterior (after the null to date)", np.ones(len(lam), bool)), ("fast camp only (life starts easily)", S["fast"])): num = (pt * pM * (1 - T_em / 2) * (1 - T_me))[m].sum(); den = (pt * (pM + (1 - pM) * T_em))[m].sum() ind_rows.append(dict(body="Mars", prior=prior, p_native_origin=float((pt * pM)[m].sum() / pt[m].sum()), T_earth_to_target=T_em, T_target_to_earth_label=f"per-rock probability {ps:g} both ways", T_target_to_earth=T_me, p_inhabited_ever=float(den / pt[m].sum()), p_independent_given_life_found=float(num / max(den, 1e-300)), note="T_me here is the chance Martian life took hold on Earth: an upper bound on the chance that Earth's tree descends from it")) # ocean worlds: no published estimate of the establishment probability exists; the T values are ASSUMPTIONS spanning four decades for body in ("Europa", "Enceladus", "Titan"): qh, oh = terms["today"][body]["q"], terms["today"][body]["o"] ph = np.interp(ll, LG, (qh[:, None] * -np.expm1(-oh[:, None] * C.LAM_GRID[None, :])).mean(0), left=0.0) for prior, m in (("model's base posterior (after the null to date)", np.ones(len(lam), bool)), ("fast camp only (life starts easily)", S["fast"])): for T in (1e-5, 1e-4, 1e-3, 1e-2): num = (pt * ph * (1 - T / 2))[m].sum(); den = (pt * (ph + (1 - ph) * T))[m].sum() ind_rows.append(dict(body=body, prior=prior, p_native_origin=float((pt * ph)[m].sum() / pt[m].sum()), T_earth_to_target=T, T_target_to_earth_label="0", T_target_to_earth=0.0, p_inhabited_ever=float(den / pt[m].sum()), p_independent_given_life_found=float(num / max(den, 1e-300)), note="T is an assumption, not a published estimate (Worth et al. 2013: about 6 to 26 rocks to Europa, about 1 each to Titan and Enceladus in 3.5 Gyr)")) with open(os.path.join(DATA, "independence.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(ind_rows[0])); w.writeheader(); w.writerows(ind_rows) summary["independence"] = ind_rows # ---- Is the Mars record a null at all? Cheyava Falls scored as positive evidence ---- # J = 'the Bright Angel mudstone hosted microbial metabolism'. P(J | lambda) = sB * [1 - (1 - pM)(1 - T)], # sB = P(a Martian biosphere, if there was one, was living in that mud), from the factor table (row Mars / sB); T = Earth-to-Mars seeding. # The observation multiplies the likelihood by 1 + (LR - 1) * P(J | lambda), LR = P(features | J) / P(features | not J). LR = 1 is the pure null. # The factor draws are kept paired: is averaged draw by draw, not as a product of means. sB_draw = vals[("today", "Mars", "-", "-", "sB")] summary["cheyava_sB"] = dict(q=C.qs(sB_draw), mean=float(sB_draw.mean())) LP = {st: (Ldraw[("today", st)] * pM_draw).mean(0) for st in ST_MAIN} # LPs = {st: (Ldraw[("today", st)] * pM_draw * sB_draw[:, None]).mean(0) for st in ST_MAIN} # Ls_ = {st: (Ldraw[("today", st)] * sB_draw[:, None]).mean(0) for st in ST_MAIN} # def cheyava_weight(st, LR, T=0.0, sB=None): """Per-sample likelihood of (null elsewhere + Cheyava Falls features) under structure st.""" L0 = interp_L(Lbar[("today", st)], ll) if sB is None: LJ = interp_L(LPs[st], ll, left=0.0) + T * (interp_L(Ls_[st], ll, left=float(sB_draw.mean())) - interp_L(LPs[st], ll, left=0.0)) else: LJ = sB * (interp_L(LP[st], ll, left=0.0) + T * (L0 - interp_L(LP[st], ll, left=0.0))) return L0 + (LR - 1) * LJ, LJ cf_rows, be_rows = [], [] for st in ST_MAIN: for sB_lab, sB in [("factor table (0.1 / 0.4 / 0.8)", None)] + [(f"{v:g}", v) for v in CHEYAVA_SB_FIXED]: for T in (0.0, 0.1, 0.3): for LR in CHEYAVA_LR: if (sB is not None or T > 0) and LR not in (1.0, 1.5, 2.0, 3.0, 10.0): continue wgt, LJ = cheyava_weight(st, LR, T, sB) stt = stats(S, w0 * wgt) cf_rows.append(dict(structure=st, sB=sB_lab, transfer_probability_T=T, likelihood_ratio_for_past_life_at_cheyava=LR, p_fast=stt["p_fast"], p_alone=stt["p_alone"], stuck=stt["stuck"], d_p_alone_vs_no_update=stt["p_alone"] - base["p_alone"], p_past_life_at_site=float((w0 * LR * LJ).sum() / (w0 * wgt).sum()))) be = bisect(lambda LR: stats(S, w0 * cheyava_weight(st, LR, T, sB)[0])["p_alone"] - base["p_alone"], 1.0, 50.0) be_rows.append(dict(structure=st, sB=sB_lab, transfer_probability_T=T, breakeven_likelihood_ratio=be, d_p_alone_pure_null=shift(interp_L(Lbar[("today", st)], ll)))) with open(os.path.join(DATA, "cheyava_sensitivity.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(cf_rows[0])); w.writeheader(); w.writerows(cf_rows) with open(os.path.join(DATA, "cheyava_breakeven.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(be_rows[0])); w.writeheader(); w.writerows(be_rows) summary["cheyava_sensitivity"] = cf_rows summary["cheyava_breakeven"] = be_rows # ---- the envelope: net change in P(alone) from the evidence in hand, over structure and over the scoring of Cheyava Falls ---- main_rows = [r for r in cf_rows if r["sB"].startswith("factor") and r["transfer_probability_T"] == 0.0] env = {} for lr_max in (1.0, 2.0, 3.0): v = [r["d_p_alone_vs_no_update"] for r in main_rows if r["likelihood_ratio_for_past_life_at_cheyava"] <= lr_max] env[f"two_structures_LR_1_to_{lr_max:g}"] = [min(v), max(v)] mf = summary["scenarios"][f"{EP_NAME['today']}; structure: {ST_NAME['model_form']}"]["d_p_alone"] env["upper_bound_models_own_form_pure_null"] = mf env["overall_LR_1_to_2_including_models_own_form"] = [env["two_structures_LR_1_to_2"][0], max(env["two_structures_LR_1_to_2"][1], mf)] summary["net_update_envelope_d_p_alone"] = env json.dump(summary, open(os.path.join(DATA, "update_summary.json"), "w"), indent=1) # ------------------------------------------------------------------ printout ------------------------------------------------------------------ b = summary["baseline"] print(f"baseline: P(alone) {b['p_alone']:.4f} +/- {b['p_alone_mc_se']:.4f} (Monte Carlo, {len(seed_vals)} seeds; per seed {np.round(b['p_alone_per_seed'], 4).tolist()}); P(fast) {b['p_fast']:.4f} +/- {b['p_fast_mc_se']:.4f}") print(f"\n{'scenario':112s} P(fast) P(alone) d P(alone) (MC se) stuck") for r in out_rows: if r["convention"] == "site": print(f"{r['scenario'][:112]:112s} {r['p_fast']:.4f} {r['p_alone']:.4f} {r['d_p_alone']:+.4f} ({r['mc_se_d_p_alone']:.4f}) {r['stuck']:.4f}") print("\ncross-check, numbers computed inside the model (same seeds, evidence-weighted):", {k: {a: round(v_, 4) for a, v_ in v.items()} for k, v in chk.items()}) print("\nequivalent sg_frac (the value that gives the same shift in the model's own form):") for k, v in eff.items(): print(f" {k:28s} d {v['d_p_alone']:+.4f} sg_frac {v['equivalent_sg_frac']:.5f} = {100 * v['share_of_model_thorough_0p3']:.2f}% of 0.3") print("\nrange over independent factor draws only (Earth-as-given):") for r in rng_rows: if r["convention"] == "site" and r["quantity"] != "stuck": print(f" {r['structure']:15s} {r['epoch']:9s} {r['quantity']:8s} median {r['p50']:.4f} 5-95% [{r['p05']:.4f}, {r['p95']:.4f}] min {r['min']:.4f} max {r['max']:.4f}") print("\nother conventions, null to date:") for r in out_rows: if r["scenario"].startswith("Null to date; structure: per") or r["scenario"].startswith("No update"): print(f" {r['convention']:8s} {r['scenario'][:60]:60s} P(fast) {r['p_fast']:.4f} P(alone) {r['p_alone']:.4f} d {r['d_p_alone']:+.4f}") print("\ncamps:"); [print(f" {r['camp'][:52]:52s} {r['scenario'][:70]:70s} P(fast) {r['p_fast']:.4f} P(alone) {r['p_alone']:.4f}") for r in camp_rows] print("\nslope of P(alone) per unit D_fast at zero:", round(summary["p_alone_vs_Dfast"]["slope_at_zero"], 4), " P(alone) at D_fast = 1:", round(pa[-1], 4)) rnd = lambda r: {k: (round(v, 4) if isinstance(v, float) else v) for k, v in r.items()} print("\nforecast:"); [print(" ", rnd(r)) for r in fc_rows] print("\nvalue of an accepted claim:"); [print(" ", rnd(r)) for r in cl_rows] print("\nindependence (Mars, T grid):") for r in ind_rows: if r["body"] == "Mars" and r["T_earth_to_target"] in (0.01, 0.1, 0.3, 0.9) or r["body"] != "Mars" or "per-rock" in r["T_target_to_earth_label"]: print(f" {r['body']:9s} {r['prior'][:22]:22s} P(native) {r['p_native_origin']:.3f} T_em {r['T_earth_to_target']:.3g} T_me {r['T_target_to_earth_label'][:38]:38s} ({r['T_target_to_earth']:.3g}) P(indep | found) {r['p_independent_given_life_found']:.4f}") print("\nCheyava Falls (sB from the factor table:", summary["cheyava_sB"], "):") for r in cf_rows: if r["sB"].startswith("factor"): print(f" {r['structure']:15s} T {r['transfer_probability_T']:.1f} LR {r['likelihood_ratio_for_past_life_at_cheyava']:5.2f} P(alone) {r['p_alone']:.4f} net d {r['d_p_alone_vs_no_update']:+.4f} P(fast) {r['p_fast']:.4f} P(past life at the site) {r['p_past_life_at_site']:.3f}") print("break-even likelihood ratios:") for r in be_rows: if r["transfer_probability_T"] == 0.0: print(f" {r['structure']:15s} sB {r['sB']:32s} break-even LR {r['breakeven_likelihood_ratio']:.3f} (pure null d {r['d_p_alone_pure_null']:+.4f})") print("\nnet update envelope:", env) print("\ntransfer:"); [print(" ", t["route"], t["rocks"], "seed prob for T=0.5:", f"{t['seed_prob_for_T_half']:.2g}") for t in transfer] if __name__ == "__main__": main()