MachineSex/tasks/workorder-E2-analysis-addons.md

Work order: E2 analysis add-ons

Increment to the Lamarckian-Society Layer-1 build. Assumes the blueprint conventions (seeded RNG, config-driven, results in tidy DataFrames, figures are pure functions of results.parquet). All code below was written and tested against NumPy 2.x / pandas 3.x before hand-off; the "verified numbers" section gives the expected outputs so you can confirm your wiring reproduces them.

Context

E1 and E2 are validated: E1 reproduces the analytic geometric decay H0(1-1/n)^t (tail-first, support collapse, KL divergence); E2's stationary simulation points lie on the exact closed form H_eq = H* * m(2n+m-1)/(n+2nm+m^2). Three refinements make the E2 figure publication-honest. None change the dynamics; they are analysis + plotting only.

1. Operational g* with a bootstrap CI. The middle E2 panel is a smooth saturating curve, not a sharp transition, so "critical g*" must be defined, not asserted: g* = the grounding fraction at which stationary H first reaches frac×H* (default frac=0.95). Report it with a percentile-bootstrap CI over replicates.
2. Tail-mass companion curve. The right E2 panel plots fraction of tail items alive, which saturates low (~6%) because the deep Zipf tail is unrescuable at any feasible grounding. Add the frequency-weighted tail-mass-retained curve alongside it: mass is carried by the shallow tail and is rescued, so it looks healthy where item-count does not. Both are correct; showing both is the honest picture. (tail_mass is already a logged metric — this is a second series, no new storage.)
3. g=0 point is a finite-time artifact. In the middle panel it sits at H≈0.10, above the exact H_eq(0)=0, because at g=0 there is no stationary state (still collapsing at the last generation). Either exclude g=0 from the stationary fit or annotate that its true value is 0; don't let it read as a data–theory miss.

Optional but recommended (turns point 2 into a mechanism panel): a per-rarity-band tail metric showing the deep band stays dead while the shallow band recovers — the m·p*_i ≳ 1 threshold acting band-wise, and the direct motivation for E4/E6.

---

1. New module: src/inheritance/analysis.py

Post-hoc analysis of E2 results. Pure NumPy/pandas, seeded, deterministic.

import numpy as np
import pandas as pd


def reduce_to_stationary(df, value_col="heterozygosity", sweep_col="g",
                         replicate_col="seed", gen_col="generation", last_frac=0.33):
    """Per-generation frame -> one stationary value per (sweep, replicate), averaging
    `value_col` over the final `last_frac` of generations. Use for any logged metric
    (heterozygosity, tail_mass, ...)."""
    rows = []
    for (gval, rep), sub in df.groupby([sweep_col, replicate_col]):
        v = sub.sort_values(gen_col)[value_col].to_numpy()
        k = max(1, int(round(last_frac * v.size)))
        rows.append({sweep_col: gval, replicate_col: rep, value_col: v[-k:].mean()})
    return pd.DataFrame(rows)


def _interp_crossing(g, H, target):
    """First upward crossing of `target` by the (monotone-ish) curve H(g), by linear
    interpolation between grid points. Returns (g_star, status) with status in
    {'ok', 'below_grid', 'above_grid'}."""
    g = np.asarray(g, float); H = np.asarray(H, float)
    o = np.argsort(g); g, H = g[o], H[o]
    if H[0] >= target:
        return g[0], "below_grid"        # already above at smallest g swept
    idx = np.where(H >= target)[0]
    if idx.size == 0:
        return g[-1], "above_grid"       # never reaches target within swept range
    i = idx[0]
    g0, g1, H0, H1 = g[i - 1], g[i], H[i - 1], H[i]
    if H1 == H0:
        return g1, "ok"
    return g0 + (target - H0) * (g1 - g0) / (H1 - H0), "ok"


def critical_grounding(stationary_df, H_star, frac=0.95, sweep_col="g",
                       value_col="heterozygosity", n_boot=2000,
                       ci=(2.5, 97.5), seed=0):
    """Operational critical grounding fraction g*: the g at which stationary
    heterozygosity first reaches `frac` * `H_star`, with a percentile-bootstrap CI
    over replicates.

    `stationary_df`: one row per (sweep_col, replicate_col) with the stationary value
    (e.g. the output of reduce_to_stationary). `H_star`: heterozygosity of the truth,
    = metrics.heterozygosity(p_star). Returns a dict with g_star (point estimate on the
    replicate means), ci_low, ci_high, status, target_H, frac, n_boot.

    Note: 'status' flags right/left censoring. If the sweep does not bracket the target,
    widen the g grid rather than trusting a censored g*."""
    target = frac * H_star
    gs = np.sort(stationary_df[sweep_col].unique())
    by_g = {gv: stationary_df.loc[stationary_df[sweep_col] == gv, value_col].to_numpy()
            for gv in gs}
    mean_H = np.array([by_g[gv].mean() for gv in gs])
    g_star, status = _interp_crossing(gs, mean_H, target)

    rng = np.random.default_rng(seed)
    boots = np.empty(n_boot)
    for b in range(n_boot):
        Hb = np.array([rng.choice(by_g[gv], by_g[gv].size, replace=True).mean()
                       for gv in gs])
        boots[b], _ = _interp_crossing(gs, Hb, target)
    lo, hi = np.percentile(boots, ci)
    return {"g_star": float(g_star), "ci_low": float(lo), "ci_high": float(hi),
            "status": status, "target_H": float(target), "frac": frac,
            "n_boot": n_boot}

---

2. New online metric in src/inheritance/metrics.py (optional band panel)

Compute this each generation from the current p and log the per-band arrays exactly like the existing per-region metrics (e.g. columns tail_frac_alive_band{b} and tail_mass_alive_band{b}, or long form). Both returned quantities are bounded in [0,1] — do not use a raw mass ratio (tiny deep-band denominators make it explode).

def tail_band_metrics(p, p_star, tail_mask, n_bands=4, alive_eps=1e-9):
    """Stratify the tail into `n_bands` equal-count rarity bands (band 0 = rarest /
    deepest). Return two BOUNDED [0,1] arrays of length n_bands:
      frac_alive[b]        fraction of band-b items with p > alive_eps
      truth_mass_alive[b]  share of band-b's TRUE mass (sum p_star) carried by
                           still-alive items
    Averaged over replicates, deeper bands sit strictly below shallower ones and the
    gap narrows as grounding rises -- the per-item m*p*_i >~ 1 survival threshold made
    visible (blueprint prediction 4). Single-run values are noisy; always average over
    replicates before plotting."""
    import numpy as np
    idx = np.where(tail_mask)[0]
    order = idx[np.argsort(p_star[idx])]           # rarest first
    bands = np.array_split(order, n_bands)
    frac_alive = np.empty(n_bands)
    truth_mass_alive = np.empty(n_bands)
    for b, items in enumerate(bands):
        alive = p[items] > alive_eps
        frac_alive[b] = alive.mean()
        ps = p_star[items]
        truth_mass_alive[b] = ps[alive].sum() / ps.sum() if ps.sum() > 0 else np.nan
    return frac_alive, truth_mass_alive

---

3. Figure updates: figures/plot_fig2_grounding_sweep.py

• Middle panel: call critical_grounding(reduce_to_stationary(df), H_star, frac=0.95) and draw a vertical line/marker at g_star with a shaded CI band; annotate g* ≈ {g_star:.3f} (95% CI [...]). Also report frac=0.90 in the caption so the "3–5%" range is explicit. Soften the title from "critical grounding ratio" to e.g. "grounding saturates by g* ≈ 0.05 (95% of H*)".
• Right panel: add stationary tail-mass-retained vs g as a second series (reduce_to_stationary(df, value_col="tail_mass")), on a twin axis if scales differ, so mass-healthy vs items-poor is visible in one panel.
• g=0: drop it from the stationary fit or mark it hollow with a note "pre-convergence; true H_eq(0)=0".
• Optional band panel: if tail_band_metrics is logged, add a small-multiples or grouped-bar panel of frac_alive per band across the g sweep (replicate-averaged).

---

4. Tests to add (tests/test_analysis.py)

import numpy as np, pandas as pd, pytest
from inheritance.analysis import reduce_to_stationary, critical_grounding
from inheritance.metrics import tail_band_metrics


def _H_eq(n, m, Hs): return Hs * m * (2 * n + m - 1) / (n + 2 * n * m + m * m)


def _synthetic_E2(n=200, Hs=0.95, reps=100, seed=1):
    gs = [0.0, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.4]
    rng = np.random.default_rng(seed); rows = []
    for g in gs:
        m = 0 if g == 0 else int(round(g * n / (1 - g)))
        Htrue = 0.10 if g == 0 else _H_eq(n, m, Hs)   # g=0: finite-time artifact
        for s in range(reps):
            rows.append({"g": g, "seed": s, "heterozygosity": Htrue + rng.normal(0, 0.004)})
    return pd.DataFrame(rows)


def test_critical_grounding_matches_known_crossing():
    stat = _synthetic_E2()
    r = critical_grounding(stat, H_star=0.95, frac=0.95, seed=7)
    assert r["status"] == "ok"
    assert 0.03 < r["g_star"] < 0.07                   # ~0.047 for frac=0.95
    assert r["ci_low"] <= r["g_star"] <= r["ci_high"]
    r90 = critical_grounding(stat, H_star=0.95, frac=0.90, seed=7)
    assert r90["g_star"] < r["g_star"]                 # lower bar -> smaller g*


def test_reduce_to_stationary_recovers_plateau():
    # constant plateau + noise -> mean ~ plateau
    rng = np.random.default_rng(0); rows = []
    for g, plateau in [(0.02, 0.843), (0.05, 0.908)]:
        for s in range(20):
            for t in range(300):
                v = (0.95 if t < 50 else plateau) + rng.normal(0, 0.003)
                rows.append({"g": g, "seed": s, "generation": t, "heterozygosity": v})
    red = reduce_to_stationary(pd.DataFrame(rows), last_frac=0.33)
    means = red.groupby("g")["heterozygosity"].mean()
    assert means[0.02] == pytest.approx(0.843, abs=0.01)
    assert means[0.05] == pytest.approx(0.908, abs=0.01)


def test_tail_band_deep_below_shallow():
    # grounded dynamics on a Zipf tail: deepest band <= shallowest band, replicate-avg
    K = 1000; s = 1.1
    w = 1.0 / np.arange(1, K + 1) ** s; p_star = w / w.sum()
    tail_mask = p_star < 1e-3
    n, m = 200, 10
    FA = np.zeros(4)
    for r in range(20):
        rng = np.random.default_rng(1000 + r); p = p_star.copy()
        for _ in range(600):
            c = rng.multinomial(n, p) + rng.multinomial(m, p_star); p = c / c.sum()
        fa, _ = tail_band_metrics(p, p_star, tail_mask, n_bands=4)
        FA += fa
    FA /= 20
    assert FA[0] <= FA[-1]                              # deepest no better than shallowest
    assert FA[-1] > FA[0]                               # and strictly worse on average

---

Verified numbers (expected outputs — confirm your wiring reproduces these)

On the synthetic E2 set (n=200, H*=0.95, exact H_eq + N(0,0.004) noise, 100 reps):

frac	g* (point)	95% CI (tight, low-noise synthetic)
0.90	~0.026	~[0.025, 0.026]
0.95	~0.047	~[0.047, 0.048]
0.99	~0.25	~[0.22, 0.27]

(Real E2 replicate spread will widen these CIs — that is expected and correct.)

reduce_to_stationary on a plateau frame recovers 0.843 (g=0.02) and 0.908 (g=0.05).

tail_band_metrics, grounded Zipf tail (K=1000, tail=p*<1e-3, 889 tail items), 40 reps, frac_alive per band [0=deepest .. 3=shallowest]:

g	band0	band1	band2	band3
0.0	0.000	0.000	0.000	0.000
0.01	0.002	0.002	0.004	0.010
0.05	0.004	0.009	0.012	0.033
0.20	0.015	0.022	0.032	0.075

Monotone (deeper = worse) at every g>0; all bands rise with g; deep band lags throughout.