#!/usr/bin/env python3 """Single entry point: regenerates every table, number and figure of the mdwarf_early study. python3 code/run_all.py # uses cached raw downloads if present, else downloads python3 code/run_all.py --refresh # force re-download of NASA Exoplanet Archive + arXiv metadata python3 code/run_all.py --offline # no network: use the snapshots stored in data/ Runtime: about 3.5 minutes, single-threaded. Outputs: data/*.csv, data/results_summary.json, figures/*.png. (Revised after the referee report; see report.md, "Response to referee".) """ import json import os import sys import time import numpy as np import pandas as pd sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import build_table # noqa: E402 import earliness as E # noqa: E402 import fetch_data # noqa: E402 import figures # noqa: E402 import inference as inf # noqa: E402 from common import DATA_DIR, SEED, M_DWARF_MAX_MASS, F_ATM_FLOOR, A_EARTH_GYR # noqa: E402 T_START = time.time() STALE = [("data", "posterior_earliness.csv"), ("data", "posterior_earliness_sfh_sensitivity.csv"), ("figures", "fig5_posterior_early.png")] # first-pass outputs that were renamed def log(msg): print(f"[{time.time() - T_START:6.1f}s] {msg}", flush=True) def pct(x, qs=(5, 16, 50, 84, 95)): return [float(v) for v in np.percentile(x, qs)] def main(): refresh = "--refresh" in sys.argv offline = "--offline" in sys.argv R = {} # everything quoted in report.md goes in here for sub_dir, f in STALE: path = os.path.join(os.path.dirname(DATA_DIR), sub_dir, f) if os.path.exists(path): os.remove(path) # ------------------------------------------------------------------ Part 1 sel, small, meta = fetch_data.main(refresh=refresh, offline=offline) census, targets = build_table.build(sel, meta) log(f"Part 1: census of {len(census)} planets; {len(small)} small M-dwarf planets in archive sample") c = census m_rocky = (c.host_class == "M") & (c.group == "rocky") strict = m_rocky & c.outcome.isin(inf.INFORMATIVE) R["census"] = dict( n_rows=int(len(c)), n_Mdwarf=int((c.host_class == "M").sum()), n_Mdwarf_rocky=int(m_rocky.sum()), n_Mdwarf_volatile_rich=int(((c.host_class == "M") & (c.group == "volatile_rich?")).sum()), n_other_host_rocky_or_volatile=int(((c.host_class != "M") & (c.group != "sub_neptune")).sum()), n_sub_neptune=int((c.group == "sub_neptune").sum()), outcome_counts_Mdwarf_rocky=c[m_rocky].outcome.value_counts().to_dict(), outcome_counts_all=c.outcome.value_counts().to_dict(), n_decisive=int(strict.sum()), decisive_planets=list(c.planet[strict]), n_decisive_hosts=len({inf.host_of(p) for p in c.planet[strict]}), n_decisive_S_below_25=int((strict & ~c.lava_world & (c.S_earth < 25)).sum()), n_decisive_in_HZ=int((strict & (c.S_earth < 0.9)).sum()), min_S_decisive=float(c.S_earth[strict].min()), median_S_decisive=float(c.S_earth[strict].median()), n_citations=int(len(fetch_data.all_arxiv_ids())), nea_rows_selected=int(len(sel)), nea_rows_small_mdwarf=int(len(small)), msini_targets=list(targets.pl_name[targets.mass_is_lower_limit])) # ------------------------------------------------------------------ Part 2a bb = inf.beta_binomial_table(census) def _fkey(r): k = r["sample"].split(".")[0] if "(pi_" in r["sample"]: k += "_" + r["sample"].split("(")[1].split(")")[0] return k + "|" + r["prior"] R["fraction"] = {_fkey(r): dict( k=None if pd.isna(r["k"]) else int(r["k"]), n=int(r["n"]), mean=r["mean"], median=r["median"], upper95=r["upper95"]) for _, r in bb.iterrows()} log("Part 2a: atmosphere-fraction posteriors written") # ------------------------------------------------------------------ Part 2b coef, n_hosts, ml_scatter = inf.fit_mass_luminosity(small) R["mass_luminosity"] = dict(coef=[float(x) for x in coef], n_hosts=n_hosts, rms_dex=ml_scatter) fits = {} fit_specs = [ ("base", "base", inf.PRIOR_WEAK, "logit"), ("plus_subneptunes", "plus_subneptunes", inf.PRIOR_WEAK, "logit"), ("flat_slope_priors", "base", inf.PRIOR_FLAT, "logit"), ("probit_link", "base", inf.PRIOR_WEAK, "probit"), ("plus_lava", "plus_lava", inf.PRIOR_WEAK, "logit"), ("ss_no_thin", "ss_no_thin", inf.PRIOR_WEAK, "logit"), ("pi_pess", "pi_pess", inf.PRIOR_WEAK, "logit"), ("pi_opt", "pi_opt", inf.PRIOR_WEAK, "logit"), ("rocky_only", "rocky_only", inf.PRIOR_WEAK, "logit"), ("drop_trappist", "drop_trappist", inf.PRIOR_WEAK, "logit"), ("drop_ltt1445", "drop_ltt1445", inf.PRIOR_WEAK, "logit"), ("solar_system_only", "solar_system_only", inf.PRIOR_WEAK, "logit"), ] imf_w = inf.imf_weights_kroupa(inf.MASS_GRID_M) tnames = list(targets.pl_name) i_t1e, i_lhs = tnames.index("TRAPPIST-1 e"), tnames.index("LHS 1140 b") rows_fit, pred_cols = [], {} for k, (label, variant, prior, link) in enumerate(fit_specs): d = inf.make_dataset(census, variant) s, rhat = inf.run_mcmc(d, prior=prior, link=link, seed=SEED + k) f = inf.f_atm_population(s, coef, link) P = inf.predict_targets(s, targets, link) fits[label] = dict(samples=s, link=link, f_atm=f, P=P, data=d) sm = inf.summarize_params(s) favg = (f * imf_w).sum(axis=1) rows_fit.append(dict( fit=label, n_points=len(d["pi"]), rhat_max=round(rhat, 4), a=sm["a"][0], a_lo=sm["a"][1], a_hi=sm["a"][2], p=sm["p"][0], p_lo=sm["p"][1], p_hi=sm["p"][2], q=sm["q"][0], q_lo=sm["q"][1], q_hi=sm["q"][2], q_p05=float(np.percentile(s[:, 2], 5)), q_p95=float(np.percentile(s[:, 2], 95)), w_dex=sm["w_dex"][0], w_lo=sm["w_dex"][1], w_hi=sm["w_dex"][2], f_atm_HZ_imfavg_mean=favg.mean(), f_atm_HZ_imfavg_p05=np.percentile(favg, 5), f_atm_HZ_imfavg_median=np.median(favg), f_atm_HZ_imfavg_p95=np.percentile(favg, 95), f_atm_HZ_at_0p1Msun_mean=float(np.interp(0.1, inf.MASS_GRID_M, f.mean(axis=0))), f_atm_HZ_at_0p3Msun_mean=float(np.interp(0.3, inf.MASS_GRID_M, f.mean(axis=0))), f_atm_HZ_at_0p5Msun_mean=float(np.interp(0.5, inf.MASS_GRID_M, f.mean(axis=0))), P_TRAPPIST1e_mean=float(P[:, i_t1e].mean()), P_LHS1140b_mean=float(P[:, i_lhs].mean()))) pred_cols[label] = P log(f" shoreline fit '{label}': N={len(d['pi'])}, R-hat={rhat:.3f}, p={sm['p'][0]:.2f}, q={sm['q'][0]:.2f}, " f"w={sm['w_dex'][0]:.2f} dex, ={favg.mean():.3f}") fit_tab = pd.DataFrame(rows_fit) fit_tab.to_csv(os.path.join(DATA_DIR, "shoreline_fits.csv"), index=False) # how much of q and of the habitable-zone extrapolation is the prior? (referee point 4) d_base = inf.make_dataset(census, "base") qrows = [] for k, (label, prior) in enumerate(inf.q_prior_variants().items()): if label.endswith("[baseline]"): s, rhat, f, P = fits["base"]["samples"], rows_fit[0]["rhat_max"], fits["base"]["f_atm"], fits["base"]["P"] else: s, rhat = inf.run_mcmc(d_base, prior=prior, seed=SEED + 300 + k) f, P = inf.f_atm_population(s, coef), inf.predict_targets(s, targets) q = s[:, 2] qrows.append(dict(q_prior=label, rhat_max=round(float(rhat), 4), q_median=np.median(q), q_p05=np.percentile(q, 5), q_p16=np.percentile(q, 16), q_p84=np.percentile(q, 84), q_p95=np.percentile(q, 95), p_median=np.median(s[:, 1]), f_atm_HZ_imfavg_mean=float((f * imf_w).sum(axis=1).mean()), P_TRAPPIST1e_mean=float(P[:, i_t1e].mean()), P_LHS1140b_mean=float(P[:, i_lhs].mean()))) sflat = fits["flat_slope_priors"]["samples"][:, 2] qrows.append(dict(q_prior="flat on [-2,5] (and flat p)", rhat_max=rows_fit[2]["rhat_max"], q_median=np.median(sflat), q_p05=np.percentile(sflat, 5), q_p16=np.percentile(sflat, 16), q_p84=np.percentile(sflat, 84), q_p95=np.percentile(sflat, 95), p_median=np.median(fits["flat_slope_priors"]["samples"][:, 1]), f_atm_HZ_imfavg_mean=rows_fit[2]["f_atm_HZ_imfavg_mean"], P_TRAPPIST1e_mean=rows_fit[2]["P_TRAPPIST1e_mean"], P_LHS1140b_mean=rows_fit[2]["P_LHS1140b_mean"])) qtab = pd.DataFrame(qrows) qtab.to_csv(os.path.join(DATA_DIR, "shoreline_q_prior_sensitivity.csv"), index=False) R["q_prior_sensitivity"] = qtab.to_dict(orient="records") log(" q-prior sensitivity: q median %s; HZ retention %s" % ( ", ".join(f"{v:.2f}" for v in qtab.q_median), ", ".join(f"{v:.2f}" for v in qtab.f_atm_HZ_imfavg_mean))) # "reference mixture": equal mixture of three structural choices. It brackets the defensible # range; it is NOT an estimate (the three components give very different answers). POOL = ["base", "plus_subneptunes", "flat_slope_priors"] rng = np.random.default_rng(SEED) n_pool = 6000 idx = {k: rng.choice(len(fits[k]["samples"]), n_pool // len(POOL), replace=False) for k in POOL} pooled_f = np.vstack([fits[k]["f_atm"][idx[k]] for k in POOL]) pooled_P = np.vstack([fits[k]["P"][idx[k]] for k in POOL]) fits["pooled"] = dict(f_atm=pooled_f, P=pooled_P) pred_rows = [] for i, nme in enumerate(tnames): row = dict(planet=nme, S_earth=round(float(targets.S_earth[i]), 3), v_esc_kms=round(float(targets.v_esc_kms[i]), 2), host_logL=round(float(targets.st_lum[i]), 3), mass_source=targets.mass_source[i], mass_is_lower_limit=bool(targets.mass_is_lower_limit[i])) for k in ["base", "plus_subneptunes", "flat_slope_priors", "pooled", "probit_link", "pi_pess", "pi_opt", "rocky_only", "solar_system_only"]: Pi = fits[k]["P"][:, i] row[f"{k}_mean"] = Pi.mean() if k in ("base", "plus_subneptunes", "flat_slope_priors", "pooled", "solar_system_only"): row[f"{k}_p05"], row[f"{k}_p50"], row[f"{k}_p95"] = np.percentile(Pi, [5, 50, 95]) pred_rows.append(row) pred = pd.DataFrame(pred_rows) pred.to_csv(os.path.join(DATA_DIR, "shoreline_predictions_temperate.csv"), index=False) fm = pd.DataFrame({"stellar_mass_Msun": inf.MASS_GRID_M, "logL_Lsun": inf.logL_of_mass(inf.MASS_GRID_M, coef)}) for k in POOL + ["pooled", "solar_system_only"]: f = fits[k]["f_atm"] fm[f"{k}_mean"] = f.mean(axis=0) fm[f"{k}_p05"], fm[f"{k}_p50"], fm[f"{k}_p95"] = np.percentile(f, [5, 50, 95], axis=0) fm.to_csv(os.path.join(DATA_DIR, "f_atm_vs_stellar_mass.csv"), index=False) R["shoreline"] = {r["fit"]: {k: (float(v) if isinstance(v, (float, np.floating)) else v) for k, v in r.items() if k != "fit"} for r in rows_fit} pf = (pooled_f * imf_w).sum(axis=1) i01 = int(np.argmin(abs(inf.MASS_GRID_M - 0.1))) R["shoreline"]["pooled"] = dict(definition="equal-weight mixture of fits: " + ", ".join(POOL), n_samples=n_pool, f_atm_HZ_imfavg_mean=float(pf.mean()), f_atm_HZ_imfavg_pct=pct(pf), f_atm_HZ_at_0p1Msun_mean=float(np.interp(0.1, inf.MASS_GRID_M, pooled_f.mean(axis=0))), f_atm_HZ_at_0p3Msun_mean=float(np.interp(0.3, inf.MASS_GRID_M, pooled_f.mean(axis=0))), f_atm_HZ_at_0p5Msun_mean=float(np.interp(0.5, inf.MASS_GRID_M, pooled_f.mean(axis=0))), frac_draws_f_at_0p1Msun_below_1e3=float((pooled_f[:, i01] < 1e-3).mean()), frac_draws_f_at_0p1Msun_below_1e10=float((pooled_f[:, i01] < 1e-10).mean()), f_at_0p1Msun_pct=pct(pooled_f[:, i01])) R["predictions"] = {r["planet"]: {k: float(r[k]) for k in r.index if k.endswith(("_mean", "_p05", "_p50", "_p95"))} for _, r in pred.iterrows()} log("Part 2b: shoreline fits and predictions written") # ------------------------------------------------------------------ Part 2c: p_M def p_m_matrix(f_atm, seed, penalties=True, floor=F_ATM_FLOOR): """M-dwarf weight samples on the M-dwarf mass grid: floored retention [x f_water x f_other].""" f = inf.apply_floor(f_atm, floor) if not penalties: return f r = np.random.default_rng(seed) fw, fo = inf.draw_p_M_factors(len(f), r) return f * (fw * fo)[:, None] f_ss = fits["solar_system_only"]["f_atm"][rng.choice(len(fits["solar_system_only"]["f_atm"]), n_pool, replace=False)] pM_jwst = p_m_matrix(pooled_f, SEED + 101) pM_pre = p_m_matrix(f_ss, SEED + 102) pm_avg = (pM_jwst * imf_w).sum(axis=1) pm_avg_pre = (pM_pre * imf_w).sum(axis=1) pm_avg_ret = (inf.apply_floor(pooled_f) * imf_w).sum(axis=1) fw_, fo_ = inf.draw_p_M_factors(200000, np.random.default_rng(SEED + 103)) R["p_M"] = dict( floor=F_ATM_FLOOR, jwst_imfavg_pct=pct(pm_avg), jwst_imfavg_mean=float(pm_avg.mean()), jwst_retention_only_imfavg_pct=pct(pm_avg_ret), jwst_retention_only_imfavg_mean=float(pm_avg_ret.mean()), jwst_at_0p1Msun_pct=pct(pM_jwst[:, i01]), jwst_at_0p1Msun_mean=float(pM_jwst[:, i01].mean()), pre_imfavg_pct=pct(pm_avg_pre), pre_imfavg_mean=float(pm_avg_pre.mean()), penalties_only_pct=pct(fw_ * fo_), penalties_only_mean=float((fw_ * fo_).mean()), factors=f"f_water ~ logU{inf.F_WATER_RANGE}; f_other ~ logU{inf.F_OTHER_RANGE} (both assumed); " f"retention from the reference-mixture shoreline with floor {F_ATM_FLOOR}") log(f"Part 2c: p_M (IMF-averaged) median {np.median(pm_avg):.3f}, 90% [{np.percentile(pm_avg, 5):.4f}, {np.percentile(pm_avg, 95):.3f}]") # ------------------------------------------------------------------ Part 3 # 3a. replication of Hanson et al. 2021 checks = [] for (n, L, quoted) in [(6, 10, "< 10.6% (all panels)"), (3, 15, "< 10% (all panels)"), (2, 20, "< 8.4% (kappa=0); < 12.4% (kappa=3)")]: for kappa in (0, 3): vals = [E.hanson_rank(n, L, phi=p, kappa=kappa) for p in (2, 4, 8)] checks.append(dict(n=n, Lbar_Gyr=L, kappa=kappa, rank_phi2=vals[0], rank_phi4=vals[1], rank_phi8=vals[2], max_rank=max(vals), hanson_quote=quoted)) hc = pd.DataFrame(checks) hc.to_csv(os.path.join(DATA_DIR, "hanson_replication_checks.csv"), index=False) R["hanson_checks"] = hc.to_dict(orient="records") # Hanson eq. 1 with a 4 Gyr SFR peak (chi = 4) at Lbar = 45 Gyr: the K-dwarf result in their framework R["hanson_eq1_L45_chi4"] = {str(n): float(E.hanson_rank(n, 45.0, phi=4.0, chi=4.0)) for n in (1, 2, 3, 4, 6)} R["hanson_eq1_L10_chi4"] = {str(n): float(E.hanson_rank(n, 10.0, phi=4.0, chi=4.0)) for n in (1, 3, 6, 10)} n_grid_h = np.array([1, 1.5, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20]) L_grid = np.exp(np.linspace(np.log(5), np.log(2e4), 34)) hanson_grid = np.array([[E.hanson_rank(n, L) for L in L_grid] for n in n_grid_h]) hw_fig = E.host_weight(M_DWARF_MAX_MASS) fig_all = np.empty((len(n_grid_h), len(L_grid))); fig_noM = np.empty_like(fig_all) for i, n in enumerate(n_grid_h): for j, L in enumerate(L_grid): num, den = E.birth_integrals(float(n), L, "MD14", "Kroupa") fig_all[i, j] = num.sum() / den.sum() fig_noM[i, j] = (num * hw_fig).sum() / (den * hw_fig).sum() log("Part 3a: Hanson et al. replication done") # 3b. our model on a grid (n, Lmax); per-mass numerators/denominators are stored for re-weighting n_int = np.arange(1, 11) def grids(model="MD14", imf="Kroupa", fwin=1.0): N_ = np.empty((len(n_int), len(L_grid), len(E.MASS_GRID))); D_ = np.empty_like(N_) for i, n in enumerate(n_int): for j, L in enumerate(L_grid): N_[i, j], D_[i, j] = E.birth_integrals(int(n), L, model, imf, fwin) return N_, D_ NUM, DEN = grids() NUM6, DEN6 = grids(fwin=0.6) hw_m = E.host_weight(M_DWARF_MAX_MASS) # 1 above 0.6 Msun, 0 below, fractional in the boundary cell w_all = np.ones(len(E.MASS_GRID)) CUTS = (0.6, 0.65, 0.7, 0.75, 0.8) w_cut = {cut: E.host_weight(cut) for cut in CUTS} w_noM, w_07 = w_cut[0.6], w_cut[0.7] rank_all = E.rank_from_weights(NUM, DEN, w_all) rank_noM = E.rank_from_weights(NUM, DEN, w_noM) rank_07 = E.rank_from_weights(NUM, DEN, w_07) grid_rows = [] for i, n in enumerate(n_int): for j, L in enumerate(L_grid): row = dict(n=int(n), Lmax_Gyr=L, rank_all_stars=rank_all[i, j]) for cut in CUTS: row[f"rank_hosts_above_{cut:g}".replace(".", "p")] = float(E.rank_from_weights(NUM[i, j], DEN[i, j], w_cut[cut])) grid_rows.append(row) pd.DataFrame(grid_rows).to_csv(os.path.join(DATA_DIR, "rank_grid_n_Lmax.csv"), index=False) # priors over the two nuisance parameters of the earliness model (they carry no data) n_priors = { "uniform n=1..10": np.ones(10) / 10.0, "Hanson-like (lognormal, median 6, sigma 0.35)": None, "n=1 (one hard step)": np.eye(10)[0], "n=1..3 (few steps)": np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0]) / 3.0, } ln = np.exp(-0.5 * ((np.log(n_int) - np.log(6.0)) / 0.35) ** 2) / n_int n_priors["Hanson-like (lognormal, median 6, sigma 0.35)"] = ln / ln.sum() L_priors = { "log-uniform 5 Gyr - 2e4 Gyr": np.ones(len(L_grid)) / len(L_grid), "log-uniform 5 - 100 Gyr (planet-limited)": (L_grid <= 100.0) / (L_grid <= 100.0).sum(), "Lmax >= 1000 Gyr (star-limited)": (L_grid >= 1000.0) / (L_grid >= 1000.0).sum(), } LP0 = "log-uniform 5 Gyr - 2e4 Gyr" R["Lmax_prior_fraction_above_20Gyr"] = float((L_grid > 20.0).mean()) R["Lmax_prior_fraction_above_100Gyr"] = float((L_grid > 100.0).mean()) # headline table: rank for selected (n, Lmax) and the critical p_M. # The critical p_M is reported only where it is not a knife-edge: if the rank without any M dwarf # is itself within a factor 2 of the threshold, the number is an artefact of the exact cut # (referee point 2: at n = 4 the M-dwarf-free rank is 1.05%, 0.05 points above the 1% line). tab = [] for n in (1, 2, 3, 4, 5, 6, 10): for L in (10.0, 20.0, 50.0, 1e3, 1e9): num, den = E.birth_integrals(n, L, "MD14", "Kroupa") row = dict(n=n, Lmax_Gyr=("inf" if L > 1e8 else L), rank_all_stars=float(E.rank_from_weights(num, den, w_all)), rank_hosts_above_0p6=float(E.rank_from_weights(num, den, w_noM)), rank_hosts_above_0p7=float(E.rank_from_weights(num, den, w_07))) for target, lab in ((0.01, "1pct"), (0.10, "10pct")): pc, r0 = E.critical_p_m(n, L, target) knife = (0.0 < pc < 1.0) and (r0 < 2.0 * target) row[f"p_M_crit_rank{lab}"] = np.nan if knife else pc row[f"p_M_crit_rank{lab}_note"] = ( "not early even with p_M = 1" if pc >= 1.0 else "early even with p_M = 0" if pc <= 0.0 else f"knife-edge: rank without M dwarfs is {r0:.4f}, within 2x of the threshold" if knife else "") tab.append(row) rank_tab = pd.DataFrame(tab) rank_tab.to_csv(os.path.join(DATA_DIR, "rank_table.csv"), index=False) R["rank_table"] = rank_tab.to_dict(orient="records") # robustness of the critical weight: suppress all hosts below m_cut by a common factor p, with # either the whole or 60% of the main-sequence lifetime habitable (unlimited Lmax) crit = [] for fwin in (1.0, 0.6): for cut in (0.5, 0.6, 0.7, 0.8): row = dict(f_window=fwin, suppressed_below_Msun=cut) for n in (1, 2, 3): pc, r0 = E.critical_p_m(n, 1e9, 0.01, m_cut=cut, f_window=fwin) row[f"p_crit_rank1pct_n{n}"] = pc row[f"rank_if_fully_excluded_n{n}"] = r0 crit.append(row) crit = pd.DataFrame(crit) crit.to_csv(os.path.join(DATA_DIR, "critical_pM_cut_window.csv"), index=False) R["critical_pM_cut_window"] = crit.to_dict(orient="records") # rank without M dwarfs as a function of the host-mass cut and the habitable-window fraction # (referee point 2). Lmax unlimited, Madau-Dickinson SFH, Kroupa IMF. cw = [] wu0 = n_priors["uniform n=1..10"][:, None] * L_priors[LP0][None, :] for fwin, (N_, D_) in ((1.0, (NUM, DEN)), (0.6, (NUM6, DEN6))): for cut in CUTS: ranks = [E.rank_host_cut(int(n), 1e9, cut, f_window=fwin) for n in n_int] rk = E.rank_from_weights(N_, D_, w_cut[cut]) below1 = [int(n) for n, r in zip(n_int, ranks) if r < 0.01] below10 = [int(n) for n, r in zip(n_int, ranks) if r < 0.10] row = dict(f_window=fwin, host_mass_cut_Msun=cut, habitable_lifetime_at_cut_Gyr=float(fwin * E.ms_lifetime_gyr(cut))) row.update({f"rank_n{int(n)}": r for n, r in zip(n_int, ranks)}) row.update(smallest_n_with_rank_below_1pct=(below1[0] if below1 else np.nan), smallest_n_with_rank_below_10pct=(below10[0] if below10 else np.nan), frac_grid_rank_below_1pct_uniform_n=float(((rk < 0.01) * wu0).sum()), frac_grid_rank_below_10pct_uniform_n=float(((rk < 0.10) * wu0).sum())) cw.append(row) cw = pd.DataFrame(cw) cw.to_csv(os.path.join(DATA_DIR, "rank_noM_cut_window.csv"), index=False) R["rank_noM_cut_window"] = cw.to_dict(orient="records") # where the weight comes from when M dwarfs are excluded at 0.6 Msun bins = ((0.6, 0.65), (0.65, 0.7), (0.7, 0.8), (0.8, 1.0), (1.0, 3.01)) shares = [] for n in n_int: num, den = E.birth_integrals(int(n), 1e9) tot_d, tot_n = (den * hw_m).sum(), (num * hw_m).sum() row = dict(n=int(n)) for a, b in bins: selm = E.host_weight(a) - E.host_weight(b) row[f"share_all_time_{a:g}-{min(b, 3.0):g}"] = float((den * selm).sum() / tot_d) for a, b in bins: selm = E.host_weight(a) - E.host_weight(b) row[f"share_born_by_today_{a:g}-{min(b, 3.0):g}"] = float((num * selm).sum() / tot_n) shares.append(row) shares = pd.DataFrame(shares) shares.to_csv(os.path.join(DATA_DIR, "noM_weight_shares.csv"), index=False) R["noM_weight_shares"] = shares.to_dict(orient="records") # robustness of the M-dwarf-free rank to SFH, IMF and habitable-window fraction rob = [] for model in ("MD14", "MD14_exp", "Sobral13", "BP15_long", "Hanson"): for imf in ("Kroupa", "Chabrier"): for fwin in (1.0, 0.6): for n in (1, 3, 6): for L in (10.0, 1e9): num, den = E.birth_integrals(n, L, model, imf, fwin) rob.append(dict(sfh=model, imf=imf, f_window=fwin, n=n, Lmax_Gyr=("inf" if L > 1e8 else L), rank_all_stars=float(E.rank_from_weights(num, den, w_all)), rank_hosts_above_0p6=float(E.rank_from_weights(num, den, w_noM)))) rob = pd.DataFrame(rob) rob.to_csv(os.path.join(DATA_DIR, "rank_robustness.csv"), index=False) R["sfh_future_over_past"] = {m: float((E._time_grids(m)[3] - E._time_grids(m)[2]) / E._time_grids(m)[2]) for m in ("MD14", "MD14_exp", "Sobral13", "BP15_long", "Hanson")} log("Part 3b: rank grids, host-cut table and robustness tables written") # 3c. prior-weighted fraction of (n, Lmax) space in which the rank is below 1% / 10%. # NOT a posterior: n and Lmax carry no data, only the M-dwarf weight does (referee point 3). def expand_weights(pM): """(Ns, n_Mgrid) p_M(M*) on the inference grid -> (Ns, n_mass) weights on the earliness grid.""" w = np.empty((len(pM), len(E.MASS_GRID))) logm = np.log(inf.MASS_GRID_M) for k in range(len(pM)): w[k] = hw_m + (1.0 - hw_m) * np.interp(np.log(E.MASS_GRID), logm, pM[k]) return w NS = 1500 sub = lambda a, n=NS: a[:: max(1, len(a) // n)][:n] # noqa: E731 deterministic thinning ones = np.ones((1, len(inf.MASS_GRID_M))) weights = {} weights["p_M = 1 (all stars equally habitable)"] = w_all[None, :] weights["penalties only (retention = 1, as in evolution models)"] = expand_weights(p_m_matrix(np.repeat(ones, NS, axis=0), SEED + 104, floor=0.0)) weights["Solar-System-only shoreline, retention only"] = expand_weights(p_m_matrix(sub(f_ss), 0, penalties=False)) weights["Solar-System-only shoreline x penalties (pre-JWST comparator)"] = expand_weights(sub(pM_pre)) weights["JWST shoreline (reference mixture), retention only"] = expand_weights(p_m_matrix(sub(pooled_f), 0, penalties=False)) for k in POOL: weights[f"JWST shoreline fit '{k}', retention only"] = expand_weights( p_m_matrix(sub(fits[k]["f_atm"][idx[k]], 500), 0, penalties=False)) weights["JWST shoreline (reference mixture) x penalties"] = expand_weights(p_m_matrix(sub(pooled_f), SEED + 101)) for fl in (0.0, 1e-2): weights[f"JWST shoreline, retention only, floor {fl:g}"] = expand_weights(p_m_matrix(sub(pooled_f), 0, penalties=False, floor=fl)) weights[f"JWST shoreline x penalties, floor {fl:g}"] = expand_weights(p_m_matrix(sub(pooled_f), SEED + 101, floor=fl)) weights["M dwarfs excluded: p_M = 0 (hosts > 0.6 Msun)"] = w_noM[None, :] weights["hosts > 0.65 Msun only"] = w_cut[0.65][None, :] weights["hosts > 0.7 Msun only"] = w_07[None, :] weights["hosts > 0.8 Msun only"] = w_cut[0.8][None, :] # value of future observations: refit the three reference-mixture shorelines with hypothetical verdicts future = { "future: T-1 e, f, g and LHS 1140 b all airless": {"TRAPPIST-1 e": 0.03, "TRAPPIST-1 f": 0.03, "TRAPPIST-1 g": 0.03, "LHS 1140 b": 0.05}, "future: T-1 e, f and LHS 1140 b have atmospheres": {"TRAPPIST-1 e": 0.95, "TRAPPIST-1 f": 0.95, "LHS 1140 b": 0.95}, } for k, (label, newpi) in enumerate(future.items()): cc = census.copy() for pl, v in newpi.items(): cc.loc[cc.planet == pl, ["pi_base", "outcome"]] = [v, "hypothetical"] fs, by_fit = [], {} for kk, (fitname, variant, prior) in enumerate([("base", "base", inf.PRIOR_WEAK), ("plus_subneptunes", "plus_subneptunes", inf.PRIOR_WEAK), ("flat_slope_priors", "base", inf.PRIOR_FLAT)]): d = inf.make_dataset(cc, variant) s, rhat = inf.run_mcmc(d, prior=prior, seed=SEED + 50 + 10 * k + kk) ff = inf.f_atm_population(s[rng.choice(len(s), 500, replace=False)], coef) fs.append(ff) by_fit[fitname] = float((ff * imf_w).sum(axis=1).mean()) f_fut = np.vstack(fs) fits[label] = dict(f_atm=f_fut) weights[label + ", retention only"] = expand_weights(p_m_matrix(f_fut, 0, penalties=False)) weights[label + " x penalties"] = expand_weights(p_m_matrix(f_fut, SEED + 200 + k)) R.setdefault("future_f_atm", {})[label] = dict( imfavg_mean=float((f_fut * imf_w).sum(axis=1).mean()), imfavg_mean_by_fit=by_fit, at_0p1Msun_mean=float(np.interp(0.1, inf.MASS_GRID_M, f_fut.mean(axis=0)))) log(f" scenario '{label}': = {(f_fut * imf_w).sum(axis=1).mean():.3f}") def ranks_of(W, N_=None, D_=None): """W: (Ns, n_mass) weights -> rank array (Ns, n, L).""" N_ = NUM if N_ is None else N_ D_ = DEN if D_ is None else D_ return np.einsum("sm,nlm->snl", W, N_) / np.einsum("sm,nlm->snl", W, D_) scenarios = {k: ranks_of(W) for k, W in weights.items()} sum_rows, byn_rows = [], [] for sname, rk in scenarios.items(): for lpn, pl in L_priors.items(): f1 = ((rk < 0.01) * pl[None, None, :]).sum(axis=2).mean(axis=0) # by n f10 = ((rk < 0.10) * pl[None, None, :]).sum(axis=2).mean(axis=0) mr = (rk * pl[None, None, :]).sum(axis=2).mean(axis=0) gm = np.exp((np.log(rk) * pl[None, None, :]).sum(axis=2).mean(axis=0)) for i, n in enumerate(n_int): byn_rows.append(dict(scenario=sname, Lmax_prior=lpn, n=int(n), frac_rank_below_1pct=float(f1[i]), frac_rank_below_10pct=float(f10[i]), mean_rank=float(mr[i]), geometric_mean_rank=float(gm[i]))) for npn, pn in n_priors.items(): sum_rows.append(dict(scenario=sname, n_prior=npn, Lmax_prior=lpn, frac_rank_below_1pct=float((f1 * pn).sum()), frac_rank_below_10pct=float((f10 * pn).sum()), mean_rank=float((mr * pn).sum()))) post = pd.DataFrame(sum_rows) post.to_csv(os.path.join(DATA_DIR, "early_fraction_prior_weighted.csv"), index=False) byn = pd.DataFrame(byn_rows) byn.to_csv(os.path.join(DATA_DIR, "early_fraction_by_n.csv"), index=False) R["early_fraction"] = post.to_dict(orient="records") R["early_fraction_by_n_default_Lprior"] = byn[byn.Lmax_prior == LP0].drop(columns="Lmax_prior").to_dict(orient="records") a1 = (rank_all < 0.01) a0 = (rank_noM < 0.01) R["grid_bounds"] = dict( frac_cells_early_pM1=float(a1.mean()), frac_cells_early_pM0=float(a0.mean()), frac_cells_where_M_weight_can_matter=float((a1 & ~a0).mean()), note="prior-weighted fraction with rank<1% (uniform n, default Lmax prior) is confined to " "[frac_cells_early_pM0, frac_cells_early_pM1] whatever the M-dwarf evidence") # 3d. conditioning on our own datum (birth date t0, host mass 1 Msun, planet age 4.54 Gyr) within # the no-deadline model: evidence of each M-dwarf scenario and posterior over (n, Lmax). sd_rows = [] sd_scen = ["p_M = 1 (all stars equally habitable)", "penalties only (retention = 1, as in evolution models)", "JWST shoreline (reference mixture), retention only", "JWST shoreline (reference mixture) x penalties", "M dwarfs excluded: p_M = 0 (hosts > 0.6 Msun)", "hosts > 0.7 Msun only", "hosts > 0.8 Msun only"] Z_ref = {} post_L = {} for sname in sd_scen: like = E.self_datum_likelihood(DEN, weights[sname], n_int, L_grid) # (Ns, n, L) rk = scenarios[sname] for npn in ("uniform n=1..10", "n=1 (one hard step)", "n=1..3 (few steps)"): for lpn in (LP0, "log-uniform 5 - 100 Gyr (planet-limited)"): prior = n_priors[npn][:, None] * L_priors[lpn][None, :] joint = (like * prior[None]).mean(axis=0) # (n, L), averaged over p_M draws Z = float(joint.sum()) pj = joint / Z key = (npn, lpn) if sname == sd_scen[0]: Z_ref[key] = Z pl_marg = pj.sum(axis=0) order = np.argsort(L_grid) cdfL = np.cumsum(pl_marg[order]) L95 = float(L_grid[order][np.searchsorted(cdfL, 0.95)]) prk1 = float(((like * prior[None]) * (rk < 0.01)).mean(axis=0).sum() / Z) mrk = float(((like * prior[None]) * rk).mean(axis=0).sum() / Z) sd_rows.append(dict(scenario=sname, n_prior=npn, Lmax_prior=lpn, evidence_per_Gyr=Z, evidence_ratio_to_pM1=Z / Z_ref[key], post_P_Lmax_le_10Gyr=float(pl_marg[L_grid <= 10.0].sum()), post_P_Lmax_le_20Gyr=float(pl_marg[L_grid <= 20.0].sum()), post_Lmax_95pct_upper_Gyr=L95, post_mean_n=float((pj.sum(axis=1) * n_int).sum()), post_P_rank_below_1pct=prk1, post_mean_rank=mrk)) if npn == "uniform n=1..10" and lpn == LP0: post_L[sname] = pl_marg sd = pd.DataFrame(sd_rows) sd.to_csv(os.path.join(DATA_DIR, "self_datum_likelihood.csv"), index=False) R["self_datum"] = sd.to_dict(orient="records") # likelihood ratio (hosts > 0.6 only : all stars) cell by cell, for the report l1 = E.self_datum_likelihood(DEN, w_all[None, :], n_int, L_grid)[0] l0 = E.self_datum_likelihood(DEN, w_noM[None, :], n_int, L_grid)[0] jl = lambda Lsel: int(np.argmin(abs(np.log(L_grid) - np.log(Lsel)))) # noqa: E731 R["self_datum_cell_ratio_noM_over_all"] = { f"n={n}, Lmax={L_grid[jl(Ls)]:.0f} Gyr": float(l0[n - 1, jl(Ls)] / l1[n - 1, jl(Ls)]) for n in (1, 3, 6) for Ls in (5.0, 10.0, 100.0, 2e4)} R["self_datum_cell_like_rel_to_best_pM1"] = { f"n={n}, Lmax={L_grid[jl(Ls)]:.0f} Gyr": float(l1[n - 1, jl(Ls)] / l1.max()) for n in (1, 3, 6) for Ls in (5.0, 10.0, 20.0, 100.0, 2e4)} log("Part 3c/d: prior-weighted early fractions and self-datum likelihoods written") # sensitivity of the prior-weighted fractions to the star-formation history and IMF W_j = weights["JWST shoreline (reference mixture) x penalties"] W_r = weights["JWST shoreline (reference mixture), retention only"] sens = [] wgt0 = n_priors["uniform n=1..10"][:, None] * L_priors[LP0][None, :] wgt1 = n_priors["n=1 (one hard step)"][:, None] * L_priors[LP0][None, :] wgt3 = n_priors["n=1..3 (few steps)"][:, None] * L_priors[LP0][None, :] for model, imf in (("MD14", "Kroupa"), ("MD14", "Chabrier"), ("MD14_exp", "Kroupa"), ("Sobral13", "Kroupa"), ("BP15_long", "Kroupa"), ("Hanson", "Kroupa")): N2, D2 = (NUM, DEN) if (model, imf) == ("MD14", "Kroupa") else grids(model, imf) for sname, W in (("p_M = 1 (all stars equally habitable)", w_all[None, :]), ("JWST shoreline (reference mixture), retention only", W_r), ("JWST shoreline (reference mixture) x penalties", W_j), ("M dwarfs excluded: p_M = 0 (hosts > 0.6 Msun)", w_noM[None, :])): rk = ranks_of(W, N2, D2) for npn, wgt in (("uniform n=1..10", wgt0), ("n=1 (one hard step)", wgt1), ("n=1..3 (few steps)", wgt3)): sens.append(dict(sfh=model, imf=imf, scenario=sname, n_prior=npn, frac_rank_below_1pct=float(((rk < 0.01) * wgt[None]).sum(axis=(1, 2)).mean()), frac_rank_below_10pct=float(((rk < 0.10) * wgt[None]).sum(axis=(1, 2)).mean()), mean_rank=float((rk * wgt[None]).sum(axis=(1, 2)).mean()))) sens = pd.DataFrame(sens) sens.to_csv(os.path.join(DATA_DIR, "early_fraction_sfh_sensitivity.csv"), index=False) R["early_fraction_sfh_sensitivity"] = sens.to_dict(orient="records") # rank for JWST-informed weights at selected (n, Lmax): median and 90% interval rk_j = scenarios["JWST shoreline (reference mixture) x penalties"] rk_r = scenarios["JWST shoreline (reference mixture), retention only"] rk_p = scenarios["Solar-System-only shoreline x penalties (pre-JWST comparator)"] sel_rows = [] for n in (1, 2, 3, 6): for Lsel in (20.0, 100.0, 2e4): j = jl(Lsel) i = n - 1 sel_rows.append(dict(n=n, Lmax_Gyr=float(L_grid[j]), rank_pM1=float(rank_all[i, j]), rank_preJWST_median=float(np.median(rk_p[:, i, j])), rank_JWST_retention_only_p05=float(np.percentile(rk_r[:, i, j], 5)), rank_JWST_retention_only_median=float(np.median(rk_r[:, i, j])), rank_JWST_retention_only_p95=float(np.percentile(rk_r[:, i, j], 95)), rank_JWST_x_penalties_p05=float(np.percentile(rk_j[:, i, j], 5)), rank_JWST_x_penalties_median=float(np.median(rk_j[:, i, j])), rank_JWST_x_penalties_p95=float(np.percentile(rk_j[:, i, j], 95)), rank_noM=float(rank_noM[i, j]))) selt = pd.DataFrame(sel_rows) selt.to_csv(os.path.join(DATA_DIR, "rank_with_jwst_weights.csv"), index=False) R["rank_with_jwst_weights"] = selt.to_dict(orient="records") log("Part 3: sensitivity tables written") # ------------------------------------------------------------------ figures figures.fig_shoreline(census, targets, fits["base"]["samples"]) figures.fig_predictions(pred, fm) figures.fig_hanson(n_grid_h, L_grid, hanson_grid, fig_all, fig_noM) figures.fig_rank_vs_pm(pM_jwst[:, i01]) figures.fig_early_by_n(cw, byn[byn.Lmax_prior == LP0]) log("figures written") json.dump(R, open(os.path.join(DATA_DIR, "results_summary.json"), "w"), indent=1, default=float) log("done -> data/results_summary.json") if __name__ == "__main__": main()