"""
Integrated Bayesian model of the Fermi observation, version 4.1 (rebuilt after red-team round 2; corrected after round 3).

Chain:  habitable planets (three host classes) -> life (rate-based) -> hard steps or an environmental gate
        -> technological civilisation -> sustained expansion (explicit universal-barrier atom + continuous part)
        -> is that expansion evident to us? (decided per expander lineage, with two universe-level exceptions)

Every sample is one possible universe. Its weight is the probability of what we actually observe.

What changed from version 3, and which round-2 objection each change answers:
  * Evidentness is decided per lineage. A universe in which every viable form of expansion is low-footprint is a
    small explicit prior atom (q_cc); elsewhere a fraction c of lineages claim systems and each of the others has
    its own chance of being evident. Many independent lineages can no longer all go unnoticed for free.
  * The 30-decade tail of the expansion prior is replaced by an explicit atom (q_barrier: sustained interstellar
    expansion is infeasible for everyone, us included) plus a continuous part with the hard-step width.
  * w_fast is a pre-data prior (one log-uniform prior on the abiogenesis rate), so early life is counted once.
  * The biology chain is exact (phase-type completion probabilities, chain4.py), for the window total and for the
    count of civilisations to date. Worlds whose later steps are all easy are treated like worlds with no steps.
  * Worlds with no hard step have an explicit gate clock, so their timing likelihood is a real density.
  * Long-lived hosts no longer get time for free: persistence can act as a constant hazard, K-dwarf windows are
    habitable-zone residence times, and prebiotic ultraviolet chemistry is penalised on K and M dwarfs.
  * The head start of earlier civilisations comes from the model's own birth tables (it was counted twice).
  * Claimed galaxies may contain unclaimed pockets; pocket civilisations are in the SSA reference class.
  * The same survival function is used in the numerator and denominator of every reference-class column.
  * Infrared survey limits are split into three temperature bands with the curves each survey supports.
  * Early-life and late-arrival terms can be switched separately; future-observation error rates are parameters.

Columns (anthropic conventions), all from the same samples:
  site      Earth is a given site (numerically: SSA among Sun-like hosts at our epoch)
  ssa_all   we are a random draw from every civilisation that ever arises, past and future, on any host
  ssa_now   we are a random draw from the civilisations that have arisen to date
  ssa_obs   as ssa_all, but drawing observers, with sustained expanders holding obs_R times more observers
  sia       universes are also weighted by how many civilisations in our situation they contain

Version 4.1 (round 3; none of these moves a headline number by more than 0.005):
  * Prior class masses and prior statistics use an unconditioned gate delay (the posterior draws it where Earth's
    record allows); pool() averages the prior statistics over seeds.
  * One run can return every subset of the four observation terms (ablate=True) and any list of hypothetical future
    observations (futures=[...]), so the decomposition and the forecasts are pooled over the same seeds as the base.
  * "Expansion shown feasible" conditions on no universal barrier; "shown to be common" is the old p_exp > 1e-3 row.
  * The chance that anyone else is in the Galaxy now counts settlers as well as home-grown civilisations.
  * Forecasts that include an assumed false-alarm rate also store the false-alarm-free probability.
  * Chain clipping widened; far-infrared survey curve updated to the revised study and made cumulative; formation
    weights normalised to stars formed to date.
  * early_bf: optional fixed odds for a fast origin in place of the early-life likelihood (opposing-camp rows).

Nothing in this file is tuned to an answer; every prior is listed in BASE with its source or marked as a judgement.
"""
import copy
import math
import os
import numpy as np
from scipy.special import gammainc, gamma as gamma_fn

import chain4 as EC
import cosmo_tables as CT

HERE = os.path.dirname(os.path.abspath(__file__))
TAB = dict(np.load(os.path.join(HERE, "cosmo_tables.npz")))
T_GRID, ETA_GRID, I0 = TAB["t"], TAB["eta"], int(TAB["i0"])
L_GRID = TAB["L_grid"]
T0 = float(T_GRID[I0])
RHO_GAL = 2.88e5          # Milky-Way equivalents per Gly^3
T_WIN = 5.6               # Gyr, reference habitable window of an Earth-like planet around a G star
_XG, _WG = np.polynomial.legendre.leggauss(40)
_UG, _WG = 0.5 * (_XG + 1), 0.5 * _WG
# formation-time nodes for the count of civilisations to date
_TF = np.linspace(0.05, T0, 49)
_TFM = 0.5 * (_TF[1:] + _TF[:-1])
_WF = CT.sfr(_TFM) * np.diff(_TF)
_WF = _WF / _WF.sum()                                           # fraction of today's planets formed at each node (Nstars is a present-day count)
_AGE = T0 - _TFM

BASE = dict(
    n=400_000, seed=7,
    tilt=False,                  # importance-sample the efficiency factors (needed only for a converged SIA column)
    # ---- habitable planets per Milky-Way-like galaxy (evidence/astro.md; gaps_and_corrections.md) ----
    log10_Nstars=(11.0, 11.6),
    f_G=(0.10, 0.18),            # F, G and early-K stars (>~0.75 Msun), uniform
    r_K=(0.4, 0.9),              # long-lived K dwarfs (0.6-0.75 Msun) per G-class star (Kroupa IMF gives ~0.66)
    eta_earth=(0.05, 0.9),       # log-uniform, rocky habitable-zone planets per star
    f_site=(0.01, 1.0),          # log-uniform, such a planet is actually habitable for Gyr-scale biology
    p_M_class=0.30,              # P(M-dwarf planets can host biospheres at all). Judgement: the M-dwarf study leaves atmosphere retention unconstrained
    h_M=(0.01, 1.0),             # if so: per-planet habitability relative to a G-star planet, log-uniform
    f_Mdwarf=0.75,
    L_K=(20.0, 45.0),            # Gyr, K-dwarf main-sequence lifetime
    L_K_hz=(0.6, 0.9),           # fraction of that lifetime a planet spends in the habitable zone (round 2, astrobiologist)
    L_M=(50.0, 1000.0),          # Gyr, log-uniform
    p_cap_short=0.7,             # P(planets themselves stop being habitable after ~6-15 Gyr, whatever the star)   [see studies/planet_lifetime]
    L_cap_short=(6.0, 15.0), L_cap_long=(15.0, 100.0),
    K_uv=(0.1, 1.0),             # factor on the abiogenesis rate around K dwarfs (less near-ultraviolet for prebiotic chemistry), log-uniform
    M_uv=(0.01, 1.0),            # same for M dwarfs
    # ---- abiogenesis: expected independent origins per habitable window, lambda ----
    w_fast=0.2,                  # prior P(lambda >= 1). Pre-data: one log-uniform prior from 3 decades above 1 to A decades below gives E[3/(A+3)] = 0.19 (0.16-0.40 depending on the cut-offs; 0.2 is a judgement)
    fast_decades=3.0, fast_lo=0.0,   # fast camp: log10 lambda ~ U(fast_lo, fast_decades)
    early_bf=None,               # if set: replace the early-life likelihood by these fixed odds in favour of the fast camp (opposing-camp rows)
    slow_decades=(3.0, 100.0),   # slow camp: log10 lambda ~ U(-A, 0), A log-uniform on this range
    tau_over_T=(0.02, 0.20),     # delay from habitable to inhabited / window (spans LUCA and fossils-only datings)   [STUDY timing]
    # ---- life -> technological civilisation ----
    m_weights={0: 0.30, 1: 0.30, 2: 0.22, 3: 0.12, 4: 0.06},   # candidate steps after abiogenesis. Judgement; the convergence study leaves 0-4 open
    step_decades=(2.0, 8.0),     # each step: log10 p ~ U(-D, 0) per window; D log-uniform
    hard_thr=0.3,                # a step counts as hard if its per-window probability is below this. The yields are exact whatever the label, but the
                                 # label decides which worlds get the gate clock (and f_gate) and sets the exponent in the timeline likelihood.
                                 # Measured sensitivity of P(alone) to thresholds 0.1-0.5: 0.004 or less (round 3 code audit)
    gate_delay=(0.1, 0.95),      # worlds with no hard step: delay from life to civilisation, in windows, uniform (an environmental clock)
    f_gate=(0.03, 1.0),          # those worlds: chance the environmental gates open at all, log-uniform
    f_persist=(0.01, 1.0),       # log-uniform, climate and biosphere persist over one reference window (Tyrrell 2020: 1.45% of runs stayed habitable for 3 Gyr)
    p_persist_hazard=0.5,        # P(persistence is a constant hazard through time) versus a fixed fraction of intrinsically stable planets
    f_c=(0.01, 1.0),             # log-uniform, intelligence -> technological civilisation
    x_now=(0.70, 0.85),          # time of our arrival / habitable window
    use_early=True, use_late=True,
    # ---- sustained expansion: per-civilisation probability of becoming a persistent, self-propagating expander ----
    q_barrier=0.15,              # P(sustained interstellar expansion is infeasible for everyone, us included); range 0.05-0.30. My judgement: studies/expansion_analogs looked for an empirical basis and found none (no prior recommended)
    W_exp=(1.0, 8.0),            # otherwise log10 p_exp ~ U(-W, 0), W log-uniform: the same width as a hard step. p_exp is the population mean
    v_gal=(1e-4, 0.3),           # settlement-front speed inside a galaxy / c, log-uniform (sets the fill time)
    q_ig=0.6,                    # P(expanders also cross between galaxies)
    v_ig=("loguniform", 0.01, 0.99),   # their speed between galaxies / c; ("uniform", 0.1, 0.99) for the expansionist camp (round 2, Hansonian)
    # ---- would expansion be evident here? ----
    p_stealth=0.03,              # P(every expander deliberately hides and leaves living planets alone); needs enforcement, so universe-level
    q_cc=0.05,                   # P(every viable form of expansion is low-footprint, by physics or economics); range 0.02-0.10
    p_here=(1e-3, 0.3),          # chance a low-footprint presence is nevertheless evident to us, log-uniform (per lineage, or per Galaxy in a q_cc universe)
    c_claim=("uniform", 0.0, 1.0),   # elsewhere: fraction of lineages that claim systems and are evident. Uniform is neutral; ("loguniform", lo, hi) for the sceptic camp
    eps_pocket=(1e-6, 0.1),      # chance a site in a claimed Galaxy is unclaimed and sees nothing, log-uniform (round 2, expansion sceptic)
    p_harvest=0.5,               # P(claiming settlers intercept more than 1% of starlight)
    f_cov=(0.01, 1.0),           # covering fraction if harvested, log-uniform
    p_warm=0.6,                  # P(waste heat emerges at 100-600 K, where WISE looks)
    p_verycold=0.15,             # P(waste heat emerges just above the CMB, where Planck looks); the rest is far-infrared (10-100 K)
    warm_curve=((0.25, 115.0), (0.5, 4.8e3), (0.85, 1e5)),     # covering fraction -> effective galaxies excluded (Curtis et al. 2026 preprint; Griffith et al. 2015)
    mid_curve=((0.75, 117.0), (0.9, 195.0)),                   # far-infrared: optical scaling relations, disc galaxies only   [STUDY tf_kiii, revised: 2.6% and 1.5%]; cumulative
    verycold_curve=((0.01, 3e3), (0.1, 1e5), (0.85, 1e6)),     # Lacki 2016 (Planck; unrefereed preprint)
    use_local=True, use_surveys=True, use_cosmic=True, use_seti=True,
    # ---- quiet civilisations and the technosignature record ----
    L_quiet=(1e2, 1e9),          # yr detectable lifetime, log-uniform
    f_beacon=(1e-4, 0.3),        # fraction of extant civilisations running an always-on beacon, log-uniform
    d_beacon=(1e-6, 0.5),        # chance the whole search record would have caught one such beacon (1e17 W to Kardashev II)  [STUDY seti_meta]
    c_leak=(1e-17, 1e-11),       # the same for a civilisation that only leaks, as we do  [STUDY seti_meta]
    N_reach=5.5e6,               # galaxies within the volume a 0.1c expansion starting now could ever reach (cosmo_tables)
    log10_G=(10.7, 11.3),        # Milky-Way equivalents in the observable universe
    obs_R=1e6,                   # ssa_obs: observers in a sustained expander relative to a civilisation like ours
    # ---- hypothetical future observations ----
    extra_obs=None,
    ablate=False,                # also return every subset of the four observation terms (Earth-as-given column)
    futures=None,                # list of (key, extra_obs, overrides): hypothetical observations evaluated on the same samples
    sg_fp=0.15, sg_frac=0.3,     # second origin in the Solar System: false-positive rate (0.05-0.3) and share of a planet's opportunity sampled
    hwo_fp=0.10, hwo_q=0.5,      # 25 habitable-zone planets: survey-level false-positive rate (0.03-0.2) and detectability of a living planet
    art_d=0.5,                   # artefact survey: chance of finding the traces of a low-footprint presence, if there is one
)

CLASSES = ["empty_universe", "rare", "en_route", "stuck", "unnoticed", "hidden", "excluded"]
CLASS_LABELS = {
    "empty_universe": "No other civilisation has ever arisen in the observable universe",
    "rare": "Others exist in the universe, none other has ever arisen in our Galaxy, nobody is on the way",
    "en_route": "Expanding civilisations already hold over 1% of the universe and have not reached us",
    "stuck": "Civilisations have arisen in our Galaxy and none became a sustained expander",
    "unnoticed": "The Galaxy has been settled by low-footprint expanders and it is not evident from here",
    "hidden": "The Galaxy has been settled and they are deliberately hiding",
    "excluded": "The Galaxy is claimed by evident settlers and we sit in an unclaimed pocket (contradicted unless pockets exist)",
}
SCHOOLS = ["site", "ssa_all", "ssa_now", "ssa_obs", "sia"]
SCHOOL_LABELS = {
    "site": "Earth as a given (SSA among Sun-like hosts)",
    "ssa_all": "Typical among every civilisation, past and future",
    "ssa_now": "Typical among civilisations to date",
    "ssa_obs": "Typical observer (expanders hold most observers)",
    "sia": "Crowded universes count more (SIA)",
}


def lu(rng, lo, hi, n):
    return 10 ** rng.uniform(np.log10(lo), np.log10(hi), n)


def interp_curve(f, curve, zero_below=True, zero_above=False):
    xs = np.log10([c[0] for c in curve]); ys = np.log10([c[1] for c in curve])
    out = 10 ** np.interp(np.log10(np.maximum(f, 1e-9)), xs, ys)
    out = np.where(f < curve[0][0], 0.0, out) if zero_below else out
    out = np.where(f > curve[-1][0], 0.0, out) if zero_above else out
    return out


def vary(base, **kw):
    P = copy.deepcopy(base); P.update(kw); return P


def _tab_at(table, nidx, lidx, tq):
    """table[nidx, lidx] interpolated linearly at cosmic times tq (n,) -> (3, n); zero before the first grid time."""
    j = np.clip(np.searchsorted(T_GRID, tq), 1, len(T_GRID) - 1)
    w = np.clip((tq - T_GRID[j - 1]) / (T_GRID[j] - T_GRID[j - 1]), 0.0, 1.0)
    out = (1 - w) * table[nidx, lidx, j - 1] + w * table[nidx, lidx, j]
    return np.where(tq <= T_GRID[0], 0.0, out)


# --------------------------------------------------------------------------------------------------------------
# biology: exact yields per host class, and the likelihood of Earth's own timeline
# --------------------------------------------------------------------------------------------------------------
def _gated_cdf(lam, h, D, s):
    """P(civilisation by time s) with no hard step: life at rate lam, then a delay D, with hazard h throughout."""
    a = np.maximum(s - D, 0.0)
    r = lam + h
    frac = np.where(r * a < 1e-8, lam * a, lam / np.maximum(r, 1e-300) * -np.expm1(-r * a))
    return np.exp(-h * D) * frac


def yields_hard(lam3, rates, m, gated, h, L3):
    """lam3, L3: (3, n) abiogenesis rate per window and habitable window (Gyr) per host class.
    Returns log10 P(civilisation within the window) (3, n) and the sum over formation times of
    P(civilisation by the planet's age today) (3, n); both without the constant factors f_gate, f_c, f_persist.
    Worlds with no hard step (gated) are left unfilled: see yields_gated."""
    n = len(m)
    logPw = np.full((3, n), -300.0)
    now = np.zeros((3, n))
    for mm in np.unique(m[~gated]):
        sel = np.where((m == mm) & ~gated)[0]
        for a0 in range(0, len(sel), 50000):
            ii = sel[a0:a0 + 50000]
            for c in range(3):
                lp = EC.logcdf_grid(np.c_[lam3[c, ii], rates[ii, :mm]], h[ii])
                logPw[c, ii] = EC.interp_log(lp, L3[c, ii] / T_WIN)
                acc = np.zeros(len(ii))
                for k in range(len(_AGE)):
                    acc += _WF[k] * 10 ** EC.interp_log(lp, np.minimum(_AGE[k], L3[c, ii]) / T_WIN)
                now[c, ii] = acc
    return logPw, now


def yields_gated(logPw, now, lam3, gated, h, D, L3):
    """Fill the gated worlds of copies of (logPw, now) for a given gate delay D."""
    logPw = logPw.copy(); now = now.copy()
    g = np.where(gated)[0]
    for c in range(3):
        s = L3[c, g] / T_WIN
        logPw[c, g] = np.log10(np.maximum(_gated_cdf(lam3[c, g], h[g], D[g], s), 1e-300))
        acc = np.zeros(len(g))
        for k in range(len(_AGE)):
            a = np.minimum(_AGE[k], L3[c, g]) / T_WIN
            acc += _WF[k] * _gated_cdf(lam3[c, g], h[g], D[g], a)
        now[c, g] = acc
    return logPw, now


def _int_life(lam, a, f):
    """integral over [0, a] of lam * exp(-lam t) * f(t) dt, by Gauss-Legendre in the cumulative variable (exact for any lam)."""
    ea = -np.expm1(-lam * a)
    small = lam * a < 1e-6
    t = np.where(small[:, None], _UG[None, :] * a[:, None],
                 -np.log1p(-_UG[None, :] * np.minimum(ea, 1 - 1e-15)[:, None]) / lam[:, None])
    return np.where(small, lam * a, ea) * (f(t) * _WG[None, :]).sum(1)


def timeline(lam, h, k, gated, D, w_D, x, tau, use_early, use_late):
    """Likelihood of Earth's timeline for a Sun-like host: life within tau of habitability (early) and our own
    arrival at fraction x of the window (late), given that a civilisation arises at all.
    Hard-step worlds: abiogenesis at rate lam, then k hard steps (density u^(k-1)), hazard h.
    Gated worlds: abiogenesis, then the gate delay D (w_D is the importance weight of the way D was drawn)."""
    n = len(lam)
    one = np.ones(n)
    if not (use_early or use_late):
        return one
    lam = np.maximum(lam, 1e-12)               # below this the likelihood is already in its small-rate limit
    kk = np.maximum(k, 1).astype(float)[:, None]
    hh = h[:, None]

    def G(a):                                  # integral over [0, a] of u^(k-1) exp(-h u) du
        with np.errstate(divide="ignore", invalid="ignore"):
            big = gammainc(kk, np.maximum(hh * a, 0.0)) * gamma_fn(kk) / np.maximum(hh, 1e-12) ** kk
        return np.where(hh > 1e-6, big, np.maximum(a, 0.0) ** kk / kk)

    den = _int_life(lam, one, lambda t: np.exp(-hh * t) * G(1 - t))
    if use_early and use_late:
        num = np.exp(-h * x) * _int_life(lam, tau, lambda t: np.maximum(x[:, None] - t, 0.0) ** (kk - 1))
    elif use_late:
        num = np.exp(-h * x) * _int_life(lam, x, lambda t: np.maximum(x[:, None] - t, 0.0) ** (kk - 1))
    else:
        num = _int_life(lam, tau, lambda t: np.exp(-hh * t) * G(1 - t))
    L_hard = num / np.maximum(den, 1e-300)
    # gated worlds
    r = lam + h
    a1 = np.maximum(1 - D, 1e-9)
    inv = np.where(r * a1 < 1e-8, 1 / a1, r / np.maximum(-np.expm1(-r * a1), 1e-300))    # r / P(life in time)
    if use_late:                               # with use_early, D was drawn inside [x - tau, x] and w_D carries the prior mass of that interval
        L_gate = np.where(D <= x, inv * np.exp(-r * np.maximum(x - D, 0.0)), 0.0) * w_D
    else:
        L_gate = -np.expm1(-r * np.minimum(tau, a1)) / np.maximum(-np.expm1(-r * a1), 1e-300)
    return np.where(gated, L_gate, L_hard)


# --------------------------------------------------------------------------------------------------------------
# one run = n universes
# --------------------------------------------------------------------------------------------------------------
TERMS = ("early", "late", "local", "surveys")     # early life; our late arrival; nobody evidently here; galaxy surveys and radio searches
SPREAD = ("en_route", "unnoticed", "hidden", "excluded")


def subset_key(on):
    return "+".join(t for t in TERMS if t in on) or "prior"


def run(P, keep=False):
    rng = np.random.default_rng(P["seed"]); n = P["n"]
    one = np.ones(n)
    log_imp = np.zeros(n)                       # log(prior / proposal), non-zero only with tilt=True

    def lu_t(lo, hi, mask=None):
        """log-uniform prior; with tilt, drawn uniformly in the variable and re-weighted (cancels the factor in the SIA weight)."""
        nonlocal log_imp
        if not P["tilt"]:
            return lu(rng, lo, hi, n)
        xv = rng.uniform(lo, hi, n)
        li = np.log((hi - lo) / (xv * math.log(hi / lo)))
        log_imp = log_imp + (li if mask is None else np.where(mask, li, 0.0))
        return xv

    # ---------------- planets ----------------
    Nstars = 10 ** rng.uniform(*P["log10_Nstars"], n)
    eta_e = lu_t(*P["eta_earth"]); f_site = lu_t(*P["f_site"])
    NG = Nstars * rng.uniform(*P["f_G"], n) * eta_e * f_site
    NK = NG * rng.uniform(*P["r_K"], n)
    M_alive = rng.random(n) < P["p_M_class"]
    NM = np.where(M_alive, Nstars * P["f_Mdwarf"] * eta_e * f_site * lu(rng, *P["h_M"], n), 0.0)
    Ncls = np.vstack([NG, NK, NM])
    cap_short = rng.random(n) < P["p_cap_short"]
    L_cap = np.where(cap_short, rng.uniform(*P["L_cap_short"], n), lu(rng, *P["L_cap_long"], n))
    LK = np.minimum(rng.uniform(*P["L_K"], n) * rng.uniform(*P["L_K_hz"], n), L_cap)
    LM = np.minimum(lu(rng, *P["L_M"], n), L_cap)
    L3 = np.vstack([np.full(n, T_WIN), LK, LM])
    # ---------------- biology ----------------
    if P["tilt"]:
        fast = rng.random(n) < 0.9
        log_imp = log_imp + np.where(fast, math.log(P["w_fast"] / 0.9), math.log((1 - P["w_fast"]) / 0.1))
    else:
        fast = rng.random(n) < P["w_fast"]
    A = lu(rng, *P["slow_decades"], n)
    f_lo = P.get("fast_lo", 0.0)
    log_lam = np.where(fast, f_lo + (P["fast_decades"] - f_lo) * rng.random(n), -A * rng.random(n))
    lam = 10 ** np.clip(log_lam, -250, 10)
    lam3 = np.vstack([lam, lam * lu(rng, *P["K_uv"], n), lam * lu(rng, *P["M_uv"], n)])
    ms = np.array(sorted(P["m_weights"])); mw = np.array([P["m_weights"][k] for k in ms], float)
    m = rng.choice(ms, size=n, p=mw / mw.sum())
    Dd = lu(rng, *P["step_decades"], n)
    p_steps = 10 ** (-Dd[:, None] * rng.random((n, max(int(ms.max()), 1))))
    rates = -np.log1p(-np.minimum(p_steps, 1 - 1e-12))                     # per reference window
    act = np.arange(p_steps.shape[1])[None, :] < m[:, None]
    m_hard = (act & (p_steps < P["hard_thr"])).sum(1)                     # hard steps after abiogenesis
    gated = m_hard == 0                                                    # no hard step: an environmental clock sets the pace
    f_gate = np.where(gated, lu_t(*P["f_gate"], mask=gated), 1.0)
    f_pers = lu_t(*P["f_persist"]); f_c = lu_t(*P["f_c"])
    hazard = rng.random(n) < P["p_persist_hazard"]
    h = np.where(hazard, -np.log(f_pers), 0.0)                             # constant hazard per window, or a flat factor
    log_ext = np.log10(f_gate) + np.log10(f_c) + np.where(hazard, 0.0, np.log10(f_pers))
    # Earth's timeline (drawn before the yields because the gate delay enters both)
    x_now = rng.uniform(*P["x_now"], n); tau = rng.uniform(*P["tau_over_T"], n)
    ga, gb = P["gate_delay"]
    uD = rng.random(n)
    early_bf = P.get("early_bf")
    use_early = P["use_early"] and not early_bf
    cond = use_early and P["use_late"]
    D_prior = ga + uD * (gb - ga)               # the prior draw of the gate delay
    if cond:                                    # posterior: draw the delay where Earth's record allows it and carry the prior mass of that interval
        lo = np.maximum(ga, x_now - tau); hi = np.minimum(gb, x_now)
        w_D = np.maximum(hi - lo, 0.0) / (gb - ga)
        D = np.where(w_D > 0, lo + uD * (hi - lo), ga)
    else:
        w_D = one; D = D_prior
    hardPw, hardNow = yields_hard(lam3, rates, m, gated, h, L3)
    lidx = np.abs(np.log(L3[:, :, None] / L_GRID[None, None, :])).argmin(axis=2)
    life_worlds = (Ncls * -np.expm1(-lam3 * L3 / T_WIN)).sum(0)             # planets in the Galaxy on which life ever starts
    # ---------------- expansion ----------------
    barrier = rng.random(n) < P["q_barrier"]
    Wexp = lu(rng, *P["W_exp"], n)
    p_exp = np.where(barrier, 0.0, 10 ** (-Wexp * rng.random(n)))
    v = lu(rng, *P["v_gal"], n)
    kind_v, v_lo, v_hi = P["v_ig"]
    v_ig = rng.uniform(v_lo, v_hi, n) if kind_v == "uniform" else lu(rng, v_lo, v_hi, n)
    ig = rng.random(n) < P["q_ig"]
    # Gyr to settle and fill a galaxy, not the transit time: 1 Gyr at probe speed (1e-4 c); never under 0.3 Gyr however fast the ships,
    # because fill-in is paced by settlement and stellar mixing (Carroll-Nellenback et al. 2019; flyby study). Swapping the floor for a
    # cap changes P(alone) by 0.0002 (round 3).
    T_fill = np.maximum(0.3, 1e-4 / v)
    kv = RHO_GAL * v_ig ** 3 * p_exp * ig
    # ---------------- evidentness ----------------
    u = rng.random(n)
    stealth = u < P["p_stealth"]
    cc = (~stealth) & (u < P["p_stealth"] + P["q_cc"])
    mixed = ~(stealth | cc)
    kind, c_lo, c_hi = P["c_claim"]
    c_draw = rng.uniform(c_lo, c_hi, n) if kind == "uniform" else lu(rng, c_lo, c_hi, n)
    c = np.where(mixed, c_draw, 0.0)             # fraction of lineages that claim systems and are evident
    eps = np.where(mixed, lu(rng, *P["eps_pocket"], n), 0.0)
    p_here = np.where(stealth, 0.0, lu(rng, *P["p_here"], n))
    harv = mixed & (rng.random(n) < P["p_harvest"])
    f_cov = np.where(harv, lu(rng, *P["f_cov"], n), 0.0)
    ub = rng.random(n)
    band = np.where(ub < P["p_warm"], 0, np.where(ub < P["p_warm"] + P["p_verycold"], 2, 1))     # 0 warm, 1 intermediate, 2 very cold
    N_eff = np.where(band == 0, interp_curve(f_cov, P["warm_curve"]),
                     np.where(band == 1, interp_curve(f_cov, P["mid_curve"]), interp_curve(f_cov, P["verycold_curve"])))
    L_q = lu(rng, *P["L_quiet"], n)
    f_b = lu(rng, *P["f_beacon"], n)
    c_det = f_b * lu(rng, *P["d_beacon"], n) + (1 - f_b) * lu(rng, *P["c_leak"], n)
    G = 10 ** rng.uniform(*P["log10_G"], n)

    def derive(Dv):
        """Everything that follows from the biology for one choice of gate delay: counts, expander lineages, the three
        silence likelihoods and the classes."""
        logPw, now = yields_gated(hardPw, hardNow, lam3, gated, h, Dv, L3)
        logW = np.log10(np.maximum(Ncls, 1e-300)) + logPw + log_ext[None, :]
        W = np.where(Ncls > 0, 10 ** np.maximum(logW, -300.0), 0.0)            # civilisations per galaxy over each class's whole window
        f_now = np.clip(now / np.maximum(10 ** logPw, 1e-300), 0.0, 1.0)      # fraction of them born by today
        lf = np.log(np.maximum(f_now, 1e-300))
        nidx = np.zeros((3, n), int)
        best = np.full((3, n), np.inf)
        for k in range(len(TAB["n_steps"])):     # birth-date shape: the table row whose born-by-now fraction matches the exact one
            dd = np.abs(np.log(np.maximum(TAB["F"][k][lidx, I0], 1e-300)) - lf)
            upd = dd < best
            nidx[upd] = k; best[upd] = dd[upd]
        F0 = TAB["F"][nidx, lidx, I0]
        Wr = W * np.clip(f_now / np.maximum(F0, 1e-300), 0.2, 5.0)             # rescaled so that W_r * F0 is the exact count to date
        N_ever = (Wr * F0).sum(0)                # civilisations that have arisen by now, per Milky-Way-like galaxy
        N_alltime = (Wr * TAB["F"][nidx, lidx, -1]).sum(0)
        mu_loc = p_exp * (Wr * _tab_at(TAB["F"], nidx, lidx, T0 - T_fill)).sum(0)     # expanders born here early enough to have filled the Galaxy
        x = kv * (Wr * TAB["kern"][nidx, lidx, I0]).sum(0) if P["use_cosmic"] else np.zeros(n)
        mu = mu_loc + x                          # expected number of independent expander lineages whose domain includes the Galaxy today
        S_claim = np.exp(-np.minimum(mu * c, 700.0))                           # no claiming lineage has reached us
        L_claim = S_claim + (1 - S_claim) * eps
        L_low = np.where(mixed, np.exp(-np.minimum(mu * (1 - c) * p_here, 700.0)), 1 - (-np.expm1(-mu)) * p_here)
        N_now = N_ever * L_q / 1e10
        c_excl = mixed & (mu * c >= 1)
        c_unnot = (mixed & ~c_excl & (mu * (1 - c) >= 1)) | (cc & (mu >= 1))
        c_hidden = stealth & (mu >= 1)
        rest = ~(c_excl | c_unnot | c_hidden)
        c_route = rest & (x >= 0.01)
        rest2 = rest & ~c_route
        c_empty = rest2 & (N_ever * G < 1)
        c_rare = rest2 & ~c_empty & (N_ever < 1)
        c_stuck = rest2 & (N_ever >= 1)
        cls = np.full(n, -1)
        for i, mk in enumerate([c_empty, c_rare, c_route, c_stuck, c_unnot, c_hidden, c_excl]):
            cls[mk] = i
        assert (cls >= 0).all()
        alone = N_ever < 1
        nobody = c_empty | c_rare
        sub = dict(                                                             # finer splits of the headline classes
            alone_sterile=alone & (life_worlds < 1), alone_inhabited=alone & (life_worlds >= 1),
            nobody_sterile=nobody & (life_worlds < 1), nobody_inhabited=nobody & (life_worlds >= 1),
            stuck_barrier=c_stuck & barrier, stuck_many=c_stuck & ~barrier & (N_ever >= 100), stuck_few=c_stuck & ~barrier & (N_ever < 100),
            stuck_nobarrier_lt1=c_stuck & ~barrier & (N_ever * p_exp < 1),
            unnoticed_cc=c_unnot & cc, unnoticed_mixed=c_unnot & mixed,
            alone_barrier=alone & barrier, nobody_barrier=nobody & barrier, alone_settled=alone & ~nobody,
            fast_alone=fast & alone, fast_nobody=fast & nobody, fast_stuck=fast & c_stuck,
        )
        return dict(logPw=logPw, W=W, f_now=f_now, nidx=nidx, Wr=Wr, N_ever=N_ever, N_alltime=N_alltime, mu_loc=mu_loc, x=x, mu=mu,
                    S_claim=S_claim, L_local=L_claim * L_low, L_survey=np.exp(-np.minimum(N_eff * (1 - S_claim), 700.0)),
                    N_now=N_now, L_seti=np.exp(-np.minimum(N_now * c_det, 700.0)), cls=cls, alone=alone, c_unnot=c_unnot, sub=sub,
                    log_pp=np.log10(np.maximum(f_site, 1e-300)) + logPw[0] + log_ext)     # P(technological civilisation | rocky habitable-zone planet of a Sun-like star)

    d = derive(D)
    dp = derive(D_prior) if cond else d          # prior statistics never use the Earth-conditioned delay
    L_time = timeline(lam, h, m_hard, gated, D, w_D, x_now, tau, use_early, P["use_late"])
    if early_bf:
        L_time = L_time * np.where(fast, float(early_bf), 1.0)
    w_site = L_time * (d["L_local"] if P["use_local"] else one) * (d["L_survey"] if P["use_surveys"] else one) * (d["L_seti"] if P["use_seti"] else one)
    W, Wr, nidx, N_ever, N_alltime, mu, x = d["W"], d["Wr"], d["nidx"], d["N_ever"], d["N_alltime"], d["mu"], d["x"]

    # ---------------- hypothetical future observations ----------------
    p_lifeG = -np.expm1(-lam)

    def fex(ex, Q):
        """Likelihood factor of a hypothetical observation; Q carries its error-rate parameters."""
        if ex in ("second_genesis", "second_genesis_none"):     # an independent origin of life in the Solar System, community-accepted
            det = Q["sg_fp"] + (1 - Q["sg_fp"]) * -np.expm1(-Q["sg_frac"] * lam)
            return det if ex == "second_genesis" else 1 - det
        if ex in ("hwo_none", "hwo_some"):
            p_oxy = np.where(m >= 1, p_steps[:, 0], 1.0)        # oxygenic photosynthesis as the first step
            p_bio = f_site * p_lifeG * np.sqrt(f_pers) * (0.3 + 0.7 * p_oxy) * Q["hwo_q"]
            none = (1 - np.minimum(p_bio, 1.0)) ** 25 * (1 - Q["hwo_fp"])
            return none if ex == "hwo_none" else 1 - none
        if ex == "techno_null_1e6":
            return np.exp(-np.minimum(d["N_now"] * 1e6 / Nstars, 700))
        if ex == "techno_null_galaxy":
            return np.exp(-np.minimum(d["N_now"], 700))
        if ex == "artefact_null":
            low_present = np.where(mixed, -np.expm1(-mu * (1 - c)), np.where(cc, -np.expm1(-mu), 0.0))     # a low-footprint lineage is here
            return 1 - low_present * Q["art_d"]
        if ex == "exp_feasible":                                # no universal barrier (round 3: the old row required p_exp > 1e-3)
            return (~barrier).astype(float)
        if ex == "exp_common":
            return (p_exp > 1e-3).astype(float)
        if ex == "exp_infeasible":
            return barrier.astype(float)
        if ex == "mdwarf_alive":
            return M_alive.astype(float)
        if ex == "mdwarf_dead":
            return (~M_alive).astype(float)
        if ex == "planets_long_lived":
            return (~cap_short).astype(float)
        if ex == "planets_short_lived":
            return cap_short.astype(float)
        return one

    f_ex = fex(P.get("extra_obs"), P)
    # ---------------- reference classes: civilisations that arise without being pre-empted, plus pocket civilisations ----------------
    D_all = N_alltime.copy(); D_now = N_ever.copy()
    if P["use_local"]:
        E_inf = c * (kv * (Wr * TAB["kern"][nidx, lidx, -1]).sum(0) * P["use_cosmic"] + p_exp * N_alltime)
        idx = np.where(mixed & (E_inf > 1e-3))[0]
        for a in range(0, len(idx), 20000):
            ii = idx[a:a + 20000]
            Wc = Wr[:, ii][:, :, None]
            Ncum = (Wc * TAB["F"][nidx[:, ii], lidx[:, ii]]).sum(0)                       # (chunk, nt)
            tq = T_GRID[None, :] - T_fill[ii, None]
            j = np.clip(np.searchsorted(T_GRID, tq.ravel()).reshape(tq.shape), 1, len(T_GRID) - 1)
            wq_ = np.clip((tq - T_GRID[j - 1]) / (T_GRID[j] - T_GRID[j - 1]), 0.0, 1.0)
            Nlag = (1 - wq_) * np.take_along_axis(Ncum, j - 1, 1) + wq_ * np.take_along_axis(Ncum, j, 1)
            Nlag = np.where(tq <= T_GRID[0], 0.0, Nlag)                                   # the same fill delay as in mu_loc
            E = c[ii, None] * (kv[ii, None] * (Wc * TAB["kern"][nidx[:, ii], lidx[:, ii]]).sum(0) * P["use_cosmic"] + p_exp[ii, None] * Nlag)
            dN = np.diff(Ncum, axis=1); dE = np.clip(np.diff(E, axis=1), 0.0, 700.0)      # the cumulative exposure cannot decrease (rounding in the tables)
            fac = np.where(dE > 1e-8, -np.expm1(-dE) / np.where(dE > 1e-8, dE, 1.0), 1.0)  # product integration
            terms = dN * np.exp(-np.minimum(E[:, :-1], 700)) * fac
            Dt_all = terms.sum(1); Dt_now = terms[:, :I0].sum(1)
            D_all[ii] = Dt_all + eps[ii] * np.maximum(N_alltime[ii] - Dt_all, 0.0)
            D_now[ii] = Dt_now + eps[ii] * np.maximum(N_ever[ii] - Dt_now, 0.0)
    imp = np.exp(log_imp)
    share_all = W[0] / np.maximum(D_all, 1e-300)
    share_now = W[0] / np.maximum(D_now, 1e-300)
    conv = dict(site=one, ssa_all=share_all, ssa_now=share_now, ssa_obs=share_all / (1 + p_exp * P["obs_R"]), sia=W[0])
    nofx = {k: w_site * imp * f for k, f in conv.items()}        # today's posterior weights
    weights = {k: w * f_ex for k, w in nofx.items()}

    def wq(val, wn, qs=(0.05, 0.25, 0.5, 0.75, 0.95)):
        o = np.argsort(val); cw = np.cumsum(wn[o])
        return {str(q): float(val[o][min(np.searchsorted(cw, q), n - 1)]) for q in qs}

    def summarise(dd, wts):
        cls_, alone_, Ne, Nn, mu_, x_ = dd["cls"], dd["alone"], dd["N_ever"], dd["N_now"], dd["mu"], dd["x"]
        res = {}
        for name, w in wts.items():
            ws = w.sum(); wn = w / ws
            p_now = -np.expm1(-np.minimum(Nn, 700))
            res[name] = dict(
                posterior={k: float(wn[cls_ == i].sum()) for i, k in enumerate(CLASSES)},
                ess=float(1 / np.sum(wn ** 2)), evidence=float(ws / n),
                mean_like={k: float((w[cls_ == i] / imp[cls_ == i]).mean()) if (cls_ == i).any() else 0.0 for i, k in enumerate(CLASSES)},
                sub={k: float(wn[v_].sum()) for k, v_ in dd["sub"].items()},
                p_no_other_ever_in_galaxy=float(wn[alone_].sum()),
                p_no_other_ever_poisson=float((wn * np.exp(-np.minimum(Ne, 700))).sum()),
                p_other_now_in_galaxy=float((wn * p_now).sum()),
                p_anyone_now_in_galaxy=float((wn * (1 - (1 - p_now) * np.exp(-np.minimum(mu_, 700)))).sum()),     # home-grown and active, or settlers present
                p_fast=float(wn[fast].sum()), p_gated=float(wn[gated].sum()), mean_hard_steps=float((wn * m_hard).sum()),
                p_hazard=float(wn[hazard].sum()), p_M_alive=float(wn[M_alive].sum()), p_cap_short=float(wn[cap_short].sum()),
                p_life_common=float(wn[life_worlds >= 1e3].sum()),
                p_barrier=float(wn[barrier].sum()), p_pexp_lt_1e6=float(wn[p_exp < 1e-6].sum()),
                mean_p_exp=float((wn * p_exp).sum()), mean_p_exp_evident=float((wn * p_exp * c).sum()),
                p_expanders_hold_1pct=float(wn[x_ >= 0.01].sum()), p_claimers_hold_1pct=float(wn[x_ * c >= 0.01].sum()),
                p_fast_claimers_hold_1pct=float(wn[(x_ * c >= 0.01) & (v_ig >= 0.1)].sum()),
                p_nobody_in_reach=float(wn[Ne * P["N_reach"] < 1].sum()), p_no_expander_in_reach=float(wn[Ne * p_exp * P["N_reach"] < 1].sum()),
                p_nobody_in_reach_c=float(wn[Ne * P["N_reach"] * 1e3 < 1].sum()),                                   # the region reachable at light speed
                regime=dict(mixed=float(wn[mixed].sum()), low_footprint=float(wn[cc].sum()), stealth=float(wn[stealth].sum())),
                log10_N_ever_q=wq(np.log10(np.maximum(Ne, 1e-300)), wn),
                log10_lineages_unnoticed_q=wq(np.log10(np.maximum(mu_, 1e-300)), np.where(dd["c_unnot"], wn, 0) / max(wn[dd["c_unnot"]].sum(), 1e-300), (0.05, 0.5, 0.95)),
            )
        return res

    pn = imp / imp.sum()
    prior = {k: float(pn[dp["cls"] == i].sum()) for i, k in enumerate(CLASSES)}
    prior_stats = dict(p_no_other_ever_in_galaxy=float(pn[dp["alone"]].sum()), p_fast=float(pn[fast].sum()), p_gated=float(pn[gated].sum()),
                       mean_hard_steps=float((pn * m_hard).sum()), p_barrier=float(pn[barrier].sum()), p_pexp_lt_1e6=float(pn[p_exp < 1e-6].sum()),
                       p_life_common=float(pn[life_worlds >= 1e3].sum()), p_expanders_hold_1pct=float(pn[dp["x"] >= 0.01].sum()),
                       sub={k: float(pn[v_].sum()) for k, v_ in dp["sub"].items()},
                       log10_p_civ_per_planet_q=wq(dp["log_pp"], pn), log10_p_exp_q=wq(np.log10(np.maximum(p_exp, 1e-40)), pn),
                       log10_N_ever_q=wq(np.log10(np.maximum(dp["N_ever"], 1e-300)), pn))
    out = dict(prior=prior, schools=summarise(d, weights), prior_stats=prior_stats)

    def add_predictive(sch, wts, ex, Q):
        for name in wts:                          # probability of the hypothetical observation under today's posterior
            pr = float(wts[name].sum() / nofx[name].sum())
            sch[name]["predictive"] = pr
            fp = Q["sg_fp"] if ex == "second_genesis" else Q["hwo_fp"] if ex == "hwo_some" else None
            if fp is not None:                    # the part that is not the assumed false-alarm rate
                sch[name]["predictive_true"] = (pr - fp) / (1 - fp)

    add_predictive(out["schools"], weights, P.get("extra_obs"), P)
    if P.get("futures"):
        out["futures"] = {}
        for key, ex, ov in P["futures"]:
            Q = dict(P, **ov)
            wts = {k: w * fex(ex, Q) for k, w in nofx.items()}
            sch = summarise(d, wts)
            add_predictive(sch, wts, ex, Q)
            out["futures"][key] = dict(prior=prior, schools=sch, prior_stats=prior_stats)
    if P.get("ablate"):                           # every subset of the four observation terms, Earth-as-given column, same samples
        out["ablation"] = {}
        Lt = {}
        for e in (False, True):
            for l in (False, True):
                if e and l and cond:
                    Lt[(e, l)] = (d, L_time)
                else:
                    Lt[(e, l)] = (dp, timeline(lam, h, m_hard, gated, D_prior, one, x_now, tau, e, l))
        import itertools
        for k in range(len(TERMS) + 1):
            for on in itertools.combinations(TERMS, k):
                dd, lt = Lt[("early" in on, "late" in on)]
                w = lt * (dd["L_local"] if "local" in on else one) * (dd["L_survey"] * dd["L_seti"] if "surveys" in on else one) * imp
                out["ablation"][subset_key(on)] = dict(prior=prior, schools=summarise(dd, {"site": w}), prior_stats=prior_stats)
    if keep:
        out["_s"] = dict(cls=d["cls"], weights=weights, imp=imp, N_ever=N_ever, N_now=d["N_now"], N_alltime=N_alltime, p_exp=p_exp, barrier=barrier, x=x, mu=mu,
                         mu_loc=d["mu_loc"], v=v, v_ig=v_ig, lam=lam, m=m, m_hard=m_hard, fast=fast, gated=gated, hazard=hazard, W=W, c=c, eps=eps, p_here=p_here,
                         regime=np.where(mixed, 0, np.where(cc, 1, 2)), L=dict(time=L_time, local=d["L_local"], survey=d["L_survey"], seti=d["L_seti"]),
                         M_alive=M_alive, cap_short=cap_short, G=G, life_worlds=life_worlds, log_pp=d["log_pp"], sub=d["sub"], f_now=d["f_now"], nidx=nidx,
                         D_all=D_all, D_now=D_now, LK=LK, LM=LM, Nstars=Nstars, prior_cls=dp["cls"], prior_N_ever=dp["N_ever"])
    return out


def _avg(vals, wts):
    if isinstance(vals[0], dict):
        return {k: _avg([v[k] for v in vals], wts) for k in vals[0]}
    return float(np.dot(vals, wts) / np.sum(wts))


def pool(runs):
    """Combine independent runs (different seeds) into one result: posteriors weighted by each run's total weight,
    prior statistics by a plain mean."""
    eq = np.ones(len(runs))
    out = dict(prior=_avg([r["prior"] for r in runs], eq), schools={}, prior_stats=_avg([r["prior_stats"] for r in runs], eq))
    for s in runs[0]["schools"]:
        wts = np.array([r["schools"][s]["evidence"] for r in runs])
        q = {k: _avg([r["schools"][s][k] for r in runs], wts) for k in runs[0]["schools"][s] if k not in ("ess", "evidence")}
        q["evidence"] = float(wts.mean())
        q["ess"] = float(wts.sum() ** 2 / np.sum([w ** 2 / r["schools"][s]["ess"] for w, r in zip(wts, runs)]))
        out["schools"][s] = q
    for key in ("futures", "ablation"):
        if key in runs[0]:
            out[key] = {k: pool([r[key][k] for r in runs]) for k in runs[0][key]}
    return out


if __name__ == "__main__":
    import time
    t = time.time()
    r = run(vary(BASE, n=200_000))
    print("seconds", round(time.time() - t, 1))
    print("class".ljust(16) + "prior   " + "  ".join(f"{s:8s}" for s in SCHOOLS))
    for k in CLASSES:
        print(k.ljust(16) + f"{r['prior'][k]:.3f}   " + "  ".join(f"{r['schools'][s]['posterior'][k]:8.3f}" for s in SCHOOLS))
    for s in SCHOOLS:
        q = r["schools"][s]
        print(s, {k: (round(v, 4) if isinstance(v, float) else v) for k, v in q.items() if k not in ("posterior", "mean_like", "sub", "log10_N_ever_q")})
        print("   sub", {k: round(v, 3) for k, v in q["sub"].items()})
    print("mean likelihood by class (site):", {k: round(v, 3) for k, v in r["schools"]["site"]["mean_like"].items()})
    print("prior stats:", r["prior_stats"])