"""Statistician-lens checks for HUB-PAGE-v3 (pure stdlib)."""
import json
import math
from math import comb, sqrt

RECOUNT = "/workshop/bench-archive/plans-2026-09-28/decision-models/recount-v2.json"
d = json.load(open(RECOUNT))


def binom_two_sided_exact(b, c):
    """Exact McNemar: two-sided binomial on discordant pairs, doubling the smaller tail."""
    n = b + c
    k = min(b, c)
    tail = sum(comb(n, i) for i in range(0, k + 1)) / 2 ** n
    return min(1.0, 2 * tail)


def wilson(k, n, z=1.959964):
    p = k / n
    den = 1 + z * z / n
    centre = (p + z * z / (2 * n)) / den
    half = z * sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / den
    return centre - half, centre + half


def beta_inc_cp(k, n, alpha=0.05):
    """Clopper-Pearson by bisection on the binomial CDF."""
    def cdf(x, p):
        return sum(comb(n, i) * p ** i * (1 - p) ** (n - i) for i in range(0, x + 1))
    if k == 0:
        lo = 0.0
    else:
        a, b_ = 0.0, 1.0
        for _ in range(200):
            m = (a + b_) / 2
            # P(X >= k | m) = 1 - cdf(k-1)
            if 1 - cdf(k - 1, m) < alpha / 2:
                a = m
            else:
                b_ = m
        lo = (a + b_) / 2
    if k == n:
        hi = 1.0
    else:
        a, b_ = 0.0, 1.0
        for _ in range(200):
            m = (a + b_) / 2
            if cdf(k, m) > alpha / 2:
                a = m
            else:
                b_ = m
        hi = (a + b_) / 2
    return lo, hi


def newcombe_paired(b, c, n, z=1.959964):
    """Newcombe (1998) method 10 interval for a paired difference of proportions.
    b = right only in model A, c = right only in model B; a,d unknown individually but
    only marginals matter: need p1 = (a+b)/n, p2 = (a+c)/n. We need a (both right)."""
    raise NotImplementedError


def newcombe_paired_full(a, b, c, dd, z=1.959964):
    n = a + b + c + dd
    p1 = (a + b) / n
    p2 = (a + c) / n
    l1, u1 = wilson(a + b, n, z)
    l2, u2 = wilson(a + c, n, z)
    diff = p1 - p2
    # phi correlation
    num = a * dd - b * c
    den = sqrt((a + b) * (c + dd) * (a + c) * (b + dd))
    phi = 0 if den == 0 else num / den
    if num < 0:
        phi = phi  # Newcombe uses phi as-is (with a correction when num<0 not needed here)
    dl = sqrt((p1 - l1) ** 2 - 2 * phi * (p1 - l1) * (u2 - p2) + (u2 - p2) ** 2)
    du = sqrt((u1 - p1) ** 2 - 2 * phi * (u1 - p1) * (p2 - l2) + (p2 - l2) ** 2)
    return diff, diff - dl, diff + du


print("== 1. Net gap an exact McNemar needs for p < 0.05, by number of discordant lines ==")
for nd in (8, 14, 20, 24, 32, 37, 38, 42):
    need = None
    for minority in range(nd // 2, -1, -1):
        if binom_two_sided_exact(nd - minority, minority) < 0.05:
            need = (nd - minority, minority)
            break
    print(f"  discordant {nd}: first split under 0.05 is {need[0]}-{need[1]}, net {need[0]-need[1]} lines "
          f"= {100*(need[0]-need[1])/108:.1f} points of 108; p {binom_two_sided_exact(*need):.4f}")

print("\n  unpaired two-proportion 95% half-width at 63% on 108 v 108:",
      f"{100*1.959964*sqrt(2*0.63*0.37/108):.1f} points")
lo, hi = wilson(72, 108)
print(f"  single Wilson half-width for 72/108: {100*(hi-lo)/2:.1f} points")

print("\n== 2. Paired-difference intervals vs the naive Bayes (need both-right count) ==")
# both-right a = k_model - only_model_right; NB k = 68
s0 = d["s0"]
pairs = [("kev-9b", 72), ("imajev-4b-cpu", 68), ("apus-4b-high", 67), ("apus-9b-high", 66),
         ("lev-4b-gpu", 66), ("kev-4b", 60), ("mistral-small3.2-24b-generate-coi", 73)]
for m, k in pairs:
    mc = d["mcnemar"][f"{m}|naive-bayes"]
    b, c = mc["only_first_right"], mc["only_second_right"]
    a = k - b
    assert 68 - c == a, (m, a, 68 - c)
    dd = 108 - a - b - c
    diff, lo, hi = newcombe_paired_full(a, b, c, dd)
    print(f"  {m}: {k} v 68, split {b}-{c}, diff {100*diff:+.1f} pts, Newcombe paired 95% "
          f"{100*lo:+.1f} to {100*hi:+.1f}  (lines {108*lo:+.1f} to {108*hi:+.1f}); p {binom_two_sided_exact(b,c):.3f}")

print("\n== 3. Multiplicity over the NB table's 16 rows (+deem = 17) ==")
rows = ["openjev-bf16-readout-arm3", "openjev-q4km-generate-arm13", "jevify-readout-arm1",
        "gemma4-26b-q4-base-generate-arm1", "mistral-small3.2-24b-generate-coi", "kev-9b",
        "imajev-4b-cpu", "apus-4b-high", "apus-9b-high", "lev-4b-gpu", "kev-4b", "opendecider-small",
        "opendecider-base-qwen3-4b", "apus-9b-low", "apus-4b-low", "opendecider-nano", "deem-0.8-d0"]
ps = sorted((d["mcnemar"][f"{m}|naive-bayes"]["p"], m) for m in rows)
m_ = len(ps)
for i, (p, m) in enumerate(ps):
    holm = p * (m_ - i)
    print(f"  {m:40s} p {p:.2e}  Holm-adjusted {min(1,holm):.2e}  Bonferroni {min(1,p*m_):.2e}")

print("\n== 4. Held-out gap sign: s0 minus h48 ==")
for m, v in d["s0_minus_h48"].items():
    lo, hi = v["newcombe"]
    print(f"  {m:30s} s0 {v['s0']:7s} h48 {v['h48']:6s} diff {v['diff_pts']:+6.1f}  [{lo:+6.1f}, {hi:+6.1f}] width {hi-lo:5.1f}")

print("\n== 5. Door intervals ==")
for label, k, n in [("APUS 9B high planted", 32, 36), ("APUS 9B high harmless", 0, 12),
                    ("live door planted", 27, 36), ("gate floor planted", 30, 36),
                    ("field exam 40/40", 40, 40), ("Lev long pages 52/53", 52, 53),
                    ("OpenJev held-out 44/48", 44, 48)]:
    w = wilson(k, n)
    cp = beta_inc_cp(k, n)
    print(f"  {label:26s} {k}/{n}: Wilson {100*w[0]:.1f}-{100*w[1]:.1f}; Clopper-Pearson {100*cp[0]:.1f}-{100*cp[1]:.1f}")

print("\n  Live door v APUS 9B high, planted: best case for a difference (all 24 named: live 24, APUS 23;")
print("  unnamed: live 3, APUS 9 nested) -> discordant 6 v 1, exact p", round(binom_two_sided_exact(6, 1), 4))
print("  worst case (unnamed disjoint: live 3 not in APUS 9) -> 9 v 4, p",
      round(binom_two_sided_exact(9, 4), 4))

print("\n== 6. Writer split: Gemma 4 (arm 1 generate) v Kev-9B, interaction ==")
g = d["post_hoc"]["author_family"]["paired"]
gm = g["gemma|gemma4-26b-q4-base-generate-arm1|kev-9b"]
mm = g["mistral|gemma4-26b-q4-base-generate-arm1|kev-9b"]
print("  gemma lines", gm, " mistral lines", mm)


def paired_var(b, c, n):
    return ((b + c) / n - ((b - c) / n) ** 2) / n


d1 = (gm["only_first_right"] - gm["only_second_right"]) / 51
d2 = (mm["only_first_right"] - mm["only_second_right"]) / 57
se = sqrt(paired_var(gm["only_first_right"], gm["only_second_right"], 51) +
          paired_var(mm["only_first_right"], mm["only_second_right"], 57))
z = (d1 - d2) / se
p = math.erfc(abs(z) / sqrt(2))
print(f"  lead on gemma lines {100*d1:.1f} pts, on mistral lines {100*d2:.1f} pts; "
      f"difference {100*(d1-d2):.1f} pts, Wald z {z:.2f}, p {p:.3f}")
# Fisher-type conditional check on discordant direction
def fisher_2x2(a, b, c, dd):
    n = a + b + c + dd
    r1, c1 = a + b, a + c
    def pr(x):
        return comb(r1, x) * comb(n - r1, c1 - x) / comb(n, c1)
    p0 = pr(a)
    lo_, hi_ = max(0, c1 - (n - r1)), min(r1, c1)
    return sum(pr(x) for x in range(lo_, hi_ + 1) if pr(x) <= p0 * (1 + 1e-9))
print("  Fisher on discordant direction [[13,1],[6,2]]: p",
      round(fisher_2x2(gm["only_first_right"], gm["only_second_right"],
                       mm["only_first_right"], mm["only_second_right"]), 3))

print("\n  Gemma 4 on the 57 Mistral lines v naive Bayes:", g["mistral|gemma4-26b-q4-base-generate-arm1|naive-bayes"])
print("  Kev-9B on the 57 Mistral lines v naive Bayes:", g["mistral|kev-9b|naive-bayes"])
print("  five strongest v NB by writer:")
for m in ["kev-9b", "imajev-4b-cpu", "apus-4b-high", "apus-9b-high", "lev-4b-gpu"]:
    gg = g[f"gemma|{m}|naive-bayes"]; mmm = g[f"mistral|{m}|naive-bayes"]
    print(f"   {m:16s} gemma {gg['only_first_right']}-{gg['only_second_right']} p {gg['p']:.3f} | "
          f"mistral {mmm['only_first_right']}-{mmm['only_second_right']} p {mmm['p']:.3f}")
print("  number of post-hoc paired tests in the writer split:", len(g))

print("\n== 7. Option-order spread within a model v between models ==")
for m, v in d["reorder"].items():
    r = v["right_by_order"]
    print(f"  {m:28s} {r}  range {max(r)-min(r)}  mean {sum(r)/len(r):.1f}")

print("\n== 8. Held-out writer mix effect on expected share ==")
# kit: 51 gemma / 57 mistral; h48: 20 gemma / 28 mistral
for m in ["kev-9b", "apus-9b-high", "naive-bayes", "openjev-q4km-readout-arm13"]:
    cts = d["post_hoc"]["author_family"]["counts"][m]
    pg, pm = cts["gemma"] / 51, cts["mistral"] / 57
    kit = (51 * pg + 57 * pm) / 108
    h = (20 * pg + 28 * pm) / 48
    print(f"  {m:28s} gemma {100*pg:.1f}% mistral {100*pm:.1f}% -> expected kit {100*kit:.1f}, "
          f"expected h48 mix {100*h:.1f}, shift {100*(kit-h):+.1f} pts")