"""Driver for model4 (version 4.1): columns, decomposition (all orders, pooled seeds), camps, one-dimensional levers, future observations (pooled seeds). -> integrate4_results.json. Usage: python3 run4.py [n_per_run]""" import os os.environ.setdefault("OMP_NUM_THREADS", "1"); os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import itertools import json import math import sys import time from multiprocessing import Pool import numpy as np import model4 as M HERE = os.path.dirname(os.path.abspath(__file__)) B = M.BASE S = M.SCHOOLS MILLS = dict(w_fast=0.6, m_weights={0: 0.6, 1: 0.25, 2: 0.15}, step_decades=(2.0, 6.0), f_site=(0.03, 1.0), f_persist=(0.3, 1.0), f_gate=(0.1, 1.0), f_c=(0.03, 1.0)) CAMPS = { "Base (modeller's judgement)": {}, "Life starts easily (chemistry optimist)": dict(w_fast=0.6), "Life starts easily, 90% sure": dict(w_fast=0.9), "Life starts easily, as in model 3 (post-data value)": dict(w_fast=0.4), "Biology optimist": dict(w_fast=0.9, m_weights={0: 0.7, 1: 0.2, 2: 0.1}, f_site=(0.1, 1.0), f_persist=(0.3, 1.0), f_c=(0.1, 1.0)), # round 3, opposing camp: the published alternative to hard steps (Mills et al. 2025), and the classic SETI position "Gradual evolution, no hard steps (Mills et al. 2025)": dict(**MILLS), "Classic SETI: life common, expansion doubtful": dict(q_barrier=0.40, W_exp=(2.0, 20.0), **MILLS), "Origin of life no rarer than once per observable universe": dict(slow_decades=(3.0, 22.0), w_fast=0.28), "Early life at the timing study's strength": dict(fast_lo=0.75, tau_over_T=(0.018, 0.125)), "Early life credited at 6 to 1": dict(early_bf=6.0), "Early life credited at 9 to 1": dict(early_bf=9.0), "Early life credited at 13 to 1": dict(early_bf=13.0), "One or two shallow steps": dict(m_weights={0: 0.15, 1: 0.35, 2: 0.30, 3: 0.13, 4: 0.07}), "Flat over 0 to 4 steps": dict(m_weights={0: 0.2, 1: 0.2, 2: 0.2, 3: 0.2, 4: 0.2}), "Carter / Snyder-Beattie hard steps": dict(w_fast=0.1, m_weights={3: 0.3, 4: 0.3, 5: 0.25, 6: 0.15}, step_decades=(2.0, 12.0)), "Strict ignorance about rates": dict(w_fast=0.1, m_weights={1: 0.2, 2: 0.3, 3: 0.3, 4: 0.2}, step_decades=(6.0, 12.0)), "Rare Earth": dict(f_site=(1e-5, 1.0), f_persist=(1e-3, 1.0)), "Life dated from LUCA (very early)": dict(tau_over_T=(0.03, 0.06)), "Life dated from accepted fossils only": dict(tau_over_T=(0.16, 0.20)), "Expansion is easy (it only takes one)": dict(q_barrier=0.02, W_exp=(1.0, 4.0)), "Expansion sceptic": dict(q_barrier=0.40, W_exp=(2.0, 20.0)), "Settlement would be obvious (Hart)": dict(q_cc=0.01, p_stealth=0.003, c_claim=("uniform", 0.5, 1.0)), "Hart, crossed: easy and obvious": dict(q_barrier=0.02, W_exp=(1.0, 4.0), q_cc=0.01, p_stealth=0.003, c_claim=("uniform", 0.5, 1.0)), "Expansionist (fast, intergalactic)": dict(q_barrier=0.02, W_exp=(1.0, 4.0), v_ig=("uniform", 0.1, 0.99), q_ig=0.9), "Expansionist with Carter steps": dict(q_barrier=0.02, W_exp=(1.0, 4.0), v_ig=("uniform", 0.1, 0.99), q_ig=0.9, w_fast=0.1, m_weights={3: 0.3, 4: 0.3, 5: 0.25, 6: 0.15}, step_decades=(2.0, 12.0)), "Expansionist, Carter steps, long-lived planets": dict(q_barrier=0.02, W_exp=(1.0, 4.0), v_ig=("uniform", 0.1, 0.99), q_ig=0.9, w_fast=0.1, m_weights={3: 0.3, 4: 0.3, 5: 0.25, 6: 0.15}, step_decades=(2.0, 12.0), p_cap_short=0.0, L_cap_long=(15.0, 300.0), K_uv=(1.0, 1.0), L_K_hz=(1.0, 1.0), p_persist_hazard=0.0), "Low-footprint settlement likely": dict(q_cc=0.30), "Most lineages are discreet": dict(c_claim=("loguniform", 1e-6, 1.0)), "Zoo sympathiser": dict(p_stealth=0.20), "Planets stay habitable for tens of Gyr": dict(p_cap_short=0.3, L_cap_long=(15.0, 300.0)), "Long-lived hosts count fully": dict(p_cap_short=0.0, L_cap_long=(15.0, 300.0), K_uv=(1.0, 1.0), M_uv=(1.0, 1.0), L_K_hz=(1.0, 1.0), p_persist_hazard=0.0), "Stable planets are a fixed fraction": dict(p_persist_hazard=0.0), "Persistence is a constant hazard": dict(p_persist_hazard=1.0), "No ultraviolet penalty on K and M dwarfs": dict(K_uv=(1.0, 1.0), M_uv=(1.0, 1.0)), "Loud beacons impossible (radio searches see nothing)": dict(f_beacon=(1e-12, 1e-11)), } LEVERS = { "w_fast": [0.05, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 0.9, 0.95], "q_barrier": [0.0, 0.05, 0.15, 0.30, 0.50], "q_cc": [0.0, 0.02, 0.05, 0.10, 0.25, 0.50], "p_stealth": [0.003, 0.03, 0.10, 0.20], "obs_R": [1e3, 1e6, 1e10], "q_ig": [0.2, 0.6, 0.9], } FUTURE = { "Independent origin of life confirmed in the Solar System": dict(extra_obs="second_genesis"), "Thorough search finds no second origin in the Solar System": dict(extra_obs="second_genesis_none"), "HWO: 25 habitable-zone planets, at least one biosignature": dict(extra_obs="hwo_some"), "HWO: 25 habitable-zone planets, no biosignature": dict(extra_obs="hwo_none"), "Earth-level technosignatures excluded for the nearest million stars": dict(extra_obs="techno_null_1e6"), "Earth-level technosignatures excluded for the whole Galaxy": dict(extra_obs="techno_null_galaxy"), "Thorough Solar System artefact survey finds nothing": dict(extra_obs="artefact_null"), "Sustained interstellar expansion shown feasible for us": dict(extra_obs="exp_feasible"), "Sustained interstellar expansion shown to be common (over 1 civilisation in 1000)": dict(extra_obs="exp_common"), "Sustained interstellar expansion shown infeasible": dict(extra_obs="exp_infeasible"), "M-dwarf planets shown habitable": dict(extra_obs="mdwarf_alive"), "M-dwarf planets shown dead": dict(extra_obs="mdwarf_dead"), "Rocky planets shown to stay habitable beyond 15 Gyr": dict(extra_obs="planets_long_lived"), "Rocky planets shown to lose habitability within 15 Gyr": dict(extra_obs="planets_short_lived"), } FUTURE_RANGE = { # the same observation under the other end of its error-rate range "second_genesis": [dict(sg_fp=0.05), dict(sg_fp=0.30)], "hwo_some": [dict(hwo_fp=0.03), dict(hwo_fp=0.20)], "artefact_null": [dict(art_d=0.1), dict(art_d=0.3), dict(art_d=0.8), dict(art_d=0.95)], } TERMS = M.TERMS SEEDS = list(range(101, 109)) FUT_LIST = [] # (key, extra_obs, overrides) evaluated inside every base run, so forecasts pool the same seeds for _name, _kw in FUTURE.items(): FUT_LIST.append((_name, _kw["extra_obs"], {})) for _alt in FUTURE_RANGE.get(_kw["extra_obs"], []): FUT_LIST.append((_name + " || " + json.dumps(_alt), _kw["extra_obs"], _alt)) CAMP_SEEDS = (7, 8) def slim(r): o = dict(prior=r["prior"], schools=r["schools"], prior_stats=r["prior_stats"]) for k in ("futures", "ablation"): if k in r: o[k] = r[k] return o def job(args): name, kw = args return name, slim(M.run(M.vary(B, **kw))) def line(name, r, s="site"): q = r["schools"][s] print(f"{name[:58]:58s} " + " ".join(f"{q['posterior'][c]:5.3f}" for c in M.CLASSES) + f" | alone {q['p_no_other_ever_in_galaxy']:.3f} now {q['p_anyone_now_in_galaxy']:.2f} ESS {q['ess']:.0f}", flush=True) def odds(p): return p / max(1 - p, 1e-12) if __name__ == "__main__": n = int(sys.argv[1]) if len(sys.argv) > 1 else 250_000 t0 = time.time() jobs = [] for sd in SEEDS: jobs.append((("base", sd), dict(n=n, seed=sd, ablate=True, futures=FUT_LIST))) jobs.append((("sia", sd), dict(n=n, seed=sd, tilt=True, futures=FUT_LIST))) for name, kw in CAMPS.items(): for sd in CAMP_SEEDS: jobs.append((("camp", name, sd), dict(n=n, seed=sd, **kw))) for key, vals in LEVERS.items(): for v in vals: jobs.append((("lever", key, v), dict(n=n, **{key: v}))) for qb in LEVERS["q_barrier"]: for qc in LEVERS["q_cc"]: jobs.append((("grid", qb, qc), dict(n=n, q_barrier=qb, q_cc=qc))) print("jobs", len(jobs), flush=True) with Pool(int(os.environ.get("FERMI_PROCS", "4"))) as pool: res = dict(pool.imap_unordered(job, jobs, chunksize=1)) print("runs done in", round(time.time() - t0), "s", flush=True) out = dict(classes=M.CLASSES, class_labels=M.CLASS_LABELS, schools=S, school_labels=M.SCHOOL_LABELS, model_version="4.1", base_params={k: v for k, v in B.items()}, n_per_run=n, n_seeds=len(SEEDS)) reps = [res[("base", sd)] for sd in SEEDS] sreps = [res[("sia", sd)] for sd in SEEDS] base = M.pool(reps) spool = M.pool(sreps) base["schools"]["sia"] = spool["schools"]["sia"] # the SIA column comes from the tilted sampler futs = base.pop("futures"); abl = base.pop("ablation") out["base"] = base out["mc_sd"] = {s: {c: float(np.std([r["schools"][s]["posterior"][c] for r in (sreps if s == "sia" else reps)], ddof=1) / math.sqrt(len(reps))) for c in M.CLASSES} for s in S} out["mc_sd_alone"] = {s: float(np.std([r["schools"][s]["p_no_other_ever_in_galaxy"] for r in reps], ddof=1) / math.sqrt(len(reps))) for s in S if s != "sia"} print(" " * 59 + " ".join(f"{c[:5]:5s}" for c in M.CLASSES)) print("PRIOR".ljust(59) + " ".join(f"{base['prior'][c]:5.3f}" for c in M.CLASSES) + f" | alone {base['prior_stats']['p_no_other_ever_in_galaxy']:.3f}") for s in S: line("BASE " + s, base, s) # ---- decomposition: every subset of the four observation terms on the pooled base seeds (Earth-as-given column) ---- out["ablation"] = abl nobody = lambda q: q["posterior"]["empty_universe"] + q["posterior"]["rare"] measures = {"alone": lambda q: q["p_no_other_ever_in_galaxy"], "nobody": nobody, "logodds_alone": lambda q: math.log(odds(q["p_no_other_ever_in_galaxy"])), "logodds_nobody": lambda q: math.log(odds(nobody(q)))} shap = {} for mname, f in measures.items(): val = lambda on: f(abl[M.subset_key(on)]["schools"]["site"]) sh = {t: 0.0 for t in TERMS}; lo = {t: math.inf for t in TERMS}; hi = {t: -math.inf for t in TERMS} for perm in itertools.permutations(TERMS): for i, t in enumerate(perm): dlt = val(set(perm[:i + 1])) - val(set(perm[:i])) sh[t] += dlt / math.factorial(len(TERMS)); lo[t] = min(lo[t], dlt); hi[t] = max(hi[t], dlt) shap[mname] = dict(start=val(set()), end=val(set(TERMS)), contributions=sh, lowest=lo, highest=hi) out["shapley_alone"] = dict(site=shap["alone"]) # same layout as before: Shapley values of P(no other civilisation ever in the Galaxy) out["shapley"] = shap # also for the class partition, and on the log-odds scale (exp = odds multiplier) pf = lambda k: abl[k]["schools"]["site"]["p_fast"] out["early_life_factor"] = dict(alone=odds(pf("early")) / odds(pf("prior")), after_late=odds(pf("early+late")) / odds(pf("late")), after_all=odds(pf("early+late+local+surveys")) / odds(pf("late+local+surveys"))) print("\nShapley contributions (Earth as given):") for mname, d in shap.items(): print(f" {mname:15s} start {d['start']:.3f} " + " ".join(f"{t} {v:+.3f} [{d['lowest'][t]:+.3f},{d['highest'][t]:+.3f}]" for t, v in d["contributions"].items()) + f" end {d['end']:.3f}") print(" odds multipliers on P(alone):", {t: round(math.exp(v), 2) for t, v in shap["logodds_alone"]["contributions"].items()}, "total", round(math.exp(shap["logodds_alone"]["end"] - shap["logodds_alone"]["start"]), 2)) print(" early-life factor on the odds of a fast origin:", {k: round(v, 2) for k, v in out["early_life_factor"].items()}) print() for nm in ("prior", "early", "late", "early+late", "local", "early+late+local", "early+late+local+surveys"): line("terms: " + nm, abl[nm]) out["camps"] = {} print() for name, kw in CAMPS.items(): r = M.pool([res[("camp", name, sd)] for sd in CAMP_SEEDS]) out["camps"][name] = dict(overrides={k: (list(v) if isinstance(v, tuple) else v) for k, v in kw.items()}, **r) line(name, r); line(" ssa_all", r, "ssa_all") out["levers"] = {key: [dict(value=v, **res[("lever", key, v)]) for v in vals] for key, vals in LEVERS.items()} print() for key, vals in LEVERS.items(): for v in vals: line(f"{key} = {v}", res[("lever", key, v)], "ssa_obs" if key == "obs_R" else "site") out["grid_barrier_cc"] = [dict(q_barrier=qb, q_cc=qc, **{s: res[("grid", qb, qc)]["schools"][s]["p_no_other_ever_in_galaxy"] for s in ("site", "ssa_all")}, stuck=res[("grid", qb, qc)]["schools"]["site"]["posterior"]["stuck"], unnoticed=res[("grid", qb, qc)]["schools"]["site"]["posterior"]["unnoticed"]) for qb in LEVERS["q_barrier"] for qc in LEVERS["q_cc"]] out["future"] = {} print() def fut(key): r = futs[key] r["schools"]["sia"] = spool["futures"][key]["schools"]["sia"] return r for name, kw in FUTURE.items(): r = fut(name) alts = [dict(overrides=a, **fut(name + " || " + json.dumps(a))) for a in FUTURE_RANGE.get(kw["extra_obs"], [])] out["future"][name] = dict(extra_obs=kw["extra_obs"], alternatives=alts, **r) line(name, r); print(" " * 10 + "forecast probability: " + " ".join(f"{s} {r['schools'][s]['predictive']:.2f}" for s in S) + (" without false alarms: " + " ".join(f"{s} {r['schools'][s]['predictive_true']:.2f}" for s in S) if "predictive_true" in r["schools"]["site"] else "")) json.dump(out, open(os.path.join(HERE, "integrate4_results.json"), "w"), indent=1) print("done in", round(time.time() - t0), "s")