Layer 1 complete: E3-E6 + E2 analysis add-ons

Finishes the Layer 1 analytical core. All six experiments run with honest,
publication-quality figures; 71 tests green.

- E3 region-matched grounding: `grounding.exercised` knob + per-region tail
  survival. Matched holds the exercised region's tail (0.49) where uniform
  spreads thin and lets it collapse (0.07).
- E4 multi-teacher recombination: `run_coverage` runner. Union coverage matches
  U(K_T,rho,q) exactly. Finding: mean-mixture distillation shows NO surviving
  benefit (a conservation law — 1/K_T dilution cancels the union gain); a
  union-preserving max-merge (M2N2-style) does. E4 reports both operators.
- E5 QD vs greedy: greedy drives fixation (H~0.01); QD holds H at 0.48-0.88,
  rising with the novelty exponent.
- E6 re-mint gate: `arm` multi-override sweep. Re-minting a collapsed lineage
  locks in divergence of KL-to-original; gating on diversity prevents it.
- E2 analysis add-ons (from the companion work order, numbers verified): new
  analysis.py (reduce_to_stationary, critical_grounding with bootstrap CI ->
  g*=0.048, 95% CI [0.047,0.050]); tail_band_metrics + per-band logging; the
  E2 figure rebuilt as a 2x2 (defined g*+CI, g=0 flagged as a finite-time
  artifact, tail item-vs-mass, per-rarity-band panel). Uses truth-mass-weighted
  tail coverage rather than the raw (martingale) tail_mass.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

src/knowledge/analysis.py

@@ -0,0 +1,103 @@
+"""Post-hoc analysis of experiment results (pure NumPy/pandas, seeded, deterministic).
+
+Turns a per-generation results frame into stationary summaries and an operational,
+CI-bearing definition of the critical grounding fraction g*. Nothing here touches the
+dynamics; it is analysis only, so figures remain a pure function of results.parquet.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+
+def reduce_to_stationary(df, value_col="heterozygosity", sweep_col="g",
+                         replicate_col="replicate", gen_col="generation", last_frac=0.33):
+    """Per-generation frame -> one stationary value per (sweep, replicate).
+
+    Averages ``value_col`` over the final ``last_frac`` of generations. Works for any
+    logged metric (heterozygosity, tail_mass, ...).
+
+    Args:
+        df (pd.DataFrame): Long-form results with sweep, replicate, generation, value cols.
+        value_col (str): Metric to reduce.
+        sweep_col (str): Swept parameter column.
+        replicate_col (str): Replicate id column.
+        gen_col (str): Generation column.
+        last_frac (float): Fraction of the tail of the trajectory to average.
+
+    Returns:
+        pd.DataFrame: One row per (sweep_col, replicate_col) with the stationary value.
+    """
+    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).
+
+    Returns (g_star, status) with status in {'ok', 'below_grid', 'above_grid'}, by linear
+    interpolation between grid points.
+    """
+    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* with a percentile-bootstrap CI.
+
+    g* is the grounding fraction at which stationary heterozygosity first reaches
+    ``frac * H_star``. The point estimate uses the replicate means; the CI resamples
+    replicates within each g.
+
+    Args:
+        stationary_df (pd.DataFrame): One row per (sweep, replicate) with the stationary
+            value (e.g. the output of :func:`reduce_to_stationary`).
+        H_star (float): Heterozygosity of the truth (``metrics.heterozygosity(p_star)``).
+        frac (float): Fraction of H_star that defines "saturated".
+        sweep_col (str): Swept-parameter column.
+        value_col (str): Stationary value column.
+        n_boot (int): Bootstrap resamples.
+        ci (tuple): Percentile CI bounds.
+        seed (int): Bootstrap seed (deterministic).
+
+    Returns:
+        dict: g_star, ci_low, ci_high, status, target_H, frac, n_boot. ``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}

src/knowledge/config.py

@@ -37,6 +37,8 @@ class TeachersCfg:
 class GroundingCfg:
     m: int = 0
     policy: str = "uniform"      # {proportional, uniform, matched}
+    # Region indices exercised this passage (matched policy only; None = all regions).
+    exercised: Optional[list] = None
 
 
 @dataclass(frozen=True)

src/knowledge/experiment.py

@@ -28,6 +28,8 @@ import yaml
 
 from .lineage import run_lineage
 from .seeding import spawn_seeds
+from .teachers import make_correlated_teachers
+from .truth import make_true_distribution
 
 # Keys that make up a single-lineage configuration (everything else is experiment-level).
 _LINEAGE_KEYS = ("truth", "dynamics", "generations", "metrics")
@@ -55,6 +57,12 @@ def _apply_param(lineage_cfg: dict, param: str, value: Any) -> dict:
         m = 0 if g <= 0.0 else int(round(n * g / (1.0 - g)))
         lineage_cfg["dynamics"].setdefault("grounding", {})["m"] = m
         return {"g": g, "m": m}
+    if param == "arm":
+        # A named arm bundling several overrides applied together (e.g. E6 varies
+        # grounding + remint settings jointly). value = {name, set: {dotted.path: v}}.
+        for path, v in value.get("set", {}).items():
+            _set_by_path(lineage_cfg, path, v)
+        return {"arm": value["name"]}
     _set_by_path(lineage_cfg, param, value)
     return {param.split(".")[-1]: value}
 
@@ -119,6 +127,75 @@ def run_experiment(cfg: dict) -> pd.DataFrame:
     return out
 
 
+def run_coverage(cfg: dict) -> pd.DataFrame:
+    """E4 runner: multi-teacher recombination coverage (blueprint 2.5-E4 / 2.7.1).
+
+    A single distillation step, not a lineage. For each (K_T, rho) grid point and
+    replicate: build K_T correlated teachers (marginal retention q, pairwise
+    correlation rho), then measure two coverages of the tail:
+
+    * ``union_coverage`` — fraction of tail items retained by >=1 teacher (the
+      construction-level union U(K_T, rho, q); must match the closed form). This is the
+      recombination *supply*.
+    * ``surviving_mean`` / ``surviving_max`` — fraction of tail items that survive the
+      pupil's size-n resampling (+ optional grounding m) under two recombination
+      operators: ``mean`` (blueprint mean-mixture distillation) and ``max`` (union-
+      preserving model-merge, à la M2N2). Under ``mean`` the union gain is diluted by
+      1/K_T and (in the rare-tail linear regime) is exactly cancelled — expected pupil
+      tail mass is conserved at q·(tail mass) regardless of K_T, so surviving is flat.
+      Under ``max`` each item keeps its strongest teacher, so surviving rises with K_T
+      and with decorrelation. The gap ``union - surviving`` is the tail recombination
+      supplied but drift/dilution re-erased.
+
+    Matched budget: the pupil draws n samples total from the combined teachers
+    (equivalently n/K_T each), so more teachers != more data.
+    """
+    truth, cov = cfg["truth"], cfg["coverage"]
+    td = make_true_distribution(
+        truth["K"], truth["R"], truth["tail"], truth["tail_frac"], truth["zipf_s"], 0,
+        tail_threshold=truth["tail_threshold"],
+    )
+    tail_idx = np.flatnonzero(td.tail_mask)
+    n, q = int(cov["n"]), float(cov["q"])
+    retain_thresh = 1e-8  # dropped tails sit at ~tail_floor (1e-9); retained at ~p*_j
+
+    sweeps = cfg["sweep"]
+    if isinstance(sweeps, dict):
+        sweeps = [sweeps]
+    params = [s["param"] for s in sweeps]
+    value_lists = [list(s["values"]) for s in sweeps]
+    seeds = spawn_seeds(int(cfg["seed"]), int(cfg["n_replicates"]))
+
+    rows: list[dict] = []
+    for combo in itertools.product(*value_lists):
+        d = dict(zip(params, combo))
+        K_T, rho = int(d["K_T"]), float(d["rho"])
+        g = float(d.get("g", cov.get("g", 0.0)))          # g may be swept or fixed
+        m = 0 if g <= 0.0 else int(round(n * g / (1.0 - g)))
+        for rep, ss in enumerate(seeds):
+            child = int(ss.generate_state(1)[0])
+            teachers = np.asarray(make_correlated_teachers(
+                td.p_star, td.tail_mask, K_T, rho, q, seed=child))
+            retained = teachers[:, tail_idx] > retain_thresh          # (K_T, T)
+            union = float(np.mean(retained.any(axis=0)))
+            surviving = {}
+            for offset, combine in ((1, teachers.mean), (2, teachers.max)):
+                p = combine(axis=0)
+                p = p / p.sum()
+                rng = np.random.default_rng(child + offset)
+                counts = rng.multinomial(n, p)
+                if m > 0:
+                    counts = counts + rng.multinomial(m, td.p_star)
+                surviving[offset] = float(np.mean(counts[tail_idx] > 0))
+            rows.append({
+                "experiment": cfg["experiment"], "K_T": K_T, "rho": rho,
+                "replicate": rep, "union_coverage": union,
+                "surviving_mean": surviving[1], "surviving_max": surviving[2],
+                "g": g, "q": q,
+            })
+    return pd.DataFrame(rows)
+
+
 def _git_commit() -> str | None:
     try:
         return subprocess.check_output(
@@ -148,12 +225,13 @@ def save_artifacts(cfg: dict, df: pd.DataFrame, out_dir: Path) -> None:
         "experiment": cfg["experiment"],
         "seed": cfg["seed"],
         "n_replicates": cfg["n_replicates"],
-        "grid": [
-            {"label": label, "lineage_cfg": lineage_cfg}
-            for label, lineage_cfg in expand_sweeps(cfg)
-        ],
         "source_config": cfg,
     }
+    if cfg.get("kind", "lineage") == "lineage":
+        resolved["grid"] = [
+            {"label": label, "lineage_cfg": lineage_cfg}
+            for label, lineage_cfg in expand_sweeps(cfg)
+        ]
     (out_dir / "resolved_config.yaml").write_text(yaml.safe_dump(resolved, sort_keys=False))
 
     manifest = {
@@ -175,7 +253,7 @@ def run_and_save(config_path: str | Path) -> Path:
     config_path = Path(config_path)
     cfg = yaml.safe_load(config_path.read_text())
     out_dir = Path(cfg.get("output", {}).get("dir", f"results/{cfg['experiment']}"))
-    df = run_experiment(cfg)
+    df = run_coverage(cfg) if cfg.get("kind") == "coverage" else run_experiment(cfg)
     save_artifacts(cfg, df, out_dir)
     return out_dir
 

src/knowledge/lineage.py

@@ -14,10 +14,13 @@ import numpy as np
 import pandas as pd
 
 from .config import LineageCfg
-from .metrics import forward_kl, heterozygosity, per_region, support_size, tail_mass
+from .metrics import (forward_kl, heterozygosity, per_region, support_size,
+                      tail_band_metrics, tail_mass)
 from .step import StepCtx, allocate_m, generation_step
 from .truth import make_true_distribution, uniform_init
 
+N_BANDS = 4  # rarity bands for per-band tail-survival logging (blueprint pred. 4)
+
 
 def run_lineage(cfg: Mapping[str, Any] | LineageCfg, seed: int) -> pd.DataFrame:
     """Run one lineage and return per-generation metrics.
@@ -52,8 +55,10 @@ def run_lineage(cfg: Mapping[str, Any] | LineageCfg, seed: int) -> pd.DataFrame:
         raise ValueError(f"unknown init {cfg.truth.init!r} (expected uniform|truth)")
 
     p_star_eff = p_star_orig.copy()  # grounding reference; may be re-minted (E6)
+    exercised = cfg.dynamics.grounding.exercised
+    exercised = np.asarray(exercised) if exercised is not None else None
     m_vector = allocate_m(cfg.dynamics.grounding.m, cfg.truth.R,
-                          cfg.dynamics.grounding.policy)
+                          cfg.dynamics.grounding.policy, exercised)
     step_ctx = StepCtx(
         n=cfg.dynamics.n,
         m_vector=m_vector,
@@ -83,12 +88,29 @@ def run_lineage(cfg: Mapping[str, Any] | LineageCfg, seed: int) -> pd.DataFrame:
             "head_support": int(np.sum(p[head_mask] > eps)),
             "tail_frac_alive": (float(np.mean(p[tail_mask] > eps)) if n_tail else 0.0),
             "head_frac_alive": (float(np.mean(p[head_mask] > eps)) if n_head else 0.0),
+            # Truth-mass-weighted tail coverage: share of the tail's TRUE mass carried by
+            # still-alive items. Bounded, monotone in grounding, and the honest "mass
+            # rescued" companion to tail_frac_alive (raw tail_mass is a drift martingale).
+            "tail_truth_mass_alive": (
+                float(p_star_orig[tail_mask][p[tail_mask] > eps].sum()
+                      / p_star_orig[tail_mask].sum()) if n_tail else 0.0),
         }
+        if n_tail >= N_BANDS:  # per-rarity-band survival (band 0 = rarest); see metrics
+            fa, _ = tail_band_metrics(p, p_star_orig, tail_mask, n_bands=N_BANDS,
+                                      alive_eps=eps)
+            for b in range(N_BANDS):
+                row[f"band{b}_alive"] = fa[b]
         if cfg.truth.R > 1:  # per-region columns only when there is >1 region
             for r, v in per_region(heterozygosity, p, regions).items():
                 row[f"H_region_{r}"] = v
             for r, v in per_region(tail_mass, p, regions, tail_mask).items():
                 row[f"tail_region_{r}"] = v
+            # per-region tail-item survival (E3's honest metric: how much of each
+            # region's rare tail is kept alive, not just its martingale mass)
+            for r in range(cfg.truth.R):
+                region_tail = (regions == r) & tail_mask
+                row[f"tailalive_region_{r}"] = (
+                    float(np.mean(p[region_tail] > eps)) if region_tail.any() else 0.0)
         rows.append(row)
 
     record(0, p)

src/knowledge/metrics.py

@@ -79,6 +79,42 @@ def support_size(p: np.ndarray, eps: float) -> int:
     return int(np.sum(p > eps))
 
 
+def tail_band_metrics(p: np.ndarray, p_star: np.ndarray, tail_mask: np.ndarray,
+                      n_bands: int = 4, alive_eps: float = 1e-9):
+    """Stratify the tail into ``n_bands`` equal-count rarity bands (band 0 = rarest).
+
+    Makes the per-item survival threshold ``m·p*_i ≳ 1`` (blueprint prediction 4) visible
+    band-wise: deeper (rarer) bands sit strictly below shallower ones and the gap narrows
+    as grounding rises. Both returned arrays are bounded in [0, 1]; single-run values are
+    noisy, so average over replicates before plotting.
+
+    Args:
+        p (np.ndarray): Current distribution.
+        p_star (np.ndarray): True distribution.
+        tail_mask (np.ndarray): Boolean tail mask.
+        n_bands (int): Number of equal-count rarity bands.
+        alive_eps (float): An item is "alive" if ``p_i > alive_eps``.
+
+    Returns:
+        tuple[np.ndarray, np.ndarray]: ``frac_alive[b]`` (fraction of band-b items alive)
+        and ``truth_mass_alive[b]`` (share of band-b's TRUE mass carried by alive items),
+        each length ``n_bands``.
+    """
+    p = np.asarray(p, dtype=float)
+    p_star = np.asarray(p_star, dtype=float)
+    idx = np.where(np.asarray(tail_mask, dtype=bool))[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
+
+
 def per_region(func, p: np.ndarray, regions: np.ndarray, *args) -> dict[int, float]:
     """Apply a metric independently to each region's sub-vector.