Blame for protocol_engine.py

#!/usr/bin/env python3
"""
protocol_engine.py -- protocol design engine. Python standard library only.

Input : one JSON config describing a registered construct test
        (stages/arms, effect thresholds, alpha, target power, gate rules).
Output: a markdown design report (stdout) containing, per arm:
          1. required sample size, closed form;
          2. Monte Carlo power at that size (cross-check of the closed form);
          3. power at the currently registered size;
        plus a falsifier-fire audit, an input-separation gate budget,
        a scale table, and a registration skeleton filled from the config.

Deterministic: one seed in the config; every simulation stream derives
from it. Re-running the same command reproduces every number.

Run:  python3 protocol_engine.py config_summarization_bias.json
"""

import json
import math
import random
import sys
from statistics import NormalDist

ND = NormalDist()


# ----------------------------------------------------------------------
# Helpers
# ----------------------------------------------------------------------

def z_alpha(alpha, sides):
    return ND.inv_cdf(1 - alpha / sides)


def z_power(power):
    return ND.inv_cdf(power)


def cohens_h(p1, p2):
    return 2 * math.asin(math.sqrt(p1)) - 2 * math.asin(math.sqrt(p2))


def rate_from_h(p_base, h):
    """Rate p2 > p_base such that cohens_h(p2, p_base) == h."""
    return math.sin(math.asin(math.sqrt(p_base)) + h / 2.0) ** 2


def rate_from_or(p_base, odds_ratio):
    odds = odds_ratio * p_base / (1.0 - p_base)
    return odds / (1.0 + odds)


def binom_pmf(k, n, p):
    return math.comb(n, k) * p ** k * (1.0 - p) ** (n - k)


def binom_cdf(k, n, p):
    if k < 0:
        return 0.0
    if k >= n:
        return 1.0
    return sum(binom_pmf(i, n, p) for i in range(k + 1))


def binom_draw(rng, n, p):
    if hasattr(rng, "binomialvariate"):  # Python >= 3.12
        return rng.binomialvariate(n, p)
    return sum(1 for _ in range(n) if rng.random() < p)


def fmt(x, digits=3):
    return f"{x:.{digits}f}"


# ----------------------------------------------------------------------
# Closed-form sample sizes
# ----------------------------------------------------------------------

def n_two_proportions(h, alpha, sides, power):
    """Per-group n for a two-proportion test, arcsine (Cohen's h) form."""
    z = z_alpha(alpha, sides) + z_power(power)
    return math.ceil(2.0 * (z / h) ** 2)


def n_one_proportion(p0, p1, alpha, sides, power):
    za = z_alpha(alpha, sides)
    zb = z_power(power)
    num = za * math.sqrt(p0 * (1 - p0)) + zb * math.sqrt(p1 * (1 - p1))
    return math.ceil((num / abs(p1 - p0)) ** 2)


def _beta_cf(a, b, x):
    """Continued fraction for the incomplete beta function."""
    qab, qap, qam = a + b, a + 1.0, a - 1.0
    c = 1.0
    d = 1.0 - qab * x / qap
    if abs(d) < 1e-300:
        d = 1e-300
    d = 1.0 / d
    h = d
    for m in range(1, 201):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d
        if abs(d) < 1e-300:
            d = 1e-300
        c = 1.0 + aa / c
        if abs(c) < 1e-300:
            c = 1e-300
        d = 1.0 / d
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d
        if abs(d) < 1e-300:
            d = 1e-300
        c = 1.0 + aa / c
        if abs(c) < 1e-300:
            c = 1e-300
        d = 1.0 / d
        delta = d * c
        h *= delta
        if abs(delta - 1.0) < 3e-16:
            break
    return h


def _reg_incomplete_beta(x, a, b):
    if x <= 0.0:
        return 0.0
    if x >= 1.0:
        return 1.0
    bt = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
                  + a * math.log(x) + b * math.log(1.0 - x))
    if x < (a + 1.0) / (a + b + 2.0):
        return bt * _beta_cf(a, b, x) / a
    return 1.0 - bt * _beta_cf(b, a, 1.0 - x) / b


def t_cdf(t, df):
    """CDF of the central Student t distribution (incomplete beta)."""
    x = df / (df + t * t)
    ib = _reg_incomplete_beta(x, df / 2.0, 0.5)
    return 1.0 - 0.5 * ib if t >= 0 else 0.5 * ib


def t_crit(alpha, sides, df):
    target = 1.0 - alpha / sides
    lo, hi = 0.0, 1000.0
    while hi - lo > 1e-12:
        mid = (lo + hi) / 2.0
        if t_cdf(mid, df) < target:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2.0


def _chi_pdf(v, df):
    if v <= 0.0:
        return 0.0
    k = df / 2.0
    return math.exp((k - 1.0) * math.log(v) - v / 2.0
                    - k * math.log(2.0) - math.lgamma(k))


def _adaptive_simpson(f, a, b, eps, depth):
    c = (a + b) / 2.0
    fa, fb, fc = f(a), f(b), f(c)
    whole = (b - a) * (fa + 4.0 * fc + fb) / 6.0

    def step(a, b, eps, whole, fa, fb, fc, depth):
        c = (a + b) / 2.0
        lm, rm = (a + c) / 2.0, (c + b) / 2.0
        flm, frm = f(lm), f(rm)
        left = (c - a) * (fa + 4.0 * flm + fc) / 6.0
        right = (b - c) * (fc + 4.0 * frm + fb) / 6.0
        delta = left + right - whole
        if depth <= 0 or abs(delta) <= 15.0 * eps:
            return left + right + delta / 15.0
        return (step(a, c, eps / 2.0, left, fa, fc, flm, depth - 1)
                + step(c, b, eps / 2.0, right, fc, fb, frm, depth - 1))

    return step(a, b, eps, whole, fa, fb, fc, depth)


def nct_cdf(t, df, ncp):
    """CDF of the noncentral t distribution, from its definition
    T = (Z + ncp) / sqrt(V/df): P(T <= t) = E_V[Phi(t sqrt(V/df) - ncp)]
    over the chi-square density of V, by adaptive Simpson quadrature."""
    if t < 0:
        return 1.0 - nct_cdf(-t, df, -ncp)
    f = lambda v: ND.cdf(t * math.sqrt(v / df) - ncp) * _chi_pdf(v, df)
    hi = df + 20.0 * math.sqrt(2.0 * df) + 60.0
    return _adaptive_simpson(f, 0.0, hi, 1e-12, 40)


def paired_t_power_exact(n, dz, alpha, sides):
    """Exact power of the paired (one-sample) t test on differences:
    noncentral t at df = n-1, noncentrality dz*sqrt(n)."""
    df = float(n - 1)
    ncp = dz * math.sqrt(n)
    crit = t_crit(alpha, sides, df)
    power = 1.0 - nct_cdf(crit, df, ncp)
    if sides == 2:
        power += nct_cdf(-crit, df, ncp)
    return power


def n_paired_dz(dz, alpha, sides, power):
    """Smallest n whose exact paired-t power reaches the target --
    an exact search over the noncentral t power, not the normal
    closed form."""
    n = 2
    while paired_t_power_exact(n, dz, alpha, sides) < power:
        n += 1
    return n


# ----------------------------------------------------------------------
# Two-proportion z test: p-value, Monte Carlo power, exact power
# ----------------------------------------------------------------------

def two_prop_z_pvalue(x1, n1, x2, n2):
    pooled = (x1 + x2) / (n1 + n2)
    se = math.sqrt(pooled * (1 - pooled) * (1.0 / n1 + 1.0 / n2))
    if se == 0.0:
        return 1.0
    z = (x1 / n1 - x2 / n2) / se
    return 2.0 * (1.0 - ND.cdf(abs(z)))


def mc_power_two_prop(n, p1, p2, alpha, reps, seed):
    rng = random.Random(seed)
    hits = 0
    for _ in range(reps):
        x = binom_draw(rng, n, p1)
        y = binom_draw(rng, n, p2)
        if two_prop_z_pvalue(x, n, y, n) < alpha:
            hits += 1
    p = hits / reps
    return p, math.sqrt(p * (1 - p) / reps)


def exact_power_two_prop_z(n, p1, p2, alpha):
    """Exact power of the pooled z test by enumerating all outcome pairs.
    Practical for small n (used at the registered n = 20)."""
    probs1 = [binom_pmf(x, n, p1) for x in range(n + 1)]
    probs2 = [binom_pmf(y, n, p2) for y in range(n + 1)]
    total = 0.0
    for x in range(n + 1):
        for y in range(n + 1):
            if two_prop_z_pvalue(x, n, y, n) < alpha:
                total += probs1[x] * probs2[y]
    return total


def fisher_pvalue(x, n, y, m):
    """Two-sided Fisher exact p-value for [[x, n-x], [y, m-y]]."""
    s = x + y
    lo = max(0, s - m)
    hi = min(n, s)
    denom = math.comb(n + m, s)

    def hyper(k):
        return math.comb(n, k) * math.comb(m, s - k) / denom

    p_obs = hyper(x)
    return sum(hyper(k) for k in range(lo, hi + 1)
               if hyper(k) <= p_obs * (1 + 1e-12))


def exact_power_fisher(n, p1, p2, alpha):
    probs1 = [binom_pmf(x, n, p1) for x in range(n + 1)]
    probs2 = [binom_pmf(y, n, p2) for y in range(n + 1)]
    total = 0.0
    for x in range(n + 1):
        for y in range(n + 1):
            if fisher_pvalue(x, n, y, n) < alpha:
                total += probs1[x] * probs2[y]
    return total


# ----------------------------------------------------------------------
# One-sample binomial test: exact power
# ----------------------------------------------------------------------

def one_sample_critical(n, p0, alpha, sides):
    """Equal-tail exact critical counts. Reject if X <= lo or X >= hi."""
    if sides == 2:
        lo = -1
        for k in range(n + 1):
            if binom_cdf(k, n, p0) <= alpha / 2:
                lo = k
            else:
                break
        hi = n + 1
        for k in range(n, -1, -1):
            if 1.0 - binom_cdf(k - 1, n, p0) <= alpha / 2:
                hi = k
            else:
                break
        return lo, hi
    hi = n + 1
    for k in range(n, -1, -1):
        if 1.0 - binom_cdf(k - 1, n, p0) <= alpha:
            hi = k
        else:
            break
    return -1, hi


def exact_power_one_sample(n, p0, p1, alpha, sides):
    lo, hi = one_sample_critical(n, p0, alpha, sides)
    power = 0.0
    if lo >= 0:
        power += binom_cdf(lo, n, p1)
    if hi <= n:
        power += 1.0 - binom_cdf(hi - 1, n, p1)
    return power


# ----------------------------------------------------------------------
# Paired-difference test (Stage 2): Monte Carlo power
# ----------------------------------------------------------------------

def mc_power_paired_dz(n, dz, alpha, sides, reps, seed):
    """Independent check of the exact paired-t size: simulate the
    paired t test itself. Differences D ~ Normal(dz, 1) are drawn,
    the variance is ESTIMATED from each sample (sample sd, ddof=1),
    and the statistic is compared with the t critical value -- the
    t behaviour the closed form must answer to. Because the check
    shares no formula with n_paired_dz, it can falsify it: run at
    the old normal-approximation size n = 32 (dz = 0.5, alpha 0.05)
    it returns ~0.78, below the 0.80 target the old size claimed."""
    rng = random.Random(seed)
    crit = t_crit(alpha, sides, float(n - 1))
    hits = 0
    for _ in range(reps):
        diffs = [rng.gauss(dz, 1.0) for _ in range(n)]
        mean = sum(diffs) / n
        var = sum((d - mean) ** 2 for d in diffs) / (n - 1)
        if var <= 0.0:
            continue
        t = mean / math.sqrt(var / n)
        if (sides == 2 and abs(t) > crit) or (sides == 1 and t > crit):
            hits += 1
    p = hits / reps
    return p, math.sqrt(p * (1 - p) / reps)


# ----------------------------------------------------------------------
# Input-separation gate budget
# ----------------------------------------------------------------------

def gate_build_count(n_required, fail_rate, assurance):
    """Smallest build count B with P(Binomial(B, 1-fail) >= n_required)
    >= assurance. Exact search over the binomial CDF."""
    p_pass = 1.0 - fail_rate
    B = n_required
    while True:
        if 1.0 - binom_cdf(n_required - 1, B, p_pass) >= assurance:
            return B
        B += 1


def mc_gate_check(B, n_required, fail_rate, reps, seed):
    rng = random.Random(seed)
    hits = sum(1 for _ in range(reps)
               if binom_draw(rng, B, 1.0 - fail_rate) >= n_required)
    p = hits / reps
    return p, math.sqrt(p * (1 - p) / reps)


def exact_search_one_sample(p0, p1, alpha, sides, target, n_start):
    """Smallest n >= n_start whose exact binomial power reaches target.
    The exact test is conservative at small n, so this sits at or above
    the closed-form value. Search is capped at n_start + 120."""
    n = n_start
    while n <= n_start + 120:
        if exact_power_one_sample(n, p0, p1, alpha, sides) >= target:
            return n
        n += 1
    return n_start


# ----------------------------------------------------------------------
# Standards
#   standard : the config's alpha and target power.
#   decisive : alpha = 0.01, power = 0.95 -- the scale at which a verdict,
#              either way, is hard to argue with.
# ----------------------------------------------------------------------

DECISIVE_ALPHA = 0.01
DECISIVE_POWER = 0.95


def run(cfg):
    alpha = cfg["alpha"]
    sides = cfg.get("sides", 2)
    power_t = cfg["target_power"]
    reps = cfg.get("mc_reps", 20000)
    seed0 = cfg.get("seed", 1)
    out = []
    out.append(f"# Design report -- {cfg['construct_name']}")
    out.append("")
    out.append(f"Source protocol: {cfg['source_protocol']}")
    out.append("")
    out.append(f"- Standard design point: alpha = {alpha} ({sides}-sided), "
               f"power = {power_t}")
    out.append(f"- Decisive design point: alpha = {DECISIVE_ALPHA} "
               f"({sides}-sided), power = {DECISIVE_POWER}")
    out.append(f"- Monte Carlo reps per cell = {reps}; base seed = {seed0}; "
               f"one stream per row. Re-running reproduces every number.")
    out.append("")

    # ---- Table 1: per-arm sizing at both design points ----------------
    out.append("## 1. Required sample size per arm")
    out.append("")
    out.append("| Arm | Test | Effect as specified | n at standard point | "
               "MC power at that n | n at decisive point | Current n | "
               "Power at current n |")
    out.append("|---|---|---|---|---|---|---|---|")
    res = {}

    # Stage 3 (core): two-proportion, model vs human told-preference rate.
    s3 = cfg["stage3"]
    p_base = s3["human_told_preference_rate"]
    p_h = rate_from_h(p_base, s3["registered_thresholds"]["cohens_h"])
    p_or = rate_from_or(p_base, s3["registered_thresholds"]["odds_ratio"])
    h_or = abs(cohens_h(p_or, p_base))
    h_gov = min(s3["registered_thresholds"]["cohens_h"], h_or)
    p_alt = p_h if s3["registered_thresholds"]["cohens_h"] <= h_or else p_or
    n_std = n_two_proportions(h_gov, alpha, sides, power_t)
    n_dec = n_two_proportions(h_gov, DECISIVE_ALPHA, sides, DECISIVE_POWER)
    mc_std, se_std = mc_power_two_prop(n_std, p_alt, p_base, alpha,
                                       reps, seed0 + 101)
    n_cur = s3["current_pairs"]
    mc_cur, se_cur = mc_power_two_prop(n_cur, p_alt, p_base, alpha,
                                       reps, seed0 + 102)
    res["stage3"] = dict(n_std=n_std, n_dec=n_dec, p_alt=p_alt, p_base=p_base,
                         h=h_gov)
    out.append(
        f"| Stage 3, evaluative (core) | two-proportion z, pooled | "
        f"h = {s3['registered_thresholds']['cohens_h']} -> rate {p_base} vs "
        f"{fmt(p_h)}; OR = {s3['registered_thresholds']['odds_ratio']} -> rate "
        f"{fmt(p_or)} (h = {fmt(h_or)}); governing h = {fmt(h_gov)} | "
        f"{n_std} per group | {fmt(mc_std)} +/- {fmt(se_std)} | {n_dec} per group | "
        f"{n_cur} per group | {fmt(mc_cur)} +/- {fmt(se_cur)} |")

    # Stage 2: paired SI difference.
    s2 = cfg["stage2"]
    dz = s2["planning_effect_dz"]
    n2_std = n_paired_dz(dz, alpha, sides, power_t)
    n2_dec = n_paired_dz(dz, DECISIVE_ALPHA, sides, DECISIVE_POWER)
    mc2, se2 = mc_power_paired_dz(n2_std, dz, alpha, sides, reps, seed0 + 201)
    mc2c, se2c = mc_power_paired_dz(s2["current_matrices"], dz, alpha, sides,
                                    reps, seed0 + 202)
    res["stage2"] = dict(n_std=n2_std, n_dec=n2_dec)
    out.append(
        f"| Stage 2, generative | paired t on SI differences | dz = {dz} "
        f"(worked value -- the paper fixes no SI effect size) | {n2_std} pairs | "
        f"{fmt(mc2)} +/- {fmt(se2)} | {n2_dec} pairs | "
        f"{s2['current_matrices']} pairs | {fmt(mc2c)} +/- {fmt(se2c)} |")

    # Density calibration arm: one-sample vs chance.
    da = cfg["density_arm"]
    p1d = rate_from_h(da["null_preference_rate"],
                      da["planning_effect_h_vs_null"])
    nd_cf = n_one_proportion(da["null_preference_rate"], p1d, alpha, sides,
                             power_t)
    nd_std = exact_search_one_sample(da["null_preference_rate"], p1d, alpha,
                                     sides, power_t, nd_cf)
    nd_cf_dec = n_one_proportion(da["null_preference_rate"], p1d,
                                 DECISIVE_ALPHA, sides, DECISIVE_POWER)
    nd_dec = exact_search_one_sample(da["null_preference_rate"], p1d,
                                     DECISIVE_ALPHA, sides, DECISIVE_POWER,
                                     nd_cf_dec)
    exd = exact_power_one_sample(nd_std, da["null_preference_rate"], p1d,
                                 alpha, sides)
    exdd = exact_power_one_sample(nd_dec, da["null_preference_rate"], p1d,
                                  DECISIVE_ALPHA, sides)
    exdc = exact_power_one_sample(20, da["null_preference_rate"], p1d,
                                  alpha, sides)
    res["density"] = dict(n_std=nd_std, n_dec=nd_dec)
    out.append(
        f"| Density calibration arm | one-sample exact binomial vs chance | "
        f"h = {da['planning_effect_h_vs_null']} -> rate "
        f"{da['null_preference_rate']} vs {fmt(p1d)} (exact search on "
        f"the unrounded rate {fmt(p1d, 5)}) | {nd_std} pairs "
        f"(exact power {fmt(exd)}; closed form {nd_cf}) | -- | {nd_dec} pairs "
        f"(exact power {fmt(exdd)}; closed form {nd_cf_dec}) | 20 pairs "
        f"(exact power {fmt(exdc)}) |")

    # Reader validation: one-sample convergence vs chance.
    rv = cfg["reader_validation"]
    nr_cf = n_one_proportion(rv["null_convergence_rate"],
                             rv["planning_target_rate"], alpha, sides, power_t)
    nr_std = exact_search_one_sample(rv["null_convergence_rate"],
                                     rv["planning_target_rate"], alpha,
                                     sides, power_t, nr_cf)
    nr_cf_dec = n_one_proportion(rv["null_convergence_rate"],
                                 rv["planning_target_rate"],
                                 DECISIVE_ALPHA, sides, DECISIVE_POWER)
    nr_dec = exact_search_one_sample(rv["null_convergence_rate"],
                                     rv["planning_target_rate"],
                                     DECISIVE_ALPHA, sides, DECISIVE_POWER,
                                     nr_cf_dec)
    exr = exact_power_one_sample(nr_std, rv["null_convergence_rate"],
                                 rv["planning_target_rate"], alpha, sides)
    exrd = exact_power_one_sample(nr_dec, rv["null_convergence_rate"],
                                  rv["planning_target_rate"],
                                  DECISIVE_ALPHA, sides)
    exrc = exact_power_one_sample(20, rv["null_convergence_rate"],
                                  rv["planning_target_rate"], alpha, sides)
    res["reader"] = dict(n_std=nr_std, n_dec=nr_dec)
    out.append(
        f"| SI reader validation | one-sample exact binomial vs chance | "
        f"rate {rv['null_convergence_rate']} vs "
        f"{rv['planning_target_rate']} (worked values -- chance depends on the "
        f"coding scheme) | {nr_std} items (exact power {fmt(exr)}; closed "
        f"form {nr_cf}) | -- | {nr_dec} items (exact power {fmt(exrd)}; "
        f"closed form {nr_cf_dec}) | 20 items "
        f"(exact power {fmt(exrc)}) |")

    # Register arm: two-proportion across registers, per stratum.
    ra = cfg["register_arm"]
    p_alt_r = rate_from_h(ra["human_baseline_rate"],
                          ra["planning_effect_h_between_registers"])
    nreg_std = n_two_proportions(ra["planning_effect_h_between_registers"],
                                 alpha, sides, power_t)
    nreg_dec = n_two_proportions(ra["planning_effect_h_between_registers"],
                                 DECISIVE_ALPHA, sides, DECISIVE_POWER)
    mcr, ser = mc_power_two_prop(nreg_std, p_alt_r, ra["human_baseline_rate"],
                                 alpha, reps, seed0 + 301)
    mcrc, serc = mc_power_two_prop(20, p_alt_r, ra["human_baseline_rate"],
                                   alpha, reps, seed0 + 302)
    res["register"] = dict(n_std=nreg_std, n_dec=nreg_dec)
    out.append(
        f"| Register arm | two-proportion z per register stratum | "
        f"h = {ra['planning_effect_h_between_registers']} -> rate "
        f"{ra['human_baseline_rate']} vs {fmt(p_alt_r)} | {nreg_std} contents "
        f"per stratum ({nreg_std * 4} versions) | {fmt(mcr)} +/- {fmt(ser)} | "
        f"{nreg_dec} contents per stratum | 20 contents per stratum | "
        f"{fmt(mcrc)} +/- {fmt(serc)} |")
    out.append("")

    # ---- Section 2: falsifier-fire audit, core arm at registered size --
    approx = ND.cdf(h_gov * math.sqrt(n_cur / 2) - z_alpha(alpha, sides))
    mc_f, mc_f_se = mc_power_two_prop(n_cur, p_alt, p_base, alpha,
                                      reps, seed0 + 999)
    ex_z = exact_power_two_prop_z(n_cur, p_alt, p_base, alpha)
    fisher = exact_power_fisher(n_cur, p_alt, p_base, alpha)
    one_s = exact_power_one_sample(n_cur, p_base, p_or, alpha, sides)
    paired_mc, _ = mc_power_paired_dz(n_cur,
                                      s3["registered_thresholds"]["cohens_h"],
                                      alpha, sides, reps, seed0 + 998)
    out.append("## 2. Falsifier-fire audit -- Stage 3 at the registered size")
    out.append("")
    out.append(f"Setup: {n_cur} pairs per group; true effect set exactly at the "
               f"registered threshold (h = {fmt(h_gov)}, rates {p_base} vs "
               f"{fmt(p_alt)}). 'Fires' = the test is not significant, so the "
               f"registered falsifier triggers against a true effect.")
    out.append("")
    out.append("| Reading | Test | Power | P(falsifier fires) |")
    out.append("|---|---|---|---|")
    out.append(f"| Registered (model rate vs human rate) | normal approx | "
               f"{fmt(approx)} | {fmt(1 - approx)} |")
    out.append(f"| Registered | pooled z, Monte Carlo | {fmt(mc_f)} +/- "
               f"{fmt(mc_f_se)} | {fmt(1 - mc_f)} |")
    out.append(f"| Registered | pooled z, exact enumeration | {fmt(ex_z)} | "
               f"{fmt(1 - ex_z)} |")
    out.append(f"| Registered | Fisher exact, exact enumeration | "
               f"{fmt(fisher)} | {fmt(1 - fisher)} |")
    out.append(f"| Alternate (model rate vs chance, OR threshold) | exact "
               f"binomial | {fmt(one_s)} | {fmt(1 - one_s)} |")
    out.append(f"| Alternate (paired per-item contrast) | paired t, Monte "
               f"Carlo | {fmt(paired_mc)} | {fmt(1 - paired_mc)} |")
    out.append("")
    fires = [1 - approx, 1 - mc_f, 1 - ex_z, 1 - fisher]
    out.append(f"Registered reading, across test variants: P(fire) = "
               f"{fmt(min(fires))} to {fmt(max(fires))}.")
    out.append("")

    # ---- Section 3: gate budget at both design points ------------------
    gate = cfg["gate"]
    for label, target in (("standard", n_std), ("decisive", n_dec)):
        out.append(f"## 3{label == 'decisive' and 'b' or 'a'}. "
                   f"{gate['name']} -- build budget ({label} point)")
        out.append("")
        out.append(f"Target: >= {target} pairs pass the gate. Assurance: "
                   f"P(passing >= {target}) >= {gate['assurance']}. "
                   f"Exact binomial search; last column re-checks by Monte Carlo.")
        out.append("")
        out.append("| Pair failure rate | Pairs to build (exact) | "
                   "Expected passing | MC assurance |")
        out.append("|---|---|---|---|")
        for j, f in enumerate(gate["failure_rate_scenarios"]):
            B = gate_build_count(target, f, gate["assurance"])
            mc_a, mc_a_se = mc_gate_check(B, target, f, reps,
                                          seed0 + 5000 + (0 if label == "standard" else 100) + j)
            out.append(f"| {f} | {B} | {fmt(B * (1 - f), 1)} | "
                       f"{fmt(mc_a)} +/- {fmt(mc_a_se)} |")
        out.append("")

    # ---- Section 4: scale table ----------------------------------------
    sc = cfg["scale"]
    judges = sc["model_families"] + sc["human_raters_per_pair"]
    out.append("## 4. Scale table")
    out.append("")
    out.append(f"Per pair: {sc['versions_per_pair']} stimulus versions "
               f"({sc['versions_per_pair'] * 2} with the register arm), "
               f"{judges} Stage 3 judge calls ({sc['model_families']} model "
               f"families + {sc['human_raters_per_pair']} human rater), "
               f"SI coding passes = versions x 2 raters. Stage 4 = two further "
               f"Stage 3 passes (length-mismatched, framing-varied).")
    out.append("")
    out.append("| Pairs | Versions | Versions with register arm | "
               "Stage 3 judge calls | Stage 4 judge calls | SI coding passes |")
    out.append("|---|---|---|---|---|---|")
    for P in sc["pair_counts"]:
        v = P * sc["versions_per_pair"]
        calls = P * judges
        out.append(f"| {P} | {v} | {v * 2} | {calls} | {calls * 2} | {v * 2} |")
    out.append("")

    # ---- Section 5: registration skeleton -------------------------------
    out.append("## 5. Registration skeleton (filled from this config)")
    out.append("")
    out.append("1. Construct: summarization bias -- the systematic tendency of a "
               "model, in generation or in evaluation, to replace or prefer "
               "reconstructable shown-mode structure with the abstract summary "
               "label that names its content (paper S4.3).")
    out.append(f"2. Thresholds, fixed in advance: model told-preference rate "
               f"exceeds human rate with Cohen's h >= "
               f"{s3['registered_thresholds']['cohens_h']} or odds ratio >= "
               f"{s3['registered_thresholds']['odds_ratio']} (paper S5.3). "
               f"Alpha = {alpha} ({sides}-sided) at the standard point; "
               f"alpha = {DECISIVE_ALPHA}, power = {DECISIVE_POWER} at the "
               f"decisive point.")
    out.append(f"3. Stimuli: build the gate-budget count from Table 3 for the "
               f"expected failure rate (standard point) -- word count +/-5%, "
               f"grade matched (paper S5.1).")
    out.append("4. Coding gate: SI coded by two independent raters per version; "
               "no category proceeds at Cohen's kappa < 0.60 (paper S5.1).")
    out.append(f"5. Input-separation gate: a pair enters Stage 3 only if its two "
               f"members separate on SI by the registered minimum; failed pairs "
               f"are counted, reported, and replaced from the build budget "
               f"(thread post #5).")
    out.append(f"6. Stage 2 (generative): >= {sc['model_families']} model "
               f"families, shown-mode instruction, SI compared to the human "
               f"shown target; n = {n2_std} pairs at the standard point, "
               f"{n2_dec} at the decisive point (Table 1).")
    out.append(f"7. Stage 3 (evaluative, core): pairs judged by the same models "
               f"and by a human panel, order counterbalanced; n = {n_std} pairs "
               f"per group at the standard point, {n_dec} at the decisive "
               f"point (Table 1).")
    out.append("8. Stage 4 (discriminant): Stage 3 re-run (a) length-mismatched "
               "and (b) framing-varied; a preference that vanishes under either "
               "is reclassified as that known bias (paper S5.4).")
    out.append(f"9. Added arms, sized in Table 1: density calibration "
               f"({nd_std} pairs), SI reader validation ({nr_std} items), "
               f"register arm ({nreg_std} contents per stratum). No arm is "
               f"registered without its own size from this table.")
    out.append("10. Decision rules (paper S5.5): Stage 3 supports and Stage 4 "
               "survives -> confirmed at the scope tested. Stage 3 falsifies "
               "-> construct withdrawn, published with the same prominence as "
               "a confirmation. Post-hoc qualifiers are prohibited.")
    out.append("")
    return "\n".join(out) + "\n"


KNOWN_SCENARIOS = ("null", "registered_alternative")


def check_scenario(cfg):
    """Refuse an unknown truth.scenario by name, on any config that
    carries one -- an unknown string must never silently select the
    world a verdict is about. Configs without a truth block (the
    sizing configs this program normally reads) are unaffected."""
    truth = cfg.get("truth")
    if isinstance(truth, dict) and "scenario" in truth:
        if truth["scenario"] not in KNOWN_SCENARIOS:
            sys.exit(f"config refused: truth.scenario "
                     f"{truth['scenario']!r} is not a known scenario "
                     f"(expected 'null' or 'registered_alternative')")


def main():
    if len(sys.argv) != 2:
        sys.exit("usage: python3 protocol_engine.py <config.json>")
    with open(sys.argv[1], "r", encoding="utf-8") as fh:
        cfg = json.load(fh)
    check_scenario(cfg)
    if "versions_per_pair" not in cfg.get("scale", {}):
        cfg["scale"]["versions_per_pair"] = 2
    if "mc_reps" not in cfg:
        cfg["mc_reps"] = 20000
    sys.stdout.write(run(cfg))


if __name__ == "__main__":
    main()