# 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/knowledge/analysis.py`

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

```python
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/knowledge/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).

```python
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_E2.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`)

```python
import numpy as np, pandas as pd, pytest
from knowledge.analysis import reduce_to_stationary, critical_grounding
from knowledge.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.