""" Step 3. Draw posterior samples from the project's model (version 4.1) so that null results of any completeness can be evaluated by re-weighting, without touching the model. The model file is COPIED from models/ to the scratch directory and imported from there; it is never modified (its md5 is recorded). Each run uses at most 60,000 samples (the limit set for this study); eight seeds are pooled, the same seeds as the model's own driver (run4.py, seeds 101-108). Output: ../data/model_samples.npz with, per sample, log10 lambda (lambda = expected independent origins of life per 5.6-Gyr window of an Earth-like planet), the fast/slow flag, the class, N_ever, and today's posterior weight under four reference-class conventions (site, ssa_all, ssa_now, ssa_obs). The SIA column needs the tilted sampler and is not reproduced here. Usage: python3 model_samples.py [n_per_run=60000] (run under nice -n 10; about 10 s per seed) """ import hashlib import json import os import shutil import sys os.environ.setdefault("OMP_NUM_THREADS", "1") os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "..", "data") RAW = os.environ.get("SOS_RAW", "fermi_raw/second_origin_search") SRC = "models" COPY = os.path.join(RAW, "model_copy") SEEDS = list(range(101, 109)) SCHOOLS = ["site", "ssa_all", "ssa_now", "ssa_obs"] # prior choices taken from the model's own list of camps (models/run4.py), plus one that narrows the fast camp CAMPS = { "Life starts easily, 60% prior (chemistry optimist)": dict(w_fast=0.6), "Early life credited at 13 to 1 (Kipping 2025)": dict(early_bf=13.0), "Gradual evolution, no hard steps (Mills et al. 2025)": 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)), "Fast camp narrowed to 1-10 origins per window": dict(fast_decades=1.0), } def main(): n = int(sys.argv[1]) if len(sys.argv) > 1 else 60_000 assert n <= 60_000, "this study is limited to 60,000 samples per run" os.makedirs(COPY, exist_ok=True) for fn in ("model4.py", "chain4.py", "cosmo_tables.py", "cosmo_tables.npz"): shutil.copyfile(os.path.join(SRC, fn), os.path.join(COPY, fn)) md5 = hashlib.md5(open(os.path.join(COPY, "model4.py"), "rb").read()).hexdigest() sys.path.insert(0, COPY) import model4 as M out = {k: [] for k in ["log10_lam", "fast", "cls", "alone", "N_ever", "seed", "barrier"] + ["w_" + s for s in SCHOOLS]} check = {} # the model's own two hypothetical observations, evaluated inside the model, as a cross-check of the re-weighting futs = [("none_thorough", "second_genesis_none", {}), ("found", "second_genesis", {})] for sd in SEEDS: r = M.run(M.vary(M.BASE, n=n, seed=sd, futures=futs), keep=True) S = r["_s"] out["log10_lam"].append(np.log10(np.maximum(S["lam"], 1e-300)).astype(np.float32)) out["fast"].append(S["fast"].astype(np.bool_)) out["cls"].append(S["cls"].astype(np.int8)) out["alone"].append((S["N_ever"] < 1).astype(np.bool_)) out["N_ever"].append(S["N_ever"].astype(np.float32)) out["barrier"].append(S["barrier"].astype(np.bool_)) out["seed"].append(np.full(n, sd, np.int16)) for s in SCHOOLS: out["w_" + s].append(S["weights"][s].astype(np.float32)) check[sd] = dict( base={s: dict(alone=r["schools"][s]["p_no_other_ever_in_galaxy"], fast=r["schools"][s]["p_fast"], evidence=r["schools"][s]["evidence"], ess=r["schools"][s]["ess"]) for s in SCHOOLS}, none_thorough={s: dict(alone=r["futures"]["none_thorough"]["schools"][s]["p_no_other_ever_in_galaxy"], fast=r["futures"]["none_thorough"]["schools"][s]["p_fast"], evidence=r["futures"]["none_thorough"]["schools"][s]["evidence"]) for s in SCHOOLS}, found={s: dict(alone=r["futures"]["found"]["schools"][s]["p_no_other_ever_in_galaxy"], fast=r["futures"]["found"]["schools"][s]["p_fast"], evidence=r["futures"]["found"]["schools"][s]["evidence"]) for s in SCHOOLS}, prior=dict(alone=r["prior_stats"]["p_no_other_ever_in_galaxy"], fast=r["prior_stats"]["p_fast"])) print(f"seed {sd}: P(alone) {check[sd]['base']['site']['alone']:.4f} P(fast) {check[sd]['base']['site']['fast']:.4f} ESS {check[sd]['base']['site']['ess']:.0f}", flush=True) arr = {k: np.concatenate(v) for k, v in out.items()} np.savez_compressed(os.path.join(DATA, "model_samples.npz"), **arr) # ---- other camps (prior choices the model's own driver lists), Earth-as-given weights only, four seeds each ---- camp_meta = {} carr = {} for cname, kw in CAMPS.items(): parts = {k: [] for k in ("log10_lam", "fast", "alone", "stuck", "w_site")} for sd in SEEDS[:4]: r = M.run(M.vary(M.BASE, n=n, seed=sd, **kw), keep=True) S = r["_s"] parts["log10_lam"].append(np.log10(np.maximum(S["lam"], 1e-300)).astype(np.float32)) parts["fast"].append(S["fast"].astype(np.bool_)); parts["alone"].append((S["N_ever"] < 1).astype(np.bool_)) parts["stuck"].append((S["cls"] == 3).astype(np.bool_)); parts["w_site"].append(S["weights"]["site"].astype(np.float32)) for k, v in parts.items(): carr[f"{cname}__{k}"] = np.concatenate(v) w = carr[f"{cname}__w_site"].astype(np.float64) camp_meta[cname] = dict(overrides={k: (list(v) if isinstance(v, tuple) else v) for k, v in kw.items()}, p_alone=float(w[carr[f"{cname}__alone"]].sum() / w.sum()), p_fast=float(w[carr[f"{cname}__fast"]].sum() / w.sum())) print(f"camp {cname}: P(alone) {camp_meta[cname]['p_alone']:.4f} P(fast) {camp_meta[cname]['p_fast']:.4f}", flush=True) np.savez_compressed(os.path.join(DATA, "model_samples_camps.npz"), **carr) meta = dict(model_md5=md5, model_version="4.1 (file models/model4.py)", n_per_run=n, seeds=SEEDS, classes=M.CLASSES, schools=SCHOOLS, camps=camp_meta, sg_fp=M.BASE["sg_fp"], sg_frac=M.BASE["sg_frac"], w_fast=M.BASE["w_fast"], fast_decades=M.BASE["fast_decades"], slow_decades=list(M.BASE["slow_decades"]), T_WIN_Gyr=M.T_WIN, per_seed=check) json.dump(meta, open(os.path.join(DATA, "model_samples_meta.json"), "w"), indent=1) print("saved", len(arr["log10_lam"]), "samples; model md5", md5) if __name__ == "__main__": main()