neural: grounding refinement + all five Layer-1.5 figures

Grounding refinement (18 reps): forward-KL is the operative neural
collapse metric, not H or tail-survival. The RNN's smoothing keeps
spurious tail modes alive, so tail_truth_mass_alive is flat/non-monotone
in g and H stays ~0.8 of H*; only forward-KL falls monotonically (dry
2.08 -> g=0.2: 0.75, paired t up to 3.3). The sharp g* << 1 is an
exact-operator feature carried by the histogram bridge (0.047); the
trained RNN confirms the SIGN and softens the sharpness (half the KL gap
closes by g~0.04, but full recovery needs g~0.19). Blueprint 3.5's
directional claim holds; the pre-registered 95%-of-H*/tail falsifier is
not met because those are the wrong metrics for a smoothing model.

Robustness: a fully-degenerate RNN can emit only invalid codewords, so
measure_distribution now returns a terminal-collapse sentinel (fixation
on the dominant mode) instead of crashing a long sweep. Edge test added
(94 tests green).

Figures: plot_{bridge,collapse,grounding,architectures,recombination}.py,
each a pure function of its committed bundle, wired into `make figures`
(glob plot_*.py minus plot_E[1-6]/_*). bridge sits on the exact H_eq
curve (g*=0.047); recombination shows max-merge rising while mean-distill
stays flat; architectures shows the collapse/rescue signs across
histogram/GRU/MLP.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Giorgio Gilestro 2026-07-05 08:14:19 +01:00
parent d22dd9d535
commit b8da418034
23 changed files with 680 additions and 26 deletions

View file

@ -0,0 +1,82 @@
"""`architectures` figure — the Wright-Fisher collapse operator is architecture-general.
The same dry-collapse / grounding-rescue signature appears in three distinct inductive biases:
the exact histogram (multinomial), an autoregressive GRU, and a causal-masked MLP. If the
signs held only for the histogram, the effect would be an artefact of the exact operator;
seeing them in every trained architecture is the generality claim.
Three panels: (A) forward-KL trajectories per architecture, dry (solid) vs grounded (dashed);
(B) stationary forward-KL, dry vs grounded, grouped by architecture (all fall with grounding);
(C) tail-item survival, dry vs grounded, grouped by architecture (all rise). Reads only the
committed bundle.
Usage: python figures/plot_architectures.py [results/architectures]
"""
from __future__ import annotations
import sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
sys.path.insert(0, str(Path(__file__).parent))
from _figlib import load_bundle, mean_ci, savefig # noqa: E402
_ARCH_ORDER = ["histogram", "rnn", "mlp"]
_ARCH_LABEL = {"histogram": "histogram\n(exact)", "rnn": "GRU\n(autoregressive)",
"mlp": "MLP\n(causal-masked)"}
def main(results_dir: str = "results/architectures") -> None:
df, cfg = load_bundle(results_dir)
kinds = [k for k in _ARCH_ORDER if k in set(df["kind"].unique())]
g_dry, g_wet = min(df["g"].unique()), max(df["g"].unique())
last = int(cfg["generations"] * 0.6)
stat = df[df["generation"] >= last]
fig, axes = plt.subplots(1, 3, figsize=(16, 4.6))
arch_colors = dict(zip(kinds, plt.cm.tab10(np.arange(len(kinds)))))
# Panel A: forward-KL trajectories per architecture, dry (solid) vs grounded (dashed).
ax = axes[0]
for k in kinds:
for g, ls, alpha in [(g_dry, "-", 1.0), (g_wet, "--", 0.7)]:
s = df[(df["kind"] == k) & (df["g"] == g)].groupby("generation")["forward_kl"].mean()
ax.plot(s.index, s.values, ls, color=arch_colors[k], alpha=alpha, lw=1.8,
label=f"{k} (g={g:g})")
ax.set(xlabel="generation", ylabel=r"forward-KL $D(p^*\Vert\hat p)$",
title="Dry (solid) collapses; grounded (dashed) held —\nin every architecture")
ax.legend(frameon=False, fontsize=7, ncol=1)
# Panels B & C: grouped bars, dry vs grounded per architecture.
def grouped_bar(ax, metric, title, ylabel):
x = np.arange(len(kinds))
w = 0.36
for off, g, lab, col in [(-w / 2, g_dry, f"dry (g={g_dry:g})", "#d62728"),
(w / 2, g_wet, f"grounded (g={g_wet:g})", "#2ca02c")]:
means, errs = [], []
for k in kinds:
sub = stat[(stat["kind"] == k) & (stat["g"] == g)]
_, m, ci = mean_ci(sub.assign(_x=0), "_x", metric)
means.append(m[0]); errs.append(ci[0])
ax.bar(x + off, means, w, yerr=errs, capsize=3, label=lab, color=col, alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels([_ARCH_LABEL[k] for k in kinds], fontsize=8)
ax.set(ylabel=ylabel, title=title)
ax.legend(frameon=False, fontsize=8)
grouped_bar(axes[1], "forward_kl", "Stationary forward-KL falls with grounding",
r"stationary forward-KL")
grouped_bar(axes[2], "tail_frac_alive", "Tail-item survival rises with grounding",
"tail items alive")
fig.suptitle("architectures — dry collapse and grounding-rescue are architecture-general "
"(histogram, GRU, MLP)", y=1.02, fontsize=13)
fig.tight_layout()
savefig(fig, results_dir, "architectures")
if __name__ == "__main__":
main(*sys.argv[1:])

111
figures/plot_bridge.py Normal file
View file

@ -0,0 +1,111 @@
"""`bridge` figure — the histogram bridge reproduces Layer-1 E2 exactly (the HARD GATE).
The histogram model reduces Layer 1.5 to Layer 1 (MLE histogram + multinomial resampling),
so running it through the *neural* runner must reproduce E2's grounding phase boundary and its
exact `H_eq` closed form. Recovering g*0.047 here (vs Layer-1's 0.048) is what licenses every
later trained-model result to be read against the analytic core.
Four panels: (A) H trajectories (g=0 collapses, g>0 plateau); (B) the phase boundary
stationary H vs g on the exact `H_eq` curve, with the operational g* + bootstrap CI; (C) tail
survival vs g; (D) per-rarity-band survival (the m·p*_i1 threshold). Reads only the bundle.
Usage: python figures/plot_bridge.py [results/bridge]
"""
from __future__ import annotations
import sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
sys.path.insert(0, str(Path(__file__).parent))
from _figlib import load_bundle, mean_ci, savefig # noqa: E402
sys.path.insert(0, str(Path(__file__).parents[1] / "src"))
from knowledge.analysis import critical_grounding, reduce_to_stationary # noqa: E402
from knowledge.metrics import heterozygosity # noqa: E402
from neural.config import SyntheticCfg # noqa: E402
from neural.synthetic import make_mode_truth # noqa: E402
def main(results_dir: str = "results/bridge") -> None:
df, cfg = load_bundle(results_dir)
n = cfg["dynamics"]["n"]
syn = SyntheticCfg(**cfg["synthetic"])
H_star = heterozygosity(make_mode_truth(syn).p_star)
def H_eq(m):
m = np.asarray(m, dtype=float)
return np.where(m <= 0, 0.0, H_star * m * (2 * n + m - 1) / (n + 2 * n * m + m * m))
g_values = sorted(df["g"].unique())
last = int(cfg["generations"] * 0.8)
stat = df[df["generation"] >= last]
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
colors = plt.cm.viridis(np.linspace(0, 0.9, len(g_values)))
# Panel A: H trajectories, one line per g.
ax = axes[0, 0]
for g, c in zip(g_values, colors):
s = df[df["g"] == g].groupby("generation")["heterozygosity"].mean()
ax.plot(s.index, s.values, color=c, label=f"g={g:g}")
ax.axhline(H_star, ls=":", color="gray", lw=1)
ax.set(xlabel="generation", ylabel="heterozygosity $H$",
title="Trajectories: g=0 collapses, g>0 plateau")
ax.legend(frameon=False, fontsize=8, ncol=2)
# Panel B: stationary H vs g on the exact H_eq curve + operational g* with bootstrap CI.
ax = axes[0, 1]
st = reduce_to_stationary(stat, value_col="heterozygosity", last_frac=1.0)
gg, Hm, Hci = mean_ci(stat, "g", "heterozygosity")
m_of_g = stat.groupby("g")["m"].first().to_numpy()
nz = gg > 0
ax.errorbar(gg[nz], Hm[nz], yerr=Hci[nz], fmt="o", color="#1f77b4", capsize=3, zorder=3,
label="neural histogram runner")
ax.plot(gg[~nz], Hm[~nz], "o", mfc="white", mec="#1f77b4", zorder=3)
m_grid = np.linspace(0, m_of_g.max(), 400)
ax.plot(m_grid / (n + m_grid), H_eq(m_grid), "k--", zorder=2, label=r"exact $H_{eq}$ (Layer 1)")
ax.axhline(H_star, ls=":", color="gray", lw=1, label="$H^*$ (truth)")
r = critical_grounding(st, H_star=H_star, frac=0.95, seed=7)
ax.axvspan(r["ci_low"], r["ci_high"], color="#d62728", alpha=0.15)
ax.axvline(r["g_star"], color="#d62728", lw=1.2,
label=f"$g^*$={r['g_star']:.3f} (95% CI [{r['ci_low']:.3f},{r['ci_high']:.3f}])")
ax.set(xlabel="grounding fraction $g=m/(n+m)$", ylabel="stationary $H$",
title=r"Bridge reproduces the exact $H_{eq}$ and $g^\star$")
ax.legend(frameon=False, fontsize=8)
# Panel C: tail survival vs g (item-count and truth-mass).
ax = axes[1, 0]
ig, Im, Ici = mean_ci(stat, "g", "tail_frac_alive")
mg, Mm, Mci = mean_ci(stat, "g", "tail_truth_mass_alive")
ax.errorbar(ig, Im, yerr=Ici, fmt="s-", color="#d62728", capsize=3, label="tail items alive")
ax.errorbar(mg, Mm, yerr=Mci, fmt="o-", color="#9467bd", capsize=3, label="tail truth-mass alive")
ax.set(xlabel="grounding fraction $g$", ylabel="fraction of tail retained",
title="Tail survival rises with g (deep tail lags)")
ax.legend(frameon=False, fontsize=9)
# Panel D: per-rarity-band survival across g.
ax = axes[1, 1]
band_cols = sorted(c for c in df.columns if c.startswith("band") and c.endswith("_alive"))
band_colors = plt.cm.plasma(np.linspace(0.1, 0.85, len(band_cols)))
for col, c in zip(band_cols, band_colors):
s = stat.groupby("g")[col].mean()
depth = col.replace("band", "").replace("_alive", "")
lab = f"band {depth}" + (" (deepest)" if col == band_cols[-1] else
" (shallowest)" if depth == "0" else "")
ax.plot(s.index, s.values, "-o", color=c, ms=4, label=lab)
ax.set(xlabel="grounding fraction $g$", ylabel="fraction of band alive",
title=r"Per-rarity band: the $m\,p^*_i\gtrsim1$ threshold")
ax.legend(frameon=False, fontsize=8)
fig.suptitle("bridge — the histogram model reproduces Layer-1 E2 through the neural runner "
f"($g^*$={r['g_star']:.3f} vs Layer-1 0.048; HARD GATE passed)", y=1.0, fontsize=13)
fig.tight_layout()
savefig(fig, results_dir, "bridge")
if __name__ == "__main__":
main(*sys.argv[1:])

84
figures/plot_collapse.py Normal file
View file

@ -0,0 +1,84 @@
"""`collapse` figure — model collapse in REAL RNN weights, arrested by grounding (↔ E1/C1).
The existence proof: a GRU trained each generation on the previous generation's own samples
loses the rare tail and drifts from truth (forward-KL climbs), and even a little grounding
arrests it. Forward-KL is the operative neural collapse metric (the RNN's smoothing keeps
spurious tail support alive, so H barely moves see the `grounding` finding).
Four panels: (A) forward-KL trajectories (dry climbs, grounded suppressed); (B) H trajectories
(barely moves smoothing resists H-collapse); (C) stationary forward-KL vs g; (D) tail-item
survival vs g. Reads only the committed bundle.
Usage: python figures/plot_collapse.py [results/collapse]
"""
from __future__ import annotations
import sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
sys.path.insert(0, str(Path(__file__).parent))
from _figlib import load_bundle, mean_ci, savefig # noqa: E402
sys.path.insert(0, str(Path(__file__).parents[1] / "src"))
from knowledge.metrics import heterozygosity # noqa: E402
from neural.config import SyntheticCfg # noqa: E402
from neural.synthetic import make_mode_truth # noqa: E402
def main(results_dir: str = "results/collapse") -> None:
df, cfg = load_bundle(results_dir)
syn = SyntheticCfg(**cfg["synthetic"])
H_star = heterozygosity(make_mode_truth(syn).p_star)
g_values = sorted(df["g"].unique())
last = int(cfg["generations"] * 0.6)
stat = df[df["generation"] >= last]
colors = plt.cm.viridis(np.linspace(0, 0.85, len(g_values)))
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
# Panel A: forward-KL trajectories (dry climbs, grounded suppressed).
ax = axes[0, 0]
for g, c in zip(g_values, colors):
s = df[df["g"] == g].groupby("generation")["forward_kl"].mean()
ax.plot(s.index, s.values, "-o", color=c, ms=3, label=f"g={g:g}")
ax.set(xlabel="generation", ylabel=r"forward-KL $D(p^*\Vert\hat p)$",
title="Collapse in weights: dry KL climbs, grounding holds it")
ax.legend(frameon=False, fontsize=9)
# Panel B: H trajectories (barely moves — smoothing resists H-collapse).
ax = axes[0, 1]
for g, c in zip(g_values, colors):
s = df[df["g"] == g].groupby("generation")["heterozygosity"].mean()
ax.plot(s.index, s.values, "-o", color=c, ms=3, label=f"g={g:g}")
ax.axhline(H_star, ls=":", color="gray", lw=1, label="$H^*$")
ax.set(xlabel="generation", ylabel="heterozygosity $H$",
title="H barely moves (RNN smoothing resists H-collapse)")
ax.legend(frameon=False, fontsize=9)
# Panel C: stationary forward-KL vs g.
ax = axes[1, 0]
kg, Km, Kci = mean_ci(stat, "g", "forward_kl")
ax.errorbar(kg, Km, yerr=Kci, fmt="o-", color="#1f77b4", capsize=3)
ax.set(xlabel="grounding fraction $g$", ylabel=r"stationary forward-KL",
title="Grounding lowers stationary divergence")
# Panel D: tail-item survival vs g.
ax = axes[1, 1]
tg, Tm, Tci = mean_ci(stat, "g", "tail_frac_alive")
ax.errorbar(tg, Tm, yerr=Tci, fmt="s-", color="#d62728", capsize=3)
ax.set(xlabel="grounding fraction $g$", ylabel="tail items alive",
title="Grounding lifts tail survival")
fig.suptitle("collapse — a trained GRU collapses under dry self-training; grounding arrests it "
f"($K$={syn.K}, $n$={cfg['dynamics']['n']})", y=1.0, fontsize=13)
fig.tight_layout()
savefig(fig, results_dir, "collapse")
if __name__ == "__main__":
main(*sys.argv[1:])

142
figures/plot_grounding.py Normal file
View file

@ -0,0 +1,142 @@
"""`grounding` figure — the neural grounding response in REAL weights (↔ Layer-1 E2).
Honest reframing (see the progress log): the RNN's smoothing inductive bias makes H and
tail-SURVIVAL the *wrong* neural collapse metrics the model keeps spurious tail support
alive even while its distribution drifts far from truth, so tail_truth_mass_alive is flat /
non-monotone in g. The operative neural collapse metric is FORWARD-KL, on which grounding's
effect is monotone and significant. The sharp Layer-1 threshold (g*=0.048) is an *exact-
operator* feature reproduced quantitatively by the histogram bridge (g*=0.047); the trained
RNN confirms it in SIGN and softens it in sharpness.
Four panels: (A) forward-KL trajectories (dry climbs, grounded suppressed); (B) the phase
boundary stationary forward-KL vs g, monotone down; (C) recovery fraction with the median-
recovery grounding (Layer-1's 0.048) and the note that full recovery needs much more g in a
smoothing model; (D) the metric-choice panel H and tail-survival are flat/non-monotone
while forward-KL responds. Reads only the committed bundle.
Usage: python figures/plot_grounding.py [results/grounding]
"""
from __future__ import annotations
import sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
sys.path.insert(0, str(Path(__file__).parent))
from _figlib import load_bundle, mean_ci, savefig # noqa: E402
sys.path.insert(0, str(Path(__file__).parents[1] / "src"))
from knowledge.analysis import reduce_to_stationary # noqa: E402
from knowledge.metrics import heterozygosity # noqa: E402
from neural.config import SyntheticCfg # noqa: E402
from neural.synthetic import make_mode_truth # noqa: E402
_LAYER1_GSTAR = 0.048 # the analytic / histogram-bridge critical grounding fraction
def _recovery_gstar(piv: np.ndarray, gs: np.ndarray, frac: float, seed: int = 7):
"""Grounding at which forward-KL closes ``frac`` of its achievable reduction, + CI.
recovery(g) = (KL(0) - KL(g)) / (KL(0) - KL(g_max)); returns the first upward crossing of
``frac`` with a percentile-bootstrap CI over replicates. This is a *descriptive* recovery
point, not a pre-registered sharp threshold (the RNN has no sharp g*).
"""
def crossing(mat):
mk = mat.mean(axis=0)
rec = (mk[0] - mk) / (mk[0] - mk[-1])
idx = np.where(rec >= frac)[0]
if idx.size == 0 or idx[0] == 0:
return np.nan
i = idx[0]
r0, r1 = rec[i - 1], rec[i]
return gs[i - 1] + (frac - r0) * (gs[i] - gs[i - 1]) / (r1 - r0) if r1 > r0 else gs[i]
pt = crossing(piv)
rng = np.random.default_rng(seed)
boots = np.array([crossing(piv[rng.integers(0, piv.shape[0], piv.shape[0])])
for _ in range(4000)])
boots = boots[~np.isnan(boots)]
lo, hi = (np.percentile(boots, [2.5, 97.5]) if boots.size else (np.nan, np.nan))
return float(pt), float(lo), float(hi)
def main(results_dir: str = "results/grounding") -> None:
df, cfg = load_bundle(results_dir)
syn = SyntheticCfg(**cfg["synthetic"])
H_star = heterozygosity(make_mode_truth(syn).p_star)
g_values = np.array(sorted(df["g"].unique()))
last = int(cfg["generations"] * 0.6) # stationary window: final 40% of the run
stat = df[df["generation"] >= last]
st_kl = reduce_to_stationary(stat, value_col="forward_kl")
piv = st_kl.pivot(index="replicate", columns="g", values="forward_kl")[g_values].to_numpy()
g50, lo50, hi50 = _recovery_gstar(piv, g_values, 0.5)
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
colors = plt.cm.viridis(np.linspace(0, 0.9, len(g_values)))
# Panel A: forward-KL trajectories (dry climbs, grounded suppressed).
ax = axes[0, 0]
for g, c in zip(g_values, colors):
s = df[df["g"] == g].groupby("generation")["forward_kl"].mean()
ax.plot(s.index, s.values, color=c, label=f"g={g:g}")
ax.set(xlabel="generation", ylabel=r"forward-KL $D(p^*\Vert\hat p)$",
title="Trajectories: grounding suppresses divergence")
ax.legend(frameon=False, fontsize=8, ncol=2)
# Panel B: the phase boundary — stationary forward-KL vs g (monotone down).
ax = axes[0, 1]
kg, Km, Kci = mean_ci(st_kl, "g", "forward_kl")
ax.errorbar(kg, Km, yerr=Kci, fmt="o-", color="#1f77b4", capsize=3, zorder=3)
ax.set(xlabel="grounding fraction $g=m/(n+m)$",
ylabel=r"stationary forward-KL $D(p^*\Vert\hat p)$",
title="Neural phase boundary: KL falls monotonically\n"
"(sign confirmed; paired $t$=3.3 at g=0.2)")
# Panel C: recovery fraction with the median-recovery grounding vs Layer-1's g*.
ax = axes[1, 0]
mk = piv.mean(axis=0)
rec = (mk[0] - mk) / (mk[0] - mk[-1])
ax.plot(g_values, rec, "o-", color="#2ca02c")
ax.axhline(0.5, ls=":", color="gray", lw=1)
ax.axvspan(lo50, hi50, color="#d62728", alpha=0.15)
ax.axvline(g50, color="#d62728", lw=1.2,
label=f"median-recovery $g$={g50:.3f}\n(95% CI [{lo50:.3f},{hi50:.3f}])")
ax.axvline(_LAYER1_GSTAR, ls="--", color="k", lw=1, label=f"Layer-1 $g^*$={_LAYER1_GSTAR}")
ax.set(xlabel="grounding fraction $g$", ylabel="forward-KL recovery fraction",
title="Half the divergence gap closes by $g\\approx0.04$\n"
"(full recovery needs more g: smoothing softens the threshold)")
ax.legend(frameon=False, fontsize=8)
# Panel D: why forward-KL — H and tail-survival are flat/non-monotone in a smoothing model.
ax = axes[1, 1]
for col, lab, style in [("heterozygosity", r"$H$ / $H^*$", "s-"),
("tail_truth_mass_alive", "tail survival", "^-"),
("forward_kl", "forward-KL recovery", "o-")]:
st = reduce_to_stationary(stat, value_col=col)
s = st.groupby("g")[col].mean().reindex(g_values)
if col == "heterozygosity":
y = s.to_numpy() / H_star
elif col == "forward_kl":
y = (s.iloc[0] - s.to_numpy()) / (s.iloc[0] - s.iloc[-1]) # recovery, 0..1
else:
y = s.to_numpy()
ax.plot(g_values, y, style, ms=4, label=lab)
ax.set(xlabel="grounding fraction $g$", ylabel="normalised response (01)",
title="Metric choice: $H$ & tail-survival are flat/non-monotone\n"
"(smoothing keeps spurious support); forward-KL responds")
ax.legend(frameon=False, fontsize=8)
fig.suptitle("grounding — grounding arrests collapse in trained RNN weights (SIGN confirmed); "
f"the sharp $g^*\\ll1$ is carried by the histogram bridge ($g^*$=0.047, $H^*$={H_star:.2f})",
y=1.0, fontsize=12)
fig.tight_layout()
savefig(fig, results_dir, "grounding")
if __name__ == "__main__":
main(*sys.argv[1:])

View file

@ -0,0 +1,90 @@
"""`recombination` figure — the E4 "merge, don't average" finding in REAL weights (↔ E4/C4).
K_T specialist GRUs are trained on assignments from the exact shared-switch retention
construction, so union coverage matches the closed form U(K_T,rho,q) exactly. The pupil then
recombines the *measured* teacher distributions two ways: mean (naive pooling) vs oracle-guided
max-merge (M2N2-style union). The Layer-1 conservation law mean-mixture is flat in K_T while
max-merge rises is what must survive the move to trained weights.
Four panels: (A) union coverage vs K_T with the closed-form overlay (recombination *supply*);
(B) analytic teachers max-merge rises, mean stays flat (the conservation law); (C) the same
on *trained* teacher weights (same signs, compressed/noisier the smoothing caveat); (D) the
rho=1 control identical teachers buy nothing under either operator. Reads only the bundle.
Usage: python figures/plot_recombination.py [results/recombination]
"""
from __future__ import annotations
import sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
sys.path.insert(0, str(Path(__file__).parent))
from _figlib import load_bundle, mean_ci, savefig # noqa: E402
def U_closed(K_T, rho, q):
return rho * q + (1 - rho) * (1 - (1 - q) ** K_T)
def main(results_dir: str = "results/recombination") -> None:
df, cfg = load_bundle(results_dir)
q = cfg["coverage"]["q"]
K_Ts = np.array(sorted(df["K_T"].unique()))
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
# Panel A: union coverage vs K_T (rho=0 and rho=1) + closed-form overlay.
ax = axes[0, 0]
for rho, col in [(0.0, "#1f77b4"), (1.0, "#ff7f0e")]:
sub = df[df["rho"] == rho]
_, m, ci = mean_ci(sub, "K_T", "union_coverage")
ax.errorbar(sorted(sub["K_T"].unique()), m, yerr=ci, fmt="o", color=col, capsize=3,
label=fr"union ($\rho$={rho:g})")
ax.plot(K_Ts, [U_closed(K, rho, q) for K in K_Ts], "-", color=col, lw=1)
ax.set(xlabel="number of teachers $K_T$", ylabel="union tail coverage",
title=r"Supply: union matches closed form $U(K_T,\rho,q)$")
ax.legend(frameon=False, fontsize=9)
# Panels B/C: mean vs max surviving coverage at rho=0, analytic then trained.
def operator_panel(ax, mean_col, max_col, title):
r0 = df[df["rho"] == 0.0]
for col, style, lab, c in [(max_col, "-o", "max-merge (union)", "#2ca02c"),
(mean_col, "--s", "mean-distill (pool)", "#d62728")]:
_, m, ci = mean_ci(r0, "K_T", col)
ax.errorbar(sorted(r0["K_T"].unique()), m, yerr=ci, fmt=style, color=c, capsize=3,
label=lab)
ax.set(xlabel="number of teachers $K_T$", ylabel="surviving tail coverage", title=title)
ax.legend(frameon=False, fontsize=9)
operator_panel(axes[0, 1], "surviving_mean_target", "surviving_max_target",
r"Analytic teachers ($\rho$=0): max rises, mean flat"
"\n(the conservation law)")
operator_panel(axes[1, 0], "surviving_mean", "surviving_max",
r"Trained GRU teachers ($\rho$=0): same signs"
"\n(compressed + noisier — smoothing caveat)")
# Panel D: rho=1 control — identical teachers, both operators flat.
ax = axes[1, 1]
r1 = df[df["rho"] == 1.0]
for col, style, lab, c in [("surviving_max", "-o", "max-merge", "#2ca02c"),
("surviving_mean", "--s", "mean-distill", "#d62728"),
("union_coverage", ":^", "union", "#1f77b4")]:
_, m, ci = mean_ci(r1, "K_T", col)
ax.errorbar(sorted(r1["K_T"].unique()), m, yerr=ci, fmt=style, color=c, capsize=3, label=lab)
ax.set(xlabel="number of teachers $K_T$", ylabel="tail coverage",
title=r"Control ($\rho$=1, identical teachers):"
"\nmore teachers buy nothing")
ax.legend(frameon=False, fontsize=9)
fig.suptitle("recombination — 'merge, don't average' holds in real weights: max-merge realises "
"the multi-teacher tail benefit, mean-distill conserves collapse", y=1.0, fontsize=12)
fig.tight_layout()
savefig(fig, results_dir, "recombination")
if __name__ == "__main__":
main(*sys.argv[1:])