"""Bombardment model: lunar basin record -> large impacts on Earth. Expected number of impactors with diameter >= D striking Earth more recently than age T: mu_D(T) = R * S(D)/S(d_ref) * C(T) R Earth/Moon ratio of impact numbers for the same impactor population (gravitational focusing, derived in earth_moon_ratio) S(D) cumulative impactor size-frequency distribution (three families, see SFD_FAMILIES) d_ref impactor diameter that makes a lunar basin of main-ring diameter D_ref (300 or 600 km) C(T) number of lunar "reference basins" (main ring >= D_ref) younger than T: within the record: sum over catalogued basins of w_j 1[t_j < T], w_j ~ Exp(1) (Bayesian bootstrap = Poisson uncertainty of the lunar sample) before the record: K [N1(T) - N1(T_SPA)], K = sum(w) / N1_SPA after the last basin: K N1(min(T, T_youngest)) (small; see report) with N1(T) the chronology model. In the "smooth" variant the first term is replaced by K N1(T) throughout, i.e. basins are assumed to form in proportion to small craters. Estimates of P(no impact >= D younger than T) returned on an age grid: surv E[exp(-mu_D(T))] over all sampled parameters (analytic Poisson, no mass constraint) surv_s the same from sampled realisations (cross-check of the sampling code) surv_by[c] the same probability conditional on the realised late-accreted mass lying inside the bounds of treatment c of the highly-siderophile-element (HSE) constraint (config.HSE_CONDS: upper bound only, both bounds, none, ...). For every parameter draw the sizes of all Earth impactors >= D_REAL (150 km) since 4.51 Ga are sampled; the realised mass is their sum plus the expected mass of the smaller bodies; given the sizes, the ages are independent with distribution G(T) = C(T)/C(4.51), so P(no impact >= D younger than T | sizes) = (1 - G(T))^n_D is averaged over the realisations that satisfy the mass bounds. All treatments use the same realisations. The main-belt shape rests on very few bodies at the sizes that matter (3 bodies >= 440 km). Its sampling error is propagated with a Bayesian bootstrap: in every draw each main-belt body >= D_BOOT (100 km) gets an Exp(1) weight, so the number of bodies above a size is Gamma-distributed around the observed count. The finite radius of the impactor adds (1 + D/2R_E)^2 to Earth's cross-section (7% at 440 km); this is neglected. """ from __future__ import annotations import os import unicodedata import numpy as np import pandas as pd import chronology as chron import config as cfg # age grid (Ga): coarse where nothing happens, 5 Myr steps over the bombardment era T_GRID = np.round(np.concatenate([np.arange(0.0, 3.0, 0.05), np.arange(3.0, cfg.T_MOON + 1e-9, 0.005)]), 4) D_REF_BASIN = 300.0 # km (main-ring diameter) D_REAL = 150.0 # km: Earth impactors at least this large are sampled individually D_BOOT = 100.0 # km: main-belt bodies at least this large get Bayesian-bootstrap weights D_MOON_CHECK = 250.0 # km: impactor size used for the "nothing larger than SPA on the Moon" check # Published impactor-size estimates for lunar basins / craters used to calibrate d = a D^b. # (D km, d km, source) SCALING_POINTS = [ (160.0, 9.7, "Johnson et al. 2016, eq. 4 at 22 km/s, 45 deg"), (300.0, 19.7, "Johnson et al. 2016, eq. 4 at 22 km/s, 45 deg"), (326.0, 20.0, "Schrodinger: Miljkovic et al. 2016 / Johnson et al. 2016b via Nesvorny et al. 2023"), (937.0, 57.0, "Orientale: 50-64 km, Miljkovic et al. 2016 / Johnson et al. 2016b via Nesvorny et al. 2023"), (1321.0, 100.0, "Imbrium: >~100 km, Miljkovic et al. 2013, 2016 via Nesvorny et al. 2023"), (2400.0, 170.0, "SPA: Potter et al. 2012 via Marchi et al. 2014"), ] SFD_FAMILIES = ["mainbelt", "lunar_pl300", "lunar_pl600"] SFD_LABEL = {"mainbelt": "main-belt shape (JPL SBDB), nothing larger than Ceres", "lunar_pl300": "power law fitted to all 43 basins >= 300 km (shallow), continued", "lunar_pl600": "power law fitted to the 17 basins >= 600 km (steep), continued"} # main-ring diameter of the lunar reference basins for each family SFD_DREF = {"mainbelt": 300.0, "lunar_pl300": 300.0, "lunar_pl600": 600.0} # --------------------------------------------------------------------------------------------- def _key(s): s = unicodedata.normalize("NFKD", str(s)).encode("ascii", "ignore").decode() return s.lower().replace(" ", "").replace("-", "").replace("'", "") def load_basins(): """Neumann et al. 2015 basins merged with Fassett et al. 2012 N(20) and Orgel et al. 2018 ages.""" n = pd.read_csv(os.path.join(cfg.DATA, "neumann2015_basins.csv")) f = pd.read_csv(os.path.join(cfg.DATA, "fassett2012_table1.csv")) o = pd.read_csv(os.path.join(cfg.DATA, "orgel2018_lpsc_ages.csv")) n["key"], f["key"], o["key"] = n.name.map(_key), f.name.map(_key), o.name.map(_key) assert not (set(f.key) - set(n.key)) and not (set(o.key) - set(n.key)) m = (n.merge(f[["key", "period", "N20", "N20_err", "may_be_low", "bracketed"]], on="key", how="left") .merge(o[["key", "age_neukum_Ga", "err_plus", "err_minus"]], on="key", how="left")) m = m.sort_values("D_main_km", ascending=False).reset_index(drop=True) a, b = scaling_fit() m["d_imp_km"] = a * m.D_main_km ** b m["dated"] = m.age_neukum_Ga.notna() nk = chron.MODELS["neukum1983"] m["N1_orgel"] = nk.n1(m.age_neukum_Ga) return m.drop(columns="key") def scaling_fit(): """Least-squares power law d = a D^b through the published (basin, impactor) pairs.""" D = np.array([p[0] for p in SCALING_POINTS]) d = np.array([p[1] for p in SCALING_POINTS]) b, lna = np.polyfit(np.log(D), np.log(d), 1) return float(np.exp(lna)), float(b) def earth_moon_ratio(v_moon, a_re=60.3): """Ratio of impact numbers Earth : Moon for a population arriving with lunar impact speed v_moon (km/s). Two-body focusing: N ~ pi R^2 (1 + v_esc^2 / v_inf^2); at the Moon the approach speed already includes Earth's well, v_esc,E^2 R_E / a.""" v_moon = np.asarray(v_moon, float) vE2_at_moon = cfg.VESC_EARTH ** 2 / np.asarray(a_re, float) vinf2 = v_moon ** 2 - cfg.VESC_MOON ** 2 - vE2_at_moon f_earth = 1.0 + cfg.VESC_EARTH ** 2 / vinf2 f_moon = 1.0 + (cfg.VESC_MOON ** 2 + vE2_at_moon) / vinf2 return (cfg.R_EARTH / cfg.R_MOON) ** 2 * f_earth / f_moon def mainbelt_sizes(): s = pd.read_csv(os.path.join(cfg.DATA, "sbdb_asteroids_d10.csv")) d = s.loc[s["class"].isin(["MBA", "IMB", "OMB"]), "diameter"].to_numpy(float) return np.sort(d)[::-1] class MainBeltSFD: """Empirical cumulative N(>= d) of the main belt, log-log interpolated between objects; zero above the largest body (Ceres, 939 km). The plain methods (N, mass_per_ref) use the observed bodies with unit weight. The *_w methods take per-draw weights for the n_L bodies >= D_BOOT (Bayesian bootstrap of the sparse large end); with unit weights they reduce to the plain ones.""" def __init__(self): self.d = mainbelt_sizes() # descending self.n = len(self.d) self.lnd_desc = np.log(self.d) self.lnd = np.log(self.d[::-1]) # ascending self.lnN = np.log(np.arange(len(self.d), 0, -1.0)) self.dmax = self.d[0] self.mass_sorted = np.pi / 6.0 * cfg.RHO_IMP * (self.d * 1e3) ** 3 # kg, descending d self.cm = np.concatenate([[0.0], np.cumsum(self.mass_sorted)]) # mass of the j largest self.nL = int((self.d >= D_BOOT).sum()) self.n_real = int((self.d >= D_REAL).sum()) assert self.n_real + 1 <= self.nL def N(self, d): d = np.asarray(d, float) out = np.exp(np.interp(np.log(d), self.lnd, self.lnN)) return np.where(d > self.dmax, 0.0, out) # ---- weighted (bootstrap) versions ------------------------------------------------------- def weights(self, rng, n_draw, boot=True): """(n_draw, n_L) weights of the bodies >= D_BOOT and their cumulative sums.""" w = rng.exponential(1.0, (n_draw, self.nL)) if boot else np.ones((n_draw, self.nL)) return w, np.cumsum(w, axis=1) def _cum(self, Wc, j): """Weighted number of bodies at least as large as body j (0-based rank), per draw.""" base = np.take_along_axis(Wc, np.minimum(j, self.nL - 1), axis=1) return np.where(j < self.nL, base, base + (j - (self.nL - 1))) def N_w(self, d, Wc): """Weighted N(>= d) per draw; d has shape (n_draw, m) or (m,).""" d = np.asarray(d, float) if d.ndim == 1: d = np.broadcast_to(d[None, :], (Wc.shape[0], d.size)) i = np.searchsorted(-self.d, -d, side="right") # number of bodies >= d lo = np.clip(i - 1, 0, self.n - 1) # smallest body >= d hi = np.clip(i, 0, self.n - 1) # largest body < d c_lo, c_hi = self._cum(Wc, lo), self._cum(Wc, hi) dl = self.lnd_desc[lo] - self.lnd_desc[hi] frac = np.clip(np.where(dl > 0, (self.lnd_desc[lo] - np.log(d)) / np.where(dl > 0, dl, 1.0), 0.0), 0.0, 1.0) out = np.exp(np.log(c_lo) + frac * (np.log(c_hi) - np.log(c_lo))) return np.where(i == 0, 0.0, out) def mass_w(self, w, cnt): """Weighted mass (kg) of the cnt largest bodies, per draw; cnt has shape (n_draw,).""" Wm = np.concatenate([np.zeros((w.shape[0], 1)), np.cumsum(w * self.mass_sorted[None, :self.nL], axis=1)], axis=1) head = np.take_along_axis(Wm, np.minimum(cnt, self.nL)[:, None], axis=1)[:, 0] return head + (self.cm[np.maximum(cnt, self.nL)] - self.cm[self.nL]) def count_w(self, Wc, cnt): """Weighted number of the cnt largest bodies, per draw (cnt >= 1).""" return self._cum(Wc, (cnt - 1)[:, None])[:, 0] def sample_w(self, u, idx, Wc, n_lo): """Diameters >= D_REAL for impactors belonging to draws idx (u uniform on 0-1), from each draw's weighted shape; n_lo = N_w(D_REAL) per draw. The largest body (Ceres) is an atom.""" nr = self.n_real tab = Wc[:, :nr + 1] # cumulative weights down to the first body < D_REAL c = u * n_lo[idx] big = 2.0 * float(tab.max()) + 1.0 flat = (tab + big * np.arange(tab.shape[0])[:, None]).ravel() pos = np.searchsorted(flat, c + big * idx, side="left") - idx * (nr + 1) # entries of the row < c pos = np.clip(pos, 0, nr) j0 = np.maximum(pos - 1, 0) c0, c1 = tab[idx, j0], tab[idx, pos] with np.errstate(divide="ignore", invalid="ignore"): frac = np.where(pos > 0, (np.log(c) - np.log(c0)) / (np.log(c1) - np.log(c0)), 0.0) frac = np.clip(np.nan_to_num(frac), 0.0, 1.0) lnd = self.lnd_desc[j0] + frac * (self.lnd_desc[pos] - self.lnd_desc[j0]) return np.where(pos == 0, self.dmax, np.exp(lnd)) def pl_ratio(D, d_ref, q, d_max=np.inf): """S(D)/S(d_ref) for a cumulative power law of slope q truncated at d_max (no bodies above).""" D = np.asarray(D, float) top = (d_max / d_ref) ** (-q) return np.clip((D / d_ref) ** (-q) - top, 0.0, None) / (1.0 - top) def pl_mass_per_ref(d_ref, q, d_max, d_hi=None): """Mean mass (kg) per impactor >= d_ref carried by bodies d_ref <= d < d_hi (default d_max) of the truncated power law: dN = q d_ref^q d^-(q+1) dd / (1 - (d_max/d_ref)^-q).""" k = np.pi / 6.0 * cfg.RHO_IMP * 1e9 # kg per km^3 hi = d_max if d_hi is None else np.minimum(d_hi, d_max) e = 3.0 - q integ = np.where(np.abs(e) < 1e-9, np.log(hi / d_ref), (hi ** e - d_ref ** e) / np.where(np.abs(e) < 1e-9, 1.0, e)) return k * q * d_ref ** q * integ / (1.0 - (d_max / d_ref) ** (-q)) def pl_sample(u, d_lo, q, d_max): """Diameters between d_lo and d_max drawn from the truncated power law (u uniform on 0-1).""" a, b = d_lo ** (-q), d_max ** (-q) return (a - u * (a - b)) ** (-1.0 / q) # --------------------------------------------------------------------------------------------- def run_model(chron_name, sfd_name, n_draw, rng, basins=None, d_ref_basin=None, smooth=False, density_source="orgel", undated="pn", r_fixed=None, scale_sigma=0.2, scale_factor=1.0, pre_record=True, post_record=True, tau_pre_range=(0.07, 0.30), dmax_range=(1000.0, 3000.0), conds=None, post_scale=1.0, mb_boot=True, d_grid=cfg.D_GRID, return_draws=False): """Monte Carlo for one chronology x SFD family. Returns a dict with surv[D_index, T_index] = P(no impactor >= D younger than T), analytic, no mass constraint surv_s = same from the sampled realisations (cross-check of surv) surv_by[c] = same, conditional on the realised mass lying inside the bounds of HSE treatment c (conds: name -> (lower, upper) in kg; default config.HSE_CONDS) mu_mean, mu_by[c] = expected number of such impacts (unconditional, conditional) frac_by[c], ok_by[c] = fraction / mask of realisations satisfying treatment c mass (realised), mass_exp (expected), n_ref_earth, post_n, ... per-draw diagnostics Options smooth basins form in proportion to N1(T) instead of using their individual ages density_source 'orgel' (BNSC model ages -> N1) or 'fassett' (N20 x 1000) undated 'pn' undated basins get densities drawn from dated basins at least as dense as Nectaris; 'all' from all dated basins; 'drop' leaves them out pre_record False sets the flux to zero before the oldest catalogued basin post_record False sets the flux of basin-forming bodies to zero after the youngest basin post_scale multiplies the extrapolated flux after the youngest basin scale_factor multiplies every impactor diameter inferred from a basin diameter tau_pre_range e-folding time (Gyr) of the pre-SPA flux in the anchored chronologies dmax_range range (km) of the largest body of the power-law families (uniform) mb_boot False switches off the bootstrap of the main-belt bodies >= D_BOOT Diameters in d_grid smaller than d_ref (used only for the spherule-layer check) always use the main-belt shape normalised at d_ref. """ conds = cfg.HSE_CONDS if conds is None else conds B = load_basins() if basins is None else basins if d_ref_basin is None: d_ref_basin = SFD_DREF[sfd_name] B = B[B.D_main_km >= d_ref_basin].reset_index(drop=True) if undated == "drop": B = B[B.dated].reset_index(drop=True) nb = len(B) nk = chron.MODELS["neukum1983"] Tg = T_GRID nT = len(Tg) # ---- crater densities of the basins (Neukum-equivalent N1), one row per draw ------------- dated = B.dated.to_numpy() if density_source == "orgel": z = rng.standard_normal((n_draw, nb)) age0 = B.age_neukum_Ga.to_numpy() sig = np.where(z < 0, B.err_minus.to_numpy(), B.err_plus.to_numpy()) n1 = nk.n1(age0 + z * sig) # NaN for undated basins else: z = rng.standard_normal((n_draw, nb)) n1 = np.clip(B.N20.to_numpy() + z * B.N20_err.to_numpy(), 5.0, None) * 1e-6 * cfg.N1_OVER_N20 i_nec = int(np.where(B.name == "Nectaris")[0][0]) i_imb = int(np.where(B.name == "Imbrium")[0][0]) i_spa = int(np.where(B.name == "South Pole-Aitken")[0][0]) if density_source == "fassett": # Fassett et al. flag their SPA value as "may be low" (later basins erased its craters); # SPA is stratigraphically the oldest basin, so give it the largest density of the draw n1[:, i_spa] = np.where(dated, n1, -np.inf).max(axis=1) # no basin can be older than SPA: cap every density at SPA's n1 = np.where(dated, np.minimum(n1, n1[:, [i_spa]]), n1) # undated basins: borrow the density of a randomly chosen dated basin und = np.where(~dated)[0] if len(und): if undated == "pn": pool = np.where(dated & (B.N1_orgel.to_numpy() >= B.N1_orgel.iloc[i_nec] - 1e-12))[0] else: pool = np.where(dated)[0] pick = pool[rng.integers(0, len(pool), size=(n_draw, len(und)))] n1[:, und] = np.take_along_axis(n1, pick, axis=1) n1_spa = n1[:, i_spa] # ---- ages under the chosen chronology and the N1(T) curve of each draw ------------------- anchored = chron_name in chron.ANCHORED if not anchored: cf = chron.MODELS[chron_name] ages = cf.age(n1) n1_grid = np.broadcast_to(cf.n1(Tg), (n_draw, nT)) t_spa = ages[:, i_spa] else: t_spa0, t_spa_sig = cfg.SPA_CE6 if chron_name == "anchor_ce6" else cfg.SPA_JOY t_spa = rng.normal(t_spa0, t_spa_sig, n_draw) t_imb = rng.normal(*cfg.IMBRIUM_AGE, n_draw) t_nec = rng.uniform(3.95, np.minimum(4.20, t_spa - 0.03)) tau_pre = np.exp(rng.uniform(np.log(tau_pre_range[0]), np.log(tau_pre_range[1]), n_draw)) ages = np.empty_like(n1) n1_grid = np.empty((n_draw, nT)) for k in range(n_draw): a = (n1[k, i_imb], n1[k, i_nec], n1[k, i_spa], t_imb[k], t_nec[k], t_spa[k]) if not (a[0] < a[1] < a[2]): # rare: densities out of order in this draw s = np.sort([a[0], a[1], a[2]]) * np.array([1.0, 1.0001, 1.0002]) a = (s[0], s[1], s[2]) + a[3:] ages[k] = chron.anchored_ages(n1[k], *a) n1_grid[k] = chron.anchored_n1_grid(Tg, *a, tau_pre[k]) # ---- lunar basin-equivalents younger than T: C(T) ---------------------------------------- w = rng.exponential(1.0, (n_draw, nb)) wtot = w.sum(axis=1) K = wtot / n1_spa if smooth: C = K[:, None] * n1_grid if not pre_record: C = np.minimum(C, wtot[:, None]) if not post_record: t_y = ages.min(axis=1) n1_y = np.take_along_axis(n1_grid, np.clip(np.searchsorted(Tg, t_y), 0, nT - 1)[:, None], 1) C = np.clip(C - K[:, None] * n1_y, 0.0, None) post_n = np.full(n_draw, np.nan) else: idx = np.clip(np.searchsorted(Tg, ages, side="left"), 0, nT - 1) # first grid age >= basin age C = np.zeros((n_draw, nT)) np.add.at(C, (np.repeat(np.arange(n_draw), nb), idx.ravel()), w.ravel()) C = np.cumsum(C, axis=1) # basins with age <= T t_young = ages.min(axis=1) n1_young = np.take_along_axis(n1_grid, np.clip(np.searchsorted(Tg, t_young), 0, nT - 1)[:, None], 1) if post_record: C += post_scale * K[:, None] * np.minimum(n1_grid, n1_young) # after the last basin # lunar reference basins the extrapolation puts after the youngest catalogued basin post_n = post_scale * K * n1_young[:, 0] if post_record else np.zeros(n_draw) if pre_record: n1_at_spa = np.take_along_axis(n1_grid, np.clip(np.searchsorted(Tg, t_spa), 0, nT - 1)[:, None], 1) C += K[:, None] * np.clip(n1_grid - n1_at_spa, 0.0, None) * (Tg[None, :] > t_spa[:, None]) # ---- Earth/Moon ratio, impactor sizes, size distribution --------------------------------- if r_fixed is None: v_moon = rng.uniform(11.0, 22.0, n_draw) a_re = rng.uniform(*cfg.A_MOON_RE, n_draw) R = earth_moon_ratio(v_moon, a_re) else: R = np.full(n_draw, float(r_fixed)) a_sc, b_sc = scaling_fit() fsc = np.exp(rng.normal(0.0, scale_sigma, n_draw)) # shared scaling systematic d_ref = scale_factor * a_sc * d_ref_basin ** b_sc * fsc d_bas = a_sc * B.D_main_km.to_numpy() ** b_sc # (nb,) nominal impactor sizes # posterior of the lunar power-law slope from the basin sizes (Pareto likelihood, flat prior) slog = np.log(d_bas / (a_sc * d_ref_basin ** b_sc)).sum() q_lunar = rng.gamma(nb + 1.0, 1.0 / slog, n_draw) Dg = np.asarray(d_grid, float) mb = MainBeltSFD() w_mb, Wc_mb = mb.weights(rng, n_draw, boot=(mb_boot and sfd_name == "mainbelt")) n_at_ref = mb.N_w(d_ref[:, None], Wc_mb)[:, 0] ratio_mb = mb.N_w(Dg, Wc_mb) / n_at_ref[:, None] if sfd_name == "mainbelt": ratio = ratio_mb n_lo_mb = mb.N_w(np.array([D_REAL]), Wc_mb)[:, 0] r_real = n_lo_mb / n_at_ref cnt_ref = np.maximum(np.searchsorted(-mb.d, -d_ref, side="right"), 1) # bodies >= d_ref n_cnt = mb.count_w(Wc_mb, cnt_ref) mass_ref = mb.mass_w(w_mb, cnt_ref) / n_cnt mass_small = (mb.mass_w(w_mb, cnt_ref) - mb.mass_w(w_mb, np.full(n_draw, mb.n_real))) / n_cnt d_max = np.full(n_draw, mb.dmax) elif sfd_name in ("lunar_pl300", "lunar_pl600"): d_max = rng.uniform(dmax_range[0], dmax_range[1], n_draw) ratio = pl_ratio(Dg[None, :], d_ref[:, None], q_lunar[:, None], d_max[:, None]) ratio = np.where(Dg[None, :] < d_ref[:, None], ratio_mb, ratio) # below d_ref: main-belt shape r_real = pl_ratio(D_REAL, d_ref, q_lunar, d_max) mass_ref = pl_mass_per_ref(d_ref, q_lunar, d_max) mass_small = pl_mass_per_ref(d_ref, q_lunar, d_max, D_REAL) else: raise ValueError(sfd_name) A = R[:, None] * ratio # (n_draw, nD) # Lunar consistency: expected number of lunar impactors >= D_MOON_CHECK (clearly larger than the # SPA impactor, 172 km) since SPA formed. The Moon shows none. if np.any(np.abs(Dg - D_MOON_CHECK) < 1e-9): c_spa = np.take_along_axis(C, np.clip(np.searchsorted(Tg, t_spa), 0, nT - 1)[:, None], 1)[:, 0] lam_moon = ratio[:, int(np.argmin(np.abs(Dg - D_MOON_CHECK)))] * c_spa else: lam_moon = np.full(n_draw, np.nan) # ---- realisations of the large impactors: sizes -> realised late-accreted mass ---------- C_tot = C[:, -1] G = C / C_tot[:, None] # age distribution of an impact n_ref_earth = R * C_tot # Earth impactors >= d_ref since 4.51 Ga mass_exp = n_ref_earth * mass_ref n_big = rng.poisson(R * r_real * C_tot) idx = np.repeat(np.arange(n_draw), n_big) u = rng.uniform(size=idx.size) if sfd_name == "mainbelt": d_big = mb.sample_w(u, idx, Wc_mb, n_lo_mb) else: d_big = pl_sample(u, D_REAL, q_lunar[idx], d_max[idx]) m_big = np.pi / 6.0 * cfg.RHO_IMP * 1e9 * d_big ** 3 mass = n_ref_earth * mass_small + np.bincount(idx, weights=m_big, minlength=n_draw) ok_by = {c: (mass >= lo) & (mass <= hi) for c, (lo, hi) in conds.items()} Wm = np.stack([ok_by[c] / max(int(ok_by[c].sum()), 1) for c in conds]) # (n_cond, n_draw) j_real = int(np.argmin(np.abs(Dg - D_REAL))) assert abs(Dg[j_real] - D_REAL) < 1e-9, "d_grid must contain D_REAL" n_ge = np.stack([np.bincount(idx[d_big >= D], minlength=n_draw) if D >= D_REAL else np.zeros(n_draw, int) for D in Dg], axis=1) # (n_draw, nD) # ---- survival functions ------------------------------------------------------------------- nD = len(Dg) surv = np.zeros((nD, nT)) surv_s = np.zeros((nD, nT)) mu_mean = np.zeros((nD, nT)) surv_by = {c: np.zeros((nD, nT)) for c in conds} mu_by = {c: np.zeros((nD, nT)) for c in conds} log1mG = np.log(np.clip(1.0 - G, 1e-300, None)) for j in range(nD): mu = A[:, j, None] * C surv[j] = np.exp(-mu).mean(axis=0) mu_mean[j] = mu.mean(axis=0) if Dg[j] >= D_REAL: # given the sampled sizes: (1 - G)^n e = np.exp(n_ge[:, j, None] * log1mG) mu_r = n_ge[:, j, None] * G else: # bodies below D_REAL are not sampled: Poisson part + sampled part mu_lo = np.clip(A[:, j] - A[:, j_real], 0.0, None)[:, None] * C e = np.exp(-mu_lo + n_ge[:, j_real, None] * log1mG) mu_r = mu_lo + n_ge[:, j_real, None] * G surv_s[j] = e.mean(axis=0) sc, mc = Wm @ e, Wm @ mu_r for ci, c in enumerate(conds): surv_by[c][j], mu_by[c][j] = sc[ci], mc[ci] out = dict(chron=chron_name, sfd=sfd_name, T=Tg, D=Dg, surv=surv, surv_s=surv_s, mu_mean=mu_mean, surv_by=surv_by, mu_by=mu_by, ok_by=ok_by, frac_by={c: float(ok_by[c].mean()) for c in conds}, mass=mass, mass_exp=mass_exp, n_ref_earth=n_ref_earth, R=R, d_ref=d_ref, q_lunar=q_lunar, t_spa=t_spa, n_basins=nb, post_n=post_n, lam_moon=lam_moon, d_max=d_max, K=K, n1_spa=n1_spa, n_ge=n_ge, C_at=dict((float(t), C[:, int(np.searchsorted(Tg, t))]) for t in (3.5, 3.8, 4.0, 4.2, 4.4, 4.5))) if return_draws: out.update(A=A, C=C, ages=ages, n1=n1) return out def last_quantile(res, D, q): """Age T such that P(last impact >= D is older than T, or never happened after 4.51 Ga) = q, i.e. the q-quantile (counted from old to young) of the date of the last impact, from res["surv"]. Returns nan when even at 4.51 Ga the survival probability exceeds q (no such impact at all).""" j = int(np.argmin(np.abs(res["D"] - D))) s = res["surv"][j] # decreasing in T if s[-1] > q: return float("nan") if s[0] < q: return float(res["T"][0]) return float(np.interp(q, s[::-1], res["T"][::-1])) def view(res, cond): """A model result seen under one HSE treatment: surv and mu_mean are the conditional ones.""" return dict(res, surv=res["surv_by"][cond], mu_mean=res["mu_by"][cond], weight=res["frac_by"][cond]) def run_chunked(chron_name, sfd_name, n_draw, rng, chunk=2000, **kw): """run_model in chunks (bounded memory); concatenates per-draw diagnostics and averages the survival functions. se_by[c] is the Monte Carlo standard error of surv_by[c] from the scatter between chunks.""" parts = [] left = n_draw while left > 0: n = min(chunk, left) parts.append(run_model(chron_name, sfd_name, n, rng, **kw)) left -= n out = dict(parts[0]) conds = list(parts[0]["surv_by"]) wts = np.array([len(p["mass"]) for p in parts], float) for k in ("surv", "surv_s", "mu_mean"): out[k] = sum(p[k] * w for p, w in zip(parts, wts)) / wts.sum() out["surv_by"], out["mu_by"], out["se_by"], out["ok_by"], out["frac_by"] = {}, {}, {}, {}, {} for c in conds: oks = np.array([p["ok_by"][c].sum() for p in parts], float) if oks.sum() > 0: out["surv_by"][c] = sum(p["surv_by"][c] * o for p, o in zip(parts, oks)) / oks.sum() out["mu_by"][c] = sum(p["mu_by"][c] * o for p, o in zip(parts, oks)) / oks.sum() else: out["surv_by"][c] = np.full_like(out["surv"], np.nan) out["mu_by"][c] = np.full_like(out["surv"], np.nan) use = [p["surv_by"][c] for p, o in zip(parts, oks) if o > 0] out["se_by"][c] = (np.std(use, axis=0, ddof=1) / np.sqrt(len(use)) if len(use) > 1 else np.full_like(out["surv"], np.nan)) out["ok_by"][c] = np.concatenate([p["ok_by"][c] for p in parts]) out["frac_by"][c] = float(out["ok_by"][c].mean()) for k in ("mass", "mass_exp", "n_ref_earth", "R", "d_ref", "q_lunar", "t_spa", "d_max", "K", "n1_spa", "n_ge", "post_n", "lam_moon"): out[k] = np.concatenate([p[k] for p in parts]) out["C_at"] = {t: np.concatenate([p["C_at"][t] for p in parts]) for t in parts[0]["C_at"]} return out def mixture(results, cond, label="mixture"): """Equal-prior mixture of several model results under HSE treatment cond. Each model enters with its conditional survival function and a weight equal to the fraction of its realisations that satisfy the mass bounds (prior 1/n per model, likelihood = that fraction). With cond = "none" every model has weight 1/n.""" w = np.array([r["frac_by"][cond] for r in results]) s = sum(r["surv_by"][cond] * wi for r, wi in zip(results, w) if wi > 0) / w.sum() mu = sum(r["mu_by"][cond] * wi for r, wi in zip(results, w) if wi > 0) / w.sum() return dict(chron=label, sfd=label, cond=cond, T=results[0]["T"], D=results[0]["D"], surv=s, mu_mean=mu, weight=float(w.mean()), weights=w / w.sum()) def pareto_fit(x, xmin): """Maximum-likelihood cumulative slope q of N(>x) ~ x^-q for x >= xmin, its standard error, the KS distance and a parametric-bootstrap p-value (slope refitted in every replicate).""" from scipy import stats x = np.sort(np.asarray(x, float)) x = x[x >= xmin] n = len(x) q = n / np.log(x / xmin).sum() dist = stats.kstest(x, lambda v: 1.0 - (v / xmin) ** (-q)).statistic rng = np.random.default_rng(cfg.SEED + int(xmin)) worse, nrep = 0, 2000 for _ in range(nrep): y = xmin * rng.uniform(size=n) ** (-1.0 / q) qy = n / np.log(y / xmin).sum() worse += stats.kstest(y, lambda v: 1.0 - (v / xmin) ** (-qy)).statistic >= dist return dict(n=n, q=q, q_err=q / np.sqrt(n), ks=dist, p=worse / nrep)