"""Fast-versus-slow abiogenesis Bayes factor with an uncertain start of the habitable window. Notation (ages in Ga, durations in Gyr) t_w age at which liquid water is first present (4.404 Ga in Kipping 2020/2025) t_h age at which the window relevant to our lineage opens = min(t_w, age of the last sterilising impact). Its distribution comes from the bombardment model: P(t_h > T) = S(T) = P(no impactor >= D younger than T) for T < t_w, with an atom S(t_w) at t_h = t_w. t_ev age of the earliest evidence of life (abiogenesis happened before it) f remaining habitable time; the window measured from t_h has length T = t_h + f m number of hard steps that follow abiogenesis (m = 1 is Kipping's limiting case) Likelihoods are conditioned on success, as in the timing study. For a given t_h and all later rates -> 0 (hard steps), up to common factors, fast abiogenesis: numerator Nf = t_h^(m-1)/(m-1)! success Sf = T^m/m! slow abiogenesis: numerator Ns = [t_h^m - t_ev^m]/m! success Ss = T^(m+1)/(m+1)! (for m = 0: Nf = 1, Sf = 1, Ns = t_h - t_ev, Ss = T). Life evidence older than the last sterilising impact is impossible, so both numerators carry the indicator 1[t_h > t_ev]. t_h is a property of the planet, so it is integrated out of numerators and success probabilities separately: B = ( E[Nf] / E[Sf] ) / ( E[Ns] / E[Ss] ). With t_h fixed at t_w this is the timing study's B(m) = T m t_obs^(m-1) / ((m+1)[t_obs^m - (t_obs - t_L')^m]). """ from __future__ import annotations from math import factorial import numpy as np def th_distribution(Tg, S, t_w): """Discretise the distribution of t_h. Returns (points, masses): interval midpoints below t_w carrying mass S(T_i) - S(T_{i+1}), plus the atom S(t_w) at t_w. Masses sum to 1.""" Tg = np.asarray(Tg, float) S = np.asarray(S, float) keep = Tg < t_w T = np.concatenate([Tg[keep], [t_w]]) s = np.concatenate([S[keep], [np.interp(t_w, Tg, S)]]) pts = np.concatenate([0.5 * (T[1:] + T[:-1]), [t_w]]) mass = np.concatenate([s[:-1] - s[1:], [s[-1]]]) return pts, np.clip(mass, 0.0, None) def _terms(th, t_ev, f, m): T = th + f alive = th > t_ev if m == 0: return alive * 1.0, np.ones_like(th), alive * (th - t_ev), T Nf = alive * th ** (m - 1) / factorial(m - 1) Sf = T ** m / factorial(m) Ns = alive * (th ** m - t_ev ** m) / factorial(m) Ss = T ** (m + 1) / factorial(m + 1) return Nf, Sf, Ns, Ss def bf_sterilising(Tg, S, t_w, t_ev, f, m, convention="generative"): """Bayes factor fast:slow when impacts above the threshold reset the clock. convention 'generative': ratio of expectations (t_h integrated out of numerator and of the success probability separately). 'pointwise': the conditioned likelihood is averaged over t_h, the convention the timing study used for date uncertainties.""" th, p = th_distribution(Tg, S, t_w) Nf, Sf, Ns, Ss = _terms(th, t_ev, f, m) if convention == "generative": return float((p @ Nf) / (p @ Sf) / ((p @ Ns) / (p @ Ss))) return float((p @ (Nf / Sf)) / (p @ (Ns / Ss))) def bf_fixed(t_w, t_ev, f, m): """Closed form with the window opening at t_w (no impacts): the timing study's B(m).""" Nf, Sf, Ns, Ss = _terms(np.array([t_w]), t_ev, f, m) return float(Nf[0] / Sf[0] / (Ns[0] / Ss[0])) def window_stats(Tg, S, t_w, t_ev): """P(t_h > t_ev) and the mean and median window t_h - t_ev given t_h > t_ev.""" th, p = th_distribution(Tg, S, t_w) a = th > t_ev pa = float(p[a].sum()) if pa <= 0: return dict(p_ok=0.0, mean_window=float("nan"), median_window=float("nan"), p_atom=float(p[-1])) w = (th[a] - t_ev) c = np.cumsum(p[a]) / pa return dict(p_ok=pa, mean_window=float((p[a] @ w) / pa), median_window=float(np.interp(0.5, c, w)), p_atom=float(p[-1])) def th_full(Tg, S): """Distribution of the age of the last impact above the threshold over the whole grid (up to the Moon-forming impact, 4.51 Ga). The atom S(4.51) = "no such impact since the Moon-forming impact" is placed at the last grid age, i.e. the Moon-forming impact counts as the last one.""" Tg = np.asarray(Tg, float) S = np.asarray(S, float) pts = np.concatenate([0.5 * (Tg[1:] + Tg[:-1]), [Tg[-1]]]) mass = np.concatenate([S[:-1] - S[1:], [S[-1]]]) return pts, np.clip(mass, 0.0, None) def bf_enabled(Tg, S, t_w, t_ev, f, m, W): """Impact-enabled origin. Abiogenesis is possible only while the transient reducing atmosphere made by the last impact above the threshold persists, i.e. at ages between t_i and t_i - W (t_i = age of that impact, which also sterilises), and only once liquid water exists (age < t_w). Within that window the rate is constant; outside it is zero. opening age s = min(t_i, t_w); closing age c = max(t_i - W, -f) a planet with c >= t_w never has an origin and drops out of the success-conditioned average fast: Nf = 1[s > t_ev] s^(m-1)/(m-1)! Sf = (s + f)^m / m! slow: Ns = [s^m - max(t_ev, c)^m]_+ / m! Ss = [(s + f)^(m+1) - (c + f)^(m+1)] / (m+1)! If the evidence date is younger than the closing age, early life is certain under both hypotheses and B = 1; W -> infinity recovers the sterilising case.""" ti, p = th_full(Tg, S) s = np.minimum(ti, t_w) c = np.maximum(ti - W, -f) live = c < t_w # the planet has some window with liquid water alive = live & (s > t_ev) lo = np.maximum(t_ev, c) if m == 0: Nf = alive * 1.0 Sf = live * 1.0 Ns = alive * np.clip(s - lo, 0.0, None) Ss = live * (s - c) else: Nf = alive * s ** (m - 1) / factorial(m - 1) Sf = live * (s + f) ** m / factorial(m) Ns = alive * np.clip(s ** m - np.sign(lo) * np.abs(lo) ** m, 0.0, None) / factorial(m) Ss = live * ((s + f) ** (m + 1) - (c + f) ** (m + 1)) / factorial(m + 1) if (p @ Ns) <= 0 or (p @ Nf) <= 0: return float("nan") return float((p @ Nf) / (p @ Sf) / ((p @ Ns) / (p @ Ss))) def likelihoods(Tg, S, t_w, t_ev, f, m): """Success-conditioned likelihoods of 'life by t_ev' under the four joint hypotheses {fast, slow abiogenesis} x {impacts above the threshold sterilise, life survives them}. The slow likelihoods are per unit abiogenesis rate (the common factor lambda cancels in every ratio taken between the two slow entries or in a fast:slow ratio at fixed prior).""" th, p = th_distribution(Tg, S, t_w) Nf, Sf, Ns, Ss = _terms(th, t_ev, f, m) nf0, sf0, ns0, ss0 = (x[0] if np.ndim(x) else x for x in _terms(np.array([t_w]), t_ev, f, m)) return dict(Lf_ster=float((p @ Nf) / (p @ Sf)), Ls_ster=float((p @ Ns) / (p @ Ss)), Lf_surv=float(nf0 / sf0), Ls_surv=float(ns0 / ss0)) def bf_mixed(L, prior_ster): """Fast:slow Bayes factor when the sterilise/survive question is itself uncertain, with prior probability prior_ster that impacts above the threshold sterilise.""" num = prior_ster * L["Lf_ster"] + (1 - prior_ster) * L["Lf_surv"] den = prior_ster * L["Ls_ster"] + (1 - prior_ster) * L["Ls_surv"] return float(num / den) def bf_date_marginal(Tg, S, m, luca_ages, t_w, f): """Date-marginalised B for the LUCA case, following the timing study: the (already t_h- marginalised) conditioned likelihoods are averaged over the LUCA ages supplied. Ages within 5 Myr of t_w or older are dropped, as in that study. Vectorised over the ages.""" th, p = th_distribution(Tg, S, t_w) # th ascending T = th + f t_ev = np.asarray(luca_ages, float) t_ev = t_ev[t_ev < t_w - 0.005] k = np.searchsorted(th, t_ev, side="right") # th[k:] > t_ev def tail(v): # sum of v over th > t_ev c = np.concatenate([np.cumsum(v[::-1])[::-1], [0.0]]) return c[k] if m == 0: Nf, Sf = tail(p), 1.0 Ns, Ss = tail(p * th) - t_ev * tail(p), float(p @ T) else: Nf = tail(p * th ** (m - 1)) / factorial(m - 1) Sf = float(p @ T ** m) / factorial(m) Ns = (tail(p * th ** m) - t_ev ** m * tail(p)) / factorial(m) Ss = float(p @ T ** (m + 1)) / factorial(m + 1) return float((Nf / Sf).sum() / (Ns / Ss).sum()), int(t_ev.size)