"""All computations of the impact_window study. Called by run_all.py. Writes every table to data/*.csv and the headline numbers to data/results.json, and returns the in-memory results that figures.py draws. """ from __future__ import annotations import json import os import time import numpy as np import pandas as pd import bayes import chronology as chron import config as cfg import model # Monte Carlo draws per chronology x SFD family: main grid, and each robustness variant. The # environment variables exist only for quick tests; the published numbers use the defaults. N_MAIN = int(os.environ.get("IMPACT_WINDOW_N_MAIN", 20000)) N_ROB = int(os.environ.get("IMPACT_WINDOW_N_ROB", 2000)) D_FULL = np.array([10., 20., 40., 150., 190., 250., 300., 350., 440., 500., 600., 740., 1000., 1300., 2000., 2700.]) D_TABLE = [190., 300., 440., 500., 740., 1000., 2000.] # thresholds reported in tables D_BF = [190., 300., 440., 740., 1000., 2000.] # thresholds used as "sterilising" T_TABLE = [4.4, 4.3, 4.2, 4.1, 4.0, 3.9, 3.8, 3.7, 3.48] M_LIST = [0, 1, 2, 3, 4, 5] T_W = cfg.T0_ZIRCON F = cfg.FUTURE_BASE CONDS = list(cfg.HSE_CONDS) # every treatment of the HSE mass constraint MAIN = list(cfg.HSE_MAIN) # the two carried through every table: "upper" and "two" # like-for-like pairs for the effect of the Chang'e-6 re-dating of SPA CE6_PAIRS = [("neukum1983", "yue2026", "Neukum 1983 -> Yue+ 2026 refit"), ("anchor_joy433", "anchor_ce6", "anchors: SPA 4.33 Ga -> 4.247 Ga")] def _save(df, name): path = os.path.join(cfg.DATA, name) df.to_csv(path, index=False) return path def _dj(res, D): return int(np.argmin(np.abs(res["D"] - D))) # --------------------------------------------------------------------------------------------- # Static tables # --------------------------------------------------------------------------------------------- def basin_tables(B): """Basin table with the age of every dated basin under each chronology (central densities).""" nk = chron.MODELS["neukum1983"] out = B[["name", "D_main_km", "main_confidence", "d_imp_km", "period", "N20", "N20_err", "age_neukum_Ga", "err_plus", "err_minus", "dated"]].copy() n1 = nk.n1(B.age_neukum_Ga.to_numpy()) n1_spa = n1[B.name == "South Pole-Aitken"][0] n1 = np.minimum(n1, n1_spa) for k in chron.FUNCTIONAL + ["sawtooth_aspub", "sawtooth_uniform"]: out[f"age_{k}"] = np.where(B.dated, chron.MODELS[k].age(n1), np.nan) n1_imb = n1[B.name == "Imbrium"][0] n1_nec = n1[B.name == "Nectaris"][0] for k, (t_spa, _) in (("anchor_ce6", cfg.SPA_CE6), ("anchor_joy433", cfg.SPA_JOY)): t_nec = 0.5 * (3.95 + min(4.20, t_spa - 0.03)) a = np.full(len(B), np.nan) a[B.dated.to_numpy()] = chron.anchored_ages(n1[B.dated.to_numpy()], n1_imb, n1_nec, n1_spa, cfg.IMBRIUM_AGE[0], t_nec, t_spa) out[f"age_{k}"] = a _save(out.round(4), "basin_table.csv") _save(chron.summary_table().round(4), "chronology_table.csv") return out def ratio_table(): rows = [] for v in (11.0, 14.0, 16.5, 19.9, 22.0): for a in (60.3, 40.0): R = float(model.earth_moon_ratio(v, a)) rows.append(dict(v_moon_kms=v, a_moon_RE=a, earth_over_moon=R, per_unit_area=R / (cfg.R_EARTH / cfg.R_MOON) ** 2)) df = pd.DataFrame(rows) _save(df.round(3), "earth_moon_ratio.csv") rng = np.random.default_rng(cfg.SEED) R = model.earth_moon_ratio(rng.uniform(11, 22, 200000), rng.uniform(*cfg.A_MOON_RE, 200000)) return df, dict(median=float(np.median(R)), p05=float(np.quantile(R, 0.05)), p95=float(np.quantile(R, 0.95)), min=float(R.min()), max=float(R.max())) def scaling_table(): a, b = model.scaling_fit() rows = [dict(kind="calibration point", D_basin_km=D, d_impactor_km=d, d_fit_km=a * D ** b, source=src) for D, d, src in model.SCALING_POINTS] for name, d in cfg.THRESH.items(): rows.append(dict(kind="threshold mapped to the Moon", D_basin_km=(d / a) ** (1 / b), d_impactor_km=d, d_fit_km=d, source=name)) df = pd.DataFrame(rows) _save(df.round(2), "impactor_scaling.csv") resid = np.log([r["d_impactor_km"] / r["d_fit_km"] for r in rows[:len(model.SCALING_POINTS)]]) return dict(a=a, b=b, rms_ln=float(np.sqrt(np.mean(resid ** 2)))) def sfd_tables(B): """Power-law fits to the lunar basin diameters and the test of the main-belt shape.""" from scipy import stats a, b = model.scaling_fit() rows = [] for xmin in (200.0, 300.0, 400.0, 500.0, 600.0, 800.0): r = model.pareto_fit(B.D_main_km.to_numpy(), xmin) rows.append(dict(model="power law", D_min_km=xmin, n=r["n"], q_basin=r["q"], q_basin_err=r["q_err"], q_impactor=r["q"] / b, KS_D=r["ks"], p_boot=r["p"])) mb = model.MainBeltSFD() for xmin in (300.0, 600.0): d = a * B.D_main_km[B.D_main_km >= xmin].to_numpy() ** b dref = a * xmin ** b ks = stats.kstest(d, lambda v: 1.0 - mb.N(v) / mb.N(dref)) rows.append(dict(model="main-belt shape", D_min_km=xmin, n=len(d), q_basin=np.nan, q_basin_err=np.nan, q_impactor=np.nan, KS_D=ks.statistic, p_boot=ks.pvalue)) fits = pd.DataFrame(rows) _save(fits.round(4), "sfd_fits.csv") # predicted numbers of large lunar basins under each family, normalised to its reference count n300 = int((B.D_main_km >= 300).sum()) n600 = int((B.D_main_km >= 600).sum()) q300 = fits[(fits.model == "power law") & (fits.D_min_km == 300)].q_basin.iloc[0] q600 = fits[(fits.model == "power law") & (fits.D_min_km == 600)].q_basin.iloc[0] rows = [] for D in (600.0, 1000.0, 1321.0, 2400.0): d = a * D ** b rows.append(dict(D_basin_km=D, d_impactor_km=d, observed=int((B.D_main_km >= D).sum()), mainbelt=n300 * float(mb.N(d) / mb.N(a * 300 ** b)), lunar_pl300=n300 * (D / 300.0) ** (-q300), lunar_pl600=n600 * (D / 600.0) ** (-q600) if D >= 600 else np.nan)) pred = pd.DataFrame(rows) for fam in ("mainbelt", "lunar_pl300", "lunar_pl600"): # P(as few or fewer basins | prediction) pred[f"poisson_p_le_obs_{fam}"] = stats.poisson.cdf(pred.observed, pred[fam]) _save(pred.round(4), "sfd_predicted_lunar_counts.csv") mbn = {int(d): float(mb.N(d)) for d in (10, 20, 50, 100, 190, 300, 440, 500, 740, 900)} # deterministic cross-check of the two-sided HSE treatment: how many bodies >= 440 km and >= 190 km # come with each 1e22 kg of main-belt-shaped material (bodies >= the 300 km basin impactor) m_tot = float(mb.cm[int((mb.d >= a * 300.0 ** b).sum())]) mbn["mass_ge_dref_kg"] = m_tot mbn["n_ge440_per_1e22kg"] = float(mb.N(440.0)) / m_tot * 1e22 mbn["n_ge190_per_1e22kg"] = float(mb.N(190.0)) / m_tot * 1e22 return fits, pred, mbn # --------------------------------------------------------------------------------------------- # Monte Carlo grid and mixtures # --------------------------------------------------------------------------------------------- def run_grid(B, n_draw, seed, only_sfd=None, chron_swap=None, **kw): """One result per chronology x SFD family. Every combination gets its own seeded generator, so a single combination can be rerun in isolation. chron_swap maps a chronology name of the main set to a replacement (used for the alternative sawtooth normalisations).""" out = {} for i, c in enumerate(chron.ALL): for j, s in enumerate(model.SFD_FAMILIES): if only_sfd and s not in only_sfd: continue rng = np.random.default_rng([seed, i, j]) cc = (chron_swap or {}).get(c, c) out[(c, s)] = model.run_chunked(cc, s, n_draw, rng, basins=B, d_grid=D_FULL, **kw) return out def make_groups(grid, cond): """Mixtures under one HSE treatment: all models, by chronology, by size distribution, all models without the main-belt family, and every single chronology x SFD combination.""" g = {"ALL": model.mixture(list(grid.values()), cond, "ALL")} for c in chron.ALL: sel = [v for (cc, s), v in grid.items() if cc == c] if sel: g[f"chron:{c}"] = model.mixture(sel, cond, f"chron:{c}") for s in model.SFD_FAMILIES: sel = [v for (c, ss), v in grid.items() if ss == s] if sel: g[f"sfd:{s}"] = model.mixture(sel, cond, f"sfd:{s}") sel = [v for (c, s), v in grid.items() if s != "mainbelt"] if sel: g["ALL without main-belt shape"] = model.mixture(sel, cond, "no-mainbelt") for (c, s), r in grid.items(): g[f"{c}|{s}"] = model.view(r, cond) return g def is_combo(name): return "|" in name def model_grid_table(grid): rows = [] for (c, s), r in grid.items(): row = dict(chron=c, sfd=s, n_draw=len(r["mass"]), n_basins=r["n_basins"], t_spa_median=float(np.median(r["t_spa"])), R_median=float(np.median(r["R"])), d_ref_median_km=float(np.median(r["d_ref"])), q_impactor_median=float(np.median(r["q_lunar"])), earth_impactors_ge_dref_median=float(np.median(r["n_ref_earth"])), lunar_ref_basins_after_youngest_mean=float(np.nanmean(r["post_n"])), mass_median_kg=float(np.median(r["mass"])), mass_p16_kg=float(np.quantile(r["mass"], 0.16)), mass_p84_kg=float(np.quantile(r["mass"], 0.84)), expected_mass_median_kg=float(np.median(r["mass_exp"])), max_abs_diff_sampled_vs_analytic=float(np.abs(r["surv_s"] - r["surv"]).max())) for cond in CONDS: row[f"pass_fraction_{cond}"] = r["frac_by"][cond] for cond in MAIN: v = model.view(r, cond) for D in (440., 500., 1000.): row[f"{cond}_mean_n_ge{int(D)}_since_4p51"] = float(v["mu_mean"][_dj(v, D)][-1]) for D in (190., 440., 740., 1000.): S = v["surv"][_dj(v, D)] for T in cfg.T_REPORT: row[f"{cond}_P_none_ge{int(D)}_after_{T}"] = float(np.interp(T, v["T"], S)) row[f"{cond}_median_last_ge{int(D)}"] = model.last_quantile(v, D, 0.5) row[f"{cond}_mc_se_P_ge{int(D)}_4.2"] = float(np.interp(4.2, r["T"], r["se_by"][cond][_dj(r, D)])) rows.append(row) df = pd.DataFrame(rows) _save(df.round(5), "model_grid.csv") return df def survival_tables(G): """Long-format survival curves, P(none after T) table and quantiles of the last impact, for every HSE treatment (single chronology x SFD combinations only for the two main ones).""" Tg = model.T_GRID Tout = np.round(np.arange(3.0, 4.5101, 0.01), 2) cur, pn, lq = [], [], [] for cond in CONDS: for name, r in G[cond].items(): if cond not in MAIN and is_combo(name): continue for D in D_TABLE: j = _dj(r, D) S = r["surv"][j] mu = r["mu_mean"][j] if cond in MAIN: for T, sv in zip(Tout, np.interp(Tout, Tg, S)): cur.append((cond, name, D, T, sv)) row = dict(hse=cond, group=name, D_km=D) for T in T_TABLE: row[f"P_none_after_{T}"] = float(np.interp(T, Tg, S)) pn.append(row) q = {f"q{int(p * 100)}": model.last_quantile(r, D, p) for p in (0.1, 0.25, 0.5, 0.75, 0.9)} lq.append(dict(hse=cond, group=name, D_km=D, **q, P_none_since_4p51=float(S[-1]), P_none_since_4p404=float(np.interp(T_W, Tg, S)), expected_since_4p51=float(mu[-1]), expected_since_4p404=float(np.interp(T_W, Tg, mu)), expected_since_4p2=float(np.interp(4.2, Tg, mu)), expected_since_3p9=float(np.interp(3.9, Tg, mu)), expected_since_3p48=float(np.interp(3.48, Tg, mu)), model_weight=float(r["weight"]))) cur = pd.DataFrame(cur, columns=["hse", "group", "D_km", "T_Ga", "P_none_younger"]) _save(cur.round(5), "survival_curves.csv") pn = pd.DataFrame(pn) _save(pn.round(4), "p_none_table.csv") lq = pd.DataFrame(lq) _save(lq.round(4), "last_impact_quantiles.csv") return pn, lq # --------------------------------------------------------------------------------------------- # Bayes factors # --------------------------------------------------------------------------------------------- def bf_tables(G, t_w=T_W, f=F, fname="bf_table.csv"): Tg = model.T_GRID rows = [] for cond in MAIN: for name, r in G[cond].items(): for D in D_BF: j = _dj(r, D) S = r["surv"][j] for lab, t_ev in cfg.LIFE_DATES.items(): ws = bayes.window_stats(Tg, S, t_w, t_ev) n_after = float(np.interp(t_ev, Tg, r["mu_mean"][j])) for m in M_LIST: L = bayes.likelihoods(Tg, S, t_w, t_ev, f, m) rows.append(dict(hse=cond, group=name, D_km=D, evidence=lab, t_ev=t_ev, m=m, B_survive=L["Lf_surv"] / L["Ls_surv"], B_sterilise=L["Lf_ster"] / L["Ls_ster"], B_mixed_prior50=bayes.bf_mixed(L, 0.5), B_sterilise_pointwise=bayes.bf_sterilising(Tg, S, t_w, t_ev, f, m, "pointwise"), P_no_impact_after_evidence=ws["p_ok"], expected_impacts_after_evidence=n_after, mean_window_Gyr=ws["mean_window"], median_window_Gyr=ws["median_window"], ster_vs_surv_if_fast=L["Lf_ster"] / L["Lf_surv"], ster_vs_surv_if_slow=L["Ls_ster"] / L["Ls_surv"])) df = pd.DataFrame(rows) _save(df.round(5), fname) return df def enabled_table(G): Tg = model.T_GRID rows = [] for cond in MAIN: for name in ("ALL", "sfd:mainbelt", "sfd:lunar_pl300", "sfd:lunar_pl600"): r = G[cond][name] for D in (440., 740., 1000.): S = r["surv"][_dj(r, D)] for W in (0.01, 0.03, 0.1, 0.2, 0.5): for lab, t_ev in cfg.LIFE_DATES.items(): for m in M_LIST: rows.append(dict(hse=cond, group=name, D_km=D, W_Gyr=W, evidence=lab, t_ev=t_ev, m=m, B_enabled=bayes.bf_enabled(Tg, S, T_W, t_ev, F, m, W), B_survive=bayes.bf_fixed(T_W, t_ev, F, m))) df = pd.DataFrame(rows) _save(df.round(4), "bf_impact_enabled.csv") return df def luca_date_marginal(G): """LUCA-age uncertainty integrated out (Moody et al. 2024: 4.20 Ga, 95% 4.09-4.33).""" from scipy import stats Tg = model.T_GRID n = 20000 # stratified quantiles: no Monte Carlo noise z = stats.norm.ppf((np.arange(n) + 0.5) / n) c, lo, hi = cfg.LUCA_PAPER luca = c + np.where(z < 0, (c - lo) / 1.959964, (hi - c) / 1.959964) * z tw, ff = T_W, F rows = [] for cond in MAIN: r = G[cond]["ALL"] for D in D_BF: S = r["surv"][_dj(r, D)] p_ok = float(np.mean(np.interp(luca, Tg, S))) for m in (0, 1, 2, 3, 4): b, used = bayes.bf_date_marginal(Tg, S, m, luca, tw, ff) b0, _ = bayes.bf_date_marginal(Tg, np.ones_like(Tg), m, luca, tw, ff) rows.append(dict(hse=cond, group="ALL", D_km=D, m=m, B_sterilise_datemarg=b, B_survive_datemarg=b0, P_no_impact_after_LUCA_datemarg=p_ok, draws_used=used)) df = pd.DataFrame(rows) _save(df.round(4), "bf_luca_date_marginal.csv") return df def ce6_pairs_table(bf): """Effect of the Chang'e-6 re-dating of SPA, like for like.""" rows = [] for cond in MAIN: for old, new, lab in CE6_PAIRS: for D in (190., 440., 740.): for ev in ("LUCA 4.2 Ga (Moody+ 2024)", "Dresser 3.48 Ga (accepted)"): q = bf[(bf.hse == cond) & (bf.D_km == D) & (bf.evidence == ev) & (bf.m == 1)] b0 = float(q[q.group == f"chron:{old}"].B_sterilise.iloc[0]) b1 = float(q[q.group == f"chron:{new}"].B_sterilise.iloc[0]) rows.append(dict(hse=cond, pair=lab, D_km=D, evidence=ev, B_before=b0, B_after=b1, change_percent=100.0 * (b1 / b0 - 1.0))) df = pd.DataFrame(rows) _save(df.round(3), "ce6_like_for_like.csv") return df # --------------------------------------------------------------------------------------------- # Terrestrial checks # --------------------------------------------------------------------------------------------- def spherule_table(G): """Expected number of Earth impactors above a size in the two intervals with spherule layers.""" Tg = model.T_GRID rows = [] obs = {(2.44, 2.68): [s for s in cfg.SPHERULES if 2.44 <= s[1] <= 2.68], (3.2, 3.5): [s for s in cfg.SPHERULES if 3.2 <= s[1] <= 3.5]} names = ["ALL"] + [f"chron:{c}" for c in chron.ALL] + [f"sfd:{s}" for s in model.SFD_FAMILIES] for cond in MAIN: for name in names: r = G[cond][name] for (t1, t2), layers in obs.items(): for d in (10., 20., 40.): mu = r["mu_mean"][_dj(r, d)] pred = float(np.interp(t2, Tg, mu) - np.interp(t1, Tg, mu)) n_min = sum(1 for s in layers if s[2] >= d) # layers whose smallest size estimate >= d n_max = sum(1 for s in layers if s[3] >= d) # layers whose largest size estimate >= d rows.append(dict(hse=cond, group=name, t_young=t1, t_old=t2, d_km=d, predicted=pred, layers_total=len(layers), layers_dmin_ge=n_min, layers_dmax_ge=n_max)) df = pd.DataFrame(rows) _save(df.round(3), "spherule_check.csv") return df def hse_table(grid): """Realised late-accreted mass by size distribution and by chronology: where it falls relative to the HSE range, and how many large impactors the realisations inside each treatment contain.""" rows = [] sets = [(s, [r for (c, ss), r in grid.items() if ss == s]) for s in model.SFD_FAMILIES] sets += [(f"[chron] {c}", [r for (cc, ss), r in grid.items() if cc == c]) for c in chron.ALL] sets += [("ALL", list(grid.values()))] for lab, rs in sets: m = np.concatenate([r["mass"] for r in rs]) row = dict(sfd=lab, mass_median_kg=float(np.median(m)), mass_p05_kg=float(np.quantile(m, 0.05)), mass_p95_kg=float(np.quantile(m, 0.95)), frac_below_upper_3e22=float((m <= cfg.HSE_MASS[1]).mean()), frac_within_0p7_3e22=float(((m >= cfg.HSE_MASS[0]) & (m <= cfg.HSE_MASS[1])).mean()), frac_below_lower_0p7e22=float((m < cfg.HSE_MASS[0]).mean()), frac_above_upper_3e22=float((m > cfg.HSE_MASS[1]).mean())) for cond in MAIN: ok = np.concatenate([r["ok_by"][cond] for r in rs]) nge = np.concatenate([r["n_ge"] for r in rs]) for D in (440., 500., 1000., 2000.): j = _dj(rs[0], D) row[f"{cond}_mean_n_ge{int(D)}_since_4p51"] = float(nge[ok, j].mean()) if ok.any() else np.nan j1 = _dj(rs[0], 1000.) row[f"{cond}_frac_with_any_ge1000"] = float((nge[ok, j1] > 0).mean()) if ok.any() else np.nan # lunar check: impactors >= 250 km (clearly larger than SPA's) expected on the Moon since SPA lam = np.concatenate([r["lam_moon"] for r in rs]) row[f"{cond}_moon_expected_ge250_since_SPA"] = float(lam[ok].mean()) if ok.any() else np.nan row[f"{cond}_moon_P_none_ge250_since_SPA"] = float(np.exp(-lam[ok]).mean()) if ok.any() else np.nan rows.append(row) df = pd.DataFrame(rows) _save(df, "hse_check.csv") return df def marchi_table(G): """The quantity Marchi et al. 2014 report (date of the last impactor > 500 km; number of impactors > 500 and > 1000 km) under every HSE treatment.""" rows = [] for cond in CONDS: for name in ["ALL"] + [f"sfd:{s}" for s in model.SFD_FAMILIES]: r = G[cond][name] rows.append(dict(hse=cond, group=name, median_last_ge500=model.last_quantile(r, 500., 0.5), age_10pct_younger_ge500=model.last_quantile(r, 500., 0.9), P_last_ge500_younger_than_4p0=1.0 - float(np.interp(4.0, r["T"], r["surv"][_dj(r, 500.)])), n_ge500_since_4p51=float(r["mu_mean"][_dj(r, 500.)][-1]), n_ge1000_since_4p51=float(r["mu_mean"][_dj(r, 1000.)][-1]), mixture_weight=float(r["weight"]))) df = pd.DataFrame(rows) _save(df.round(4), "marchi2014_comparison.csv") return df # --------------------------------------------------------------------------------------------- # Robustness # --------------------------------------------------------------------------------------------- def _metrics(r, t_w=T_W, f=F, convention="generative"): """Headline quantities of one survival-function set (a mixture or a single model).""" Tg = model.T_GRID S190, S440, S740 = (r["surv"][_dj(r, D)] for D in (190., 440., 740.)) q = lambda D, p=0.5: model.last_quantile(r, D, p) bs = lambda S, t, m=1: bayes.bf_sterilising(Tg, S, t_w, t, f, m, convention) mix = lambda S, t: bayes.bf_mixed(bayes.likelihoods(Tg, S, t_w, t, f, 1), 0.5) return dict(P_none_ge440_after_4p404=float(np.interp(4.404, Tg, S440)), P_none_ge440_after_4p2=float(np.interp(4.2, Tg, S440)), P_none_ge440_after_3p9=float(np.interp(3.9, Tg, S440)), P_none_ge440_after_3p8=float(np.interp(3.8, Tg, S440)), P_none_ge190_after_4p2=float(np.interp(4.2, Tg, S190)), P_none_ge190_after_3p8=float(np.interp(3.8, Tg, S190)), median_last_ge440=q(440.), median_last_ge190=q(190.), age_10pct_younger_ge440=q(440., 0.9), age_10pct_younger_ge190=q(190., 0.9), B_LUCA_m1_ster190=bs(S190, 4.2), B_LUCA_m1_ster440=bs(S440, 4.2), B_LUCA_m1_ster740=bs(S740, 4.2), B_3p48_m1_ster440=bs(S440, 3.48), B_3p48_m1_ster190=bs(S190, 3.48), B_LUCA_m1_mixed50_440=mix(S440, 4.2), B_3p48_m1_mixed50_190=mix(S190, 3.48), B_LUCA_m1_survive=bayes.bf_fixed(t_w, 4.2, f, 1), B_3p48_m1_survive=bayes.bf_fixed(t_w, 3.48, f, 1)) def robustness(B, grid, G, say=lambda *a: None): """Every variant is evaluated under both main HSE treatments. kind: baseline the main grid hse another treatment of the HSE mass constraint (main grid, no rerun) bayes a change on the Bayes-factor side only (main grid, no rerun) mixture variant a bombardment-side change, all 21 models rerun with N_ROB draws reference the baseline rerun with N_ROB draws and another seed subset part of the model set (not comparable with the full mixture)""" rows = [] add = lambda kind, variant, cond, r, **kw: rows.append(dict(kind=kind, variant=variant, hse=cond, **_metrics(r, **kw))) for cond in MAIN: add("baseline", "main grid, all models", cond, G[cond]["ALL"]) for cond in CONDS: if cond not in MAIN: add("hse", cfg.HSE_LABEL[cond], cond, G[cond]["ALL"]) for cond in MAIN: r = G[cond]["ALL"] add("bayes", "Bayes factor: pointwise convention", cond, r, convention="pointwise") add("bayes", "habitable from 4.50 Ga instead of 4.404", cond, r, t_w=4.50) add("bayes", "habitable from 4.30 Ga instead of 4.404", cond, r, t_w=4.30) add("bayes", "1.73 Gyr of habitability left instead of 0.9", cond, r, f=1.73) add("subset", "all models without the main-belt shape (main grid)", cond, G[cond]["ALL without main-belt shape"]) add("subset", "main-belt shape only (main grid)", cond, G[cond]["sfd:mainbelt"]) variants = [ ("basins follow N1(T) smoothly (no individual basin ages)", dict(smooth=True)), ("densities from Fassett+ 2012 N(20) instead of Orgel+ ages", dict(density_source="fassett")), ("undated basins drawn from all dated basins", dict(undated="all")), ("undated basins dropped", dict(undated="drop")), ("Earth/Moon ratio fixed at 17", dict(r_fixed=17.0)), ("Earth/Moon ratio fixed at 23", dict(r_fixed=23.0)), ("impactor sizes x 0.8", dict(scale_factor=0.8)), ("impactor sizes x 1.25", dict(scale_factor=1.25)), ("no flux before the oldest basin (SPA)", dict(pre_record=False)), ("no basin-forming flux after the youngest basin", dict(post_record=False)), ("basin-forming flux after the youngest basin x 3 (spherule-layer level)", dict(post_scale=3.0)), # (small chunks: this variant samples very many early impactors per draw) ("anchored: pre-SPA e-fold 0.03-0.07 Gyr (steep)", dict(tau_pre_range=(0.03, 0.07), chunk=250)), ("anchored: pre-SPA e-fold 0.3-1.0 Gyr (flat)", dict(tau_pre_range=(0.3, 1.0))), ("power laws end at 940-1000 km (as the main belt does)", dict(dmax_range=(940.0, 1000.0))), ("power laws end at 1000-6000 km", dict(dmax_range=(1000.0, 6000.0))), ("main-belt tail not bootstrapped (observed counts taken as exact)", dict(mb_boot=False)), ("sawtooth curve as published with Orgel densities (Nectaris 4.32 Ga)", dict(chron_swap={"sawtooth2012": "sawtooth_aspub"})), ("sawtooth curve x 1.61 at all ages (Nectaris 4.1 Ga, Imbrium 3.79 Ga)", dict(chron_swap={"sawtooth2012": "sawtooth_uniform"})), ] for k, (lab, kw) in enumerate(variants): g = run_grid(B, N_ROB, cfg.SEED + 1000 + k, **kw) for cond in MAIN: add("mixture variant", lab, cond, model.mixture(list(g.values()), cond)) say(f" variant {k + 1}/{len(variants)}: {lab}") g = run_grid(B, N_ROB, cfg.SEED + 999) for cond in MAIN: add("reference", f"baseline rerun with {N_ROB} draws and another seed", cond, model.mixture(list(g.values()), cond)) g = run_grid(B, N_ROB, cfg.SEED + 998, d_ref_basin=600.0, only_sfd=["mainbelt"]) for cond in MAIN: add("subset", "main-belt shape only, normalised at 600 km basins", cond, model.mixture(list(g.values()), cond)) # the sawtooth row on its own under the three normalisations for lab, cname in (("sawtooth chronology alone: Nectaris pinned at 4.1 Ga (main grid)", None), ("sawtooth chronology alone: as published with Orgel densities", "sawtooth_aspub"), ("sawtooth chronology alone: curve x 1.61 at all ages", "sawtooth_uniform")): if cname is None: for cond in MAIN: add("subset", lab, cond, G[cond]["chron:sawtooth2012"]) continue rs = [model.run_chunked(cname, s, 3 * N_ROB, np.random.default_rng([cfg.SEED + 997, j]), basins=B, d_grid=D_FULL) for j, s in enumerate(model.SFD_FAMILIES)] for cond in MAIN: add("subset", lab, cond, model.mixture(rs, cond)) df = pd.DataFrame(rows) _save(df.round(4), "robustness.csv") return df # --------------------------------------------------------------------------------------------- def main(verbose=True): t0 = time.time() say = (lambda *a: print(f"[{time.time() - t0:6.1f}s]", *a, flush=True)) if verbose else (lambda *a: None) res = {} B = model.load_basins() res["basins"] = B res["basin_table"] = basin_tables(B) res["ratio_table"], res["ratio_stats"] = ratio_table() res["scaling"] = scaling_table() res["sfd_fits"], res["sfd_pred"], res["mainbelt_counts"] = sfd_tables(B) say("static tables written") grid = run_grid(B, N_MAIN, cfg.SEED) say(f"main grid done: {len(grid)} combinations x {N_MAIN} draws") G = {cond: make_groups(grid, cond) for cond in CONDS} res["grid"], res["G"] = grid, G res["model_grid"] = model_grid_table(grid) res["p_none"], res["last_q"] = survival_tables(G) res["bf"] = bf_tables(G) res["bf_enabled"] = enabled_table(G) res["bf_luca_dm"] = luca_date_marginal(G) res["ce6"] = ce6_pairs_table(res["bf"]) res["spherules"] = spherule_table(G) res["hse"] = hse_table(grid) res["marchi"] = marchi_table(G) say("survival, Bayes-factor and check tables written") res["robustness"] = robustness(B, grid, G, say) say("robustness variants done") # ---- headline numbers --------------------------------------------------------------------- pn, lq, bf, rob, en = res["p_none"], res["last_q"], res["bf"], res["robustness"], res["bf_enabled"] g = lambda df, **kw: df.loc[np.logical_and.reduce([df[k] == v for k, v in kw.items()])] chron_groups = [f"chron:{c}" for c in chron.ALL] sfd_groups = [f"sfd:{x}" for x in model.SFD_FAMILIES] combos = [f"{c}|{s}" for c in chron.ALL for s in model.SFD_FAMILIES] def rng_over(df, col, groupnames, **kw): v = g(df[df.group.isin(groupnames)], **kw)[col] return [float(v.min()), float(v.max())] post_n = float(np.mean([np.nanmean(r["post_n"]) for (c, s), r in grid.items() if s != "lunar_pl600"])) sp = res["sfd_pred"] summ = dict(n_draw_main=N_MAIN, n_draw_robustness=N_ROB, seed=cfg.SEED, t_w=T_W, future=F, n_basins_neumann=int(len(B)), n_basins_ge300=int((B.D_main_km >= 300).sum()), n_basins_ge600=int((B.D_main_km >= 600).sum()), n_dated=int(B.dated.sum()), n_dated_ge300=int((B.dated & (B.D_main_km >= 300)).sum()), n_dated_ge600=int((B.dated & (B.D_main_km >= 600)).sum()), earth_moon_ratio=res["ratio_stats"], scaling=res["scaling"], max_abs_diff_sampled_vs_analytic=float(res["model_grid"].max_abs_diff_sampled_vs_analytic.max()), mainbelt_counts=res["mainbelt_counts"], sawtooth_scale_at_nectaris=float(chron.MODELS["sawtooth2012"].s_nectaris), lunar_basins_ge300_expected_after_youngest=post_n, poisson_p_zero_basins_after_youngest=float(np.exp(-post_n)), # mass of the impactors of the catalogued lunar basins (nominal sizes), and the mass the # same population puts on Earth: far below the HSE veneer, which must therefore sit # in bodies larger than the SPA impactor lunar_basin_impactor_mass_kg=float((np.pi / 6.0 * cfg.RHO_IMP * 1e9 * B.d_imp_km ** 3).sum()), spa_impactor_mass_kg=float(np.pi / 6.0 * cfg.RHO_IMP * 1e9 * B.d_imp_km.max() ** 3), earth_mass_from_basin_sized_impactors_kg=float( res["ratio_stats"]["median"] * (np.pi / 6.0 * cfg.RHO_IMP * 1e9 * B.d_imp_km ** 3).sum()), moon_ge250km_since_SPA={cond: dict( expected=float(res["hse"].set_index("sfd").loc["ALL", f"{cond}_moon_expected_ge250_since_SPA"]), P_none=float(res["hse"].set_index("sfd").loc["ALL", f"{cond}_moon_P_none_ge250_since_SPA"])) for cond in MAIN}, megabasins_ge1000km=dict(observed=int(g(sp, D_basin_km=1000.0).observed.iloc[0]), mainbelt_predicted=float(g(sp, D_basin_km=1000.0).mainbelt.iloc[0]), poisson_p=float(g(sp, D_basin_km=1000.0).poisson_p_le_obs_mainbelt.iloc[0])), hse_pass_fraction_by_family={cond: {s: float(np.mean([r["frac_by"][cond] for (c, ss), r in grid.items() if ss == s])) for s in model.SFD_FAMILIES} for cond in CONDS}, mixture_weights_ALL={cond: dict(zip([f"{c}|{s}" for (c, s) in grid], [float(x) for x in G[cond]["ALL"]["weights"]])) for cond in MAIN}) by = {} for cond in CONDS: d = {} for D in D_TABLE: row = g(lq, hse=cond, group="ALL", D_km=D).iloc[0] d[f"ge{int(D)}"] = dict( median_last_ALL=float(row.q50), q25_ALL=float(row.q25), q75_ALL=float(row.q75), q90_ALL=float(row.q90), median_last_range_chron=rng_over(lq, "q50", chron_groups, hse=cond, D_km=D), median_last_range_sfd=rng_over(lq, "q50", sfd_groups, hse=cond, D_km=D), q90_range_sfd=rng_over(lq, "q90", sfd_groups, hse=cond, D_km=D), P_none_since_4p404_ALL=float(row.P_none_since_4p404), expected_since_4p51_ALL=float(row.expected_since_4p51), expected_since_4p404_ALL=float(row.expected_since_4p404), expected_since_4p2_ALL=float(row.expected_since_4p2), P_none_ALL={str(T): float(g(pn, hse=cond, group="ALL", D_km=D)[f"P_none_after_{T}"].iloc[0]) for T in T_TABLE}, P_none_range_chron={str(T): rng_over(pn, f"P_none_after_{T}", chron_groups, hse=cond, D_km=D) for T in T_TABLE}, P_none_range_sfd={str(T): rng_over(pn, f"P_none_after_{T}", sfd_groups, hse=cond, D_km=D) for T in T_TABLE}) by[cond] = d summ["by_hse"] = by bfs = {} for cond in MAIN: bfs[cond] = {} for lab in cfg.LIFE_DATES: bfs[cond][lab] = {} for D in D_BF: sub = g(bf, hse=cond, group="ALL", D_km=D, evidence=lab) bfs[cond][lab][f"ge{int(D)}"] = dict( B_survive=[float(x) for x in sub.B_survive], B_sterilise=[float(x) for x in sub.B_sterilise], B_mixed50=[float(x) for x in sub.B_mixed_prior50], P_no_impact_after=float(sub.P_no_impact_after_evidence.iloc[0]), expected_after=float(sub.expected_impacts_after_evidence.iloc[0]), mean_window=float(sub.mean_window_Gyr.iloc[0]), odds_survive_vs_sterilise_if_fast_m1=float(1.0 / sub[sub.m == 1].ster_vs_surv_if_fast.iloc[0]), odds_survive_vs_sterilise_if_slow_m1=float(1.0 / sub[sub.m == 1].ster_vs_surv_if_slow.iloc[0]), B_sterilise_m1_range_chron=rng_over(bf, "B_sterilise", chron_groups, hse=cond, D_km=D, evidence=lab, m=1), B_sterilise_m1_range_sfd=rng_over(bf, "B_sterilise", sfd_groups, hse=cond, D_km=D, evidence=lab, m=1), B_sterilise_m1_range_combos=rng_over(bf, "B_sterilise", combos, hse=cond, D_km=D, evidence=lab, m=1)) summ["bayes"] = bfs # impact-enabled origin: ranges for windows up to 200 Myr and for the 500 Myr window summ["impact_enabled"] = {} for lab in cfg.LIFE_DATES: q = en[(en.evidence == lab)] a, b = q[q.W_Gyr <= 0.2].B_enabled, q[q.W_Gyr == 0.5].B_enabled summ["impact_enabled"][lab] = dict(range_W_le_200Myr=[float(a.min()), float(a.max())], range_W_500Myr=[float(b.min()), float(b.max())]) # ranges over the bombardment-side mixture variants, per treatment mv = rob[rob.kind.isin(["mixture variant", "baseline"])] summ["robustness_ranges"] = {cond: {c: [float(mv[mv.hse == cond][c].min()), float(mv[mv.hse == cond][c].max())] for c in mv.columns if c not in ("kind", "variant", "hse")} for cond in MAIN} summ["n_mixture_variants"] = int((rob.kind == "mixture variant").sum() // len(MAIN)) summ["m_list"] = M_LIST with open(os.path.join(cfg.DATA, "results.json"), "w") as fh: json.dump(summ, fh, indent=1) res["summary"] = summ say("results.json written") return res