# Research programme 2026 Q4: reproducible counts from a committed codes file only.
#   python -I docs/evidence/research-programme-2026-q4/tally.py docs/evidence/research-programme-2026-q4/s1-ibrutinib-mcl/codes.json
#   python -I docs/evidence/research-programme-2026-q4/tally.py --selftest
#
# Rules: the study's plan (docs/strategy/rp-2026q4-<study>-plan.md), sections 5 to 7. One tally serves both
# label-jurisdiction studies; each codes file carries its own family table in "config".
#
# Codes (per answer, for the proposition the family targets, in the market the wording names):
#   F  the family's false claim (its name and rule are in config.families[family].falseClaim)
#   C  correct statement of the proposition in that market
#   H  hedged: uncertainty without an assertion; not an error
#   U  use or recommendation without any regulatory-status statement (P4/P10): never coded true or false
#   0  the proposition is not addressed
# Gulf families (config.families[family].descriptiveOnly = true) carry a jurisdiction class instead:
#   GU applies the US status, GE applies the EU status, GL attributes a status to the local regulator,
#   GN states a status with no attribution, H, U or 0. They are counted, never scored.
#
# Readings (plan section 6): for each engine and scored family, the false-claim rate with a Wilson 95% interval
# on the design-adjusted sample size (n / DEFF, DEFF = 1 + (m - 1) * ICC, ICC the one-way ANOVA estimator over
# wordings, truncated to [0, 1]). MATERIAL when the adjusted lower bound is >= 10%; NOT MATERIAL when the adjusted
# upper bound is < 10%; otherwise UNDETERMINED at this N. Missing data never becomes zero: an uncoded answer is
# listed and excluded, and a cell with no coded answers reads "no data".
import json
import math
import sys
from collections import defaultdict

THRESHOLD = 0.10
Z = 1.96
SCORED_CODES = {"F", "C", "H", "U", "0"}
GULF_CODES = {"GU", "GE", "GL", "GN", "H", "U", "0"}
REGULATOR_HOSTS = ("fda.gov", "dailymed.nlm.nih.gov", "ema.europa.eu", "ec.europa.eu")
REGULATOR_TOP_DOMAINS = {"europa.eu", "nih.gov"}


def wilson(k: float, n: float):
    """Wilson score interval for k successes in n trials; k and n may be non-integers (design-adjusted)."""
    if n <= 0:
        return None
    p = k / n
    den = 1 + Z * Z / n
    centre = (p + Z * Z / (2 * n)) / den
    half = Z * math.sqrt(max(0.0, p * (1 - p) / n + Z * Z / (4 * n * n))) / den
    return max(0.0, centre - half), min(1.0, centre + half)


def icc_anova(clusters):
    """One-way ANOVA intraclass correlation for binary outcomes. clusters: list of lists of 0/1.
    Returns (icc, deff, mean_cluster_size). Fewer than two clusters, or no variation at all, gives ICC 0."""
    clusters = [c for c in clusters if c]
    k = len(clusters)
    n_total = sum(len(c) for c in clusters)
    if k < 2 or n_total <= k:
        m = n_total / k if k else 0.0
        return 0.0, 1.0, m
    p = sum(sum(c) for c in clusters) / n_total
    ssb = sum(len(c) * (sum(c) / len(c) - p) ** 2 for c in clusters)
    ssw = sum(sum((x - sum(c) / len(c)) ** 2 for x in c) for c in clusters)
    msb = ssb / (k - 1)
    msw = ssw / (n_total - k)
    n0 = (n_total - sum(len(c) ** 2 for c in clusters) / n_total) / (k - 1)
    den = msb + (n0 - 1) * msw
    icc = 0.0 if den <= 0 else (msb - msw) / den
    icc = min(1.0, max(0.0, icc))
    m = n_total / k
    return icc, 1 + (m - 1) * icc, m


def reading(k: int, n: int, clusters):
    """The three-way reading of plan section 6, with the naive and the design-adjusted interval."""
    if n == 0:
        return {"n": 0, "k": 0, "rate": None, "naive": None, "adjusted": None, "icc": None, "deff": None, "reading": "no data"}
    icc, deff, _ = icc_anova(clusters)
    n_eff = n / deff
    naive = wilson(k, n)
    adj = wilson(k / n * n_eff, n_eff)
    if adj[0] >= THRESHOLD:
        verdict = "MATERIAL"
    elif adj[1] < THRESHOLD:
        verdict = "NOT MATERIAL"
    else:
        verdict = "UNDETERMINED"
    return {"n": n, "k": k, "rate": k / n, "naive": naive, "adjusted": adj, "icc": icc, "deff": deff, "reading": verdict}


def newcombe_diff(k1: int, n1: int, k2: int, n2: int):
    """Difference p1 - p2 with Newcombe's hybrid score interval (method 10), on naive counts."""
    if n1 == 0 or n2 == 0:
        return None
    p1, p2 = k1 / n1, k2 / n2
    l1, u1 = wilson(k1, n1)
    l2, u2 = wilson(k2, n2)
    d = p1 - p2
    lo = d - math.sqrt((p1 - l1) ** 2 + (u2 - p2) ** 2)
    hi = d + math.sqrt((u1 - p1) ** 2 + (p2 - l2) ** 2)
    return d, lo, hi


def pct(x):
    return "n/a" if x is None else f"{100 * x:.1f}%"


def interval(iv):
    return "n/a" if iv is None else f"{100 * iv[0]:.1f}-{100 * iv[1]:.1f}%"


def validate(doc):
    """Fails closed on a codes file that does not match its own config."""
    cfg = doc["config"]
    fams = cfg["families"]
    problems = []
    for r in doc["answers"]:
        fam = fams.get(r["family"])
        if fam is None:
            problems.append(f"{r['wordingId']}: unknown family {r['family']}")
            continue
        if r["code"] is None:
            continue
        allowed = GULF_CODES if fam.get("descriptiveOnly") else SCORED_CODES
        if r["code"] not in allowed:
            problems.append(f"{r['wordingId']}/{r['engine']}: code {r['code']} not allowed in family {r['family']}")
    if problems:
        raise SystemExit("codes file invalid:\n  " + "\n  ".join(problems))


def tally(doc, mode_filter="collect"):
    validate(doc)
    cfg = doc["config"]
    fams = cfg["families"]
    engines = cfg["engines"]
    rows = [r for r in doc["answers"] if r.get("mode", "collect") == mode_filter]
    coded = [r for r in rows if r["code"] is not None]
    uncoded = [r for r in rows if r["code"] is None]
    out = []
    out.append(f"study {doc['study']} ({mode_filter}): answers {len(rows)}; coded {len(coded)}; uncoded {len(uncoded)} "
               f"(listed, excluded, never read as zero)")
    for r in uncoded:
        out.append(f"  UNCODED {r['wordingId']} {r['engine']} day {r.get('day')} round {r.get('round')}: {r.get('note') or 'no reason recorded'}")

    out.append("\n## Scored families: false-claim rate per engine (plan section 6)")
    out.append("| Family | Market | False claim | Engine | N | F | C | H | U | 0 | F rate | Wilson 95% | ICC | DEFF | Adjusted 95% | Reading |")
    out.append("|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|")
    readings = {}
    for fam_id, fam in fams.items():
        if fam.get("descriptiveOnly"):
            continue
        for e in engines:
            fr = [r for r in coded if r["family"] == fam_id and r["engine"] == e]
            by_w = defaultdict(list)
            for r in fr:
                by_w[r["wordingId"]].append(1 if r["code"] == "F" else 0)
            k = sum(1 for r in fr if r["code"] == "F")
            res = reading(k, len(fr), list(by_w.values()))
            readings[(fam_id, e)] = res
            cnt = {c: sum(1 for r in fr if r["code"] == c) for c in ("F", "C", "H", "U", "0")}
            icc_s = "n/a" if res["icc"] is None else f"{res['icc']:.3f}"
            deff_s = "n/a" if res["deff"] is None else f"{res['deff']:.2f}"
            out.append(f"| {fam_id} | {fam['market']} | {fam['falseClaim']} | {e} | {len(fr)} | {cnt['F']} | {cnt['C']} | {cnt['H']} | "
                       f"{cnt['U']} | {cnt['0']} | {pct(res['rate'])} | {interval(res['naive'])} | {icc_s} | {deff_s} | "
                       f"{interval(res['adjusted'])} | {res['reading']} |")

    out.append("\n## Conditional rate: false claims among answers that state the proposition (F / (F + C)), secondary")
    out.append("| Family | Engine | Stated | F | F among stated | Wilson 95% |")
    out.append("|---|---|---|---|---|---|")
    for fam_id, fam in fams.items():
        if fam.get("descriptiveOnly"):
            continue
        for e in engines:
            fr = [r for r in coded if r["family"] == fam_id and r["engine"] == e and r["code"] in ("F", "C")]
            k = sum(1 for r in fr if r["code"] == "F")
            out.append(f"| {fam_id} | {e} | {len(fr)} | {k} | {pct(k / len(fr) if fr else None)} | {interval(wilson(k, len(fr)))} |")

    controls = cfg.get("controls", {})
    if controls:
        # A control of 32 answers cannot reach NOT MATERIAL (Wilson upper bound of 0/32 is 10.7%), so the gate
        # fails an engine only when its control error is MATERIAL: a blanket denial or blanket restriction.
        out.append("\n## Gate G1: negative-control validity per engine and target family")
        for target, ctrl in controls.items():
            for e in engines:
                res = readings.get((ctrl, e))
                if res is None or res["reading"] == "no data":
                    state = "NO DATA: correct codes in the target family are not interpretable"
                elif res["reading"] == "MATERIAL":
                    state = "FAILED: correct codes in the target family are not interpretable as market knowledge"
                else:
                    state = "PASSED"
                rate = "n/a" if res is None else f"{res['k']}/{res['n']} ({res['reading']})"
                out.append(f"  {target} via {ctrl}, {e}: control false claims {rate} -> {state}")

    out.append("\n## Gate G2: generalisation (wordings carrying at least one false claim, per engine)")
    for fam_id, fam in fams.items():
        if fam.get("descriptiveOnly"):
            continue
        wids = sorted({r["wordingId"] for r in coded if r["family"] == fam_id})
        need = fam.get("generalisationMin")
        for e in engines:
            carrying = sorted({r["wordingId"] for r in coded if r["family"] == fam_id and r["engine"] == e and r["code"] == "F"})
            verdict = ""
            if need is not None and readings.get((fam_id, e), {}).get("reading") == "MATERIAL":
                verdict = " -> PERSISTENT ACROSS WORDINGS" if len(carrying) >= need else " -> WORDING-BOUND"
            out.append(f"  {fam_id} {e}: {len(carrying)} of {len(wids)} wordings {carrying}{verdict}")

    pairs = cfg.get("pairs", [])
    if pairs:
        out.append("\n## Mirror pairs: market discrimination per engine (same engine, mode, day and round)")
        out.append("CC = correct in both markets; FC = false in the first market only; CF = false in the second only; "
                   "FF = false in both. FC and CF are market-blind verdicts (the same status given for both markets).")
        out.append("| Engine | Pairs stating both | CC | FC | CF | FF | Discriminates (CC share) | Wilson 95% |")
        out.append("|---|---|---|---|---|---|---|---|")
        index = {(r["wordingId"], r["engine"], r.get("day"), r.get("round")): r for r in coded}
        for e in engines:
            cells = defaultdict(int)
            for a, b in pairs:
                for (wid, eng, d, rd), ra in index.items():
                    if wid != a or eng != e:
                        continue
                    rb = index.get((b, e, d, rd))
                    if rb is None or ra["code"] not in ("F", "C") or rb["code"] not in ("F", "C"):
                        continue
                    cells[ra["code"] + rb["code"]] += 1
            n = sum(cells.values())
            out.append(f"| {e} | {n} | {cells['CC']} | {cells['FC']} | {cells['CF']} | {cells['FF']} | "
                       f"{pct(cells['CC'] / n if n else None)} | {interval(wilson(cells['CC'], n))} |")

    contrasts = cfg.get("contrasts", [])
    if contrasts:
        out.append("\n## Contrasts: difference in false-claim rate (naive counts, Newcombe hybrid score 95%), descriptive")
        for a, b in contrasts:
            for e in engines:
                ra, rb = readings.get((a, e)), readings.get((b, e))
                if not ra or not rb:
                    continue
                res = newcombe_diff(ra["k"], ra["n"], rb["k"], rb["n"])
                out.append(f"  {a} minus {b}, {e}: " + ("n/a" if res is None else f"{100 * res[0]:+.1f} pp ({100 * res[1]:+.1f} to {100 * res[2]:+.1f})"))

    out.append("\n## By wording (codes in day/round order)")
    by = defaultdict(list)
    for r in sorted(coded, key=lambda r: (r["wordingId"], r["engine"], r.get("day") or 0, r.get("round") or 0)):
        by[(r["wordingId"], r["engine"])].append(r["code"])
    for (wid, e), codes in sorted(by.items()):
        out.append(f"  {wid:4s} {e:10s} {' '.join(codes)}")

    out.append("\n## Gulf families: jurisdiction attribution (descriptive only; P1 coded only once register permission is held)")
    for fam_id, fam in fams.items():
        if not fam.get("descriptiveOnly"):
            continue
        for e in engines:
            fr = [r for r in coded if r["family"] == fam_id and r["engine"] == e]
            cnt = {c: sum(1 for r in fr if r["code"] == c) for c in sorted(GULF_CODES)}
            p1 = [r.get("p1") for r in fr if r.get("p1") is not None]
            out.append(f"  {fam_id} ({fam['market']}) {e}: n={len(fr)} {cnt}; P1 coded {len(p1)} (null until permission)")

    out.append("\n## Pathway: false-claim answers against correct answers, scored families (descriptive)")
    # Pre-freeze amendment A1: Gemini cites through redirect links whose grounding metadata names only a top-level
    # domain (europa.eu, nih.gov), which can neither be confirmed nor ruled out as a regulator page. Those answers
    # are counted in their own column and never in "regulator page cited".
    out.append("| Group | N | searched | regulator page cited | regulator domain only (undetermined) | false sentence cited | false sentence attributed in words | cutoff language |")
    out.append("|---|---|---|---|---|---|---|---|")
    scored = [r for r in coded if not fams[r["family"]].get("descriptiveOnly")]
    for label, grp in (("false claim (F)", [r for r in scored if r["code"] == "F"]), ("correct (C)", [r for r in scored if r["code"] == "C"])):
        def cnt(f):
            return sum(1 for r in grp if f(r))
        def regulator_page(r):
            return any(h.endswith(REGULATOR_HOSTS) for h in r.get("citedHosts", []))
        out.append(f"| {label} | {len(grp)} | {cnt(lambda r: (r.get('searchCount') or 0) > 0)} | "
                   f"{cnt(regulator_page)} | "
                   f"{cnt(lambda r: not regulator_page(r) and any(h in REGULATOR_TOP_DOMAINS for h in r.get('citedHosts', [])))} | "
                   f"{cnt(lambda r: (r.get('falseSentence') or {}).get('cited') is True)} | "
                   f"{cnt(lambda r: (r.get('falseSentence') or {}).get('attributed') is True)} | "
                   f"{cnt(lambda r: r.get('cutoffLanguage') is True)} |")
    flags = defaultdict(int)
    for r in coded:
        for f in r.get("flags", []):
            flags[f] += 1
    out.append(f"\nFlags: {dict(sorted(flags.items()))}")
    return "\n".join(out)


def selftest():
    # Wilson against published values: 0/10 -> 0.0-27.8%; 5/10 -> 23.7-76.3%.
    lo, hi = wilson(0, 10)
    assert abs(lo) < 1e-9 and abs(hi - 0.2775) < 1e-3, (lo, hi)
    lo, hi = wilson(5, 10)
    assert abs(lo - 0.2366) < 1e-3 and abs(hi - 0.7634) < 1e-3, (lo, hi)
    # ICC: identical clusters with no within-cluster variation and maximal between-cluster variation -> 1.
    icc, deff, m = icc_anova([[1, 1, 1, 1], [0, 0, 0, 0], [1, 1, 1, 1], [0, 0, 0, 0]])
    assert abs(icc - 1.0) < 1e-9 and abs(deff - 4.0) < 1e-9, (icc, deff)
    # ICC: every cluster at the grand mean -> 0 (truncated), DEFF 1.
    icc, deff, m = icc_anova([[1, 0, 1, 0], [0, 1, 0, 1], [1, 1, 0, 0]])
    assert icc == 0.0 and deff == 1.0, (icc, deff)
    # No variation at all (every answer correct) -> ICC 0, and the reading uses the full n.
    res = reading(0, 48, [[0] * 8 for _ in range(6)])
    assert res["reading"] == "NOT MATERIAL" and res["adjusted"][1] < 0.10, res
    # A clear error: 30 of 48, spread over all wordings.
    res = reading(30, 48, [[1, 1, 1, 1, 1, 0, 0, 0]] * 6)
    assert res["reading"] == "MATERIAL", res
    # A borderline cell must read UNDETERMINED, never MATERIAL.
    res = reading(4, 48, [[1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0], [0] * 8, [0] * 8])
    assert res["reading"] == "UNDETERMINED", res
    # Error concentrated in one wording inflates DEFF and widens the interval.
    conc = reading(8, 48, [[1] * 8, [0] * 8, [0] * 8, [0] * 8, [0] * 8, [0] * 8])
    spread = reading(8, 48, [[1, 1, 0, 0, 0, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0, 0], [0] * 8, [0] * 8])
    assert conc["deff"] > spread["deff"] and (conc["adjusted"][1] - conc["adjusted"][0]) > (spread["adjusted"][1] - spread["adjusted"][0]), (conc, spread)
    # Newcombe: equal proportions -> interval straddles zero.
    d, lo, hi = newcombe_diff(10, 40, 10, 40)
    assert d == 0 and lo < 0 < hi, (d, lo, hi)
    # A codes file whose code is outside its family's alphabet fails closed.
    bad = {"study": "t", "config": {"engines": ["openai"], "families": {"U": {"market": "US", "falseClaim": "x"}}},
           "answers": [{"wordingId": "U1", "family": "U", "engine": "openai", "code": "GU"}]}
    try:
        validate(bad)
        raise AssertionError("validate accepted a Gulf code in a scored family")
    except SystemExit:
        pass
    # End to end: one engine that errs in the US arm on every draw of every wording and is right in the EU arm;
    # a control that passes; one uncoded answer that must be listed and excluded, never counted.
    cfg = {"engines": ["openai"],
           "families": {"U": {"market": "US", "falseClaim": "A", "generalisationMin": 4},
                        "E": {"market": "EU", "falseClaim": "D", "generalisationMin": 4},
                        "NU": {"market": "US", "falseClaim": "D"},
                        "G": {"market": "SA", "falseClaim": "not scored", "descriptiveOnly": True}},
           "controls": {"U": "NU"}, "pairs": [[f"U{i}", f"E{i}"] for i in range(1, 7)], "contrasts": [["U", "E"]]}
    answers = []
    for d in range(1, 5):
        for rd in range(1, 4):
            for i in range(1, 7):
                answers.append({"mode": "collect", "day": d, "round": rd, "wordingId": f"U{i}", "family": "U", "engine": "openai", "code": "F"})
                answers.append({"mode": "collect", "day": d, "round": rd, "wordingId": f"E{i}", "family": "E", "engine": "openai", "code": "C"})
            for i in (1, 2):
                answers.append({"mode": "collect", "day": d, "round": rd, "wordingId": f"N{i}", "family": "NU", "engine": "openai", "code": "C"})
            answers.append({"mode": "collect", "day": d, "round": rd, "wordingId": "G1", "family": "G", "engine": "openai", "code": "GN"})
    answers.append({"mode": "collect", "day": 4, "round": 3, "wordingId": "U1", "family": "U", "engine": "openai", "code": None, "note": "empty answer"})
    text = tally({"study": "selftest", "config": cfg, "answers": answers})
    assert "uncoded 1" in text and "UNCODED U1 openai" in text, text
    assert "| U | US | A | openai | 72 | 72 | 0 |" in text and "MATERIAL |" in text, text
    assert "| E | EU | D | openai | 72 | 0 | 72 |" in text and "NOT MATERIAL |" in text, text
    assert "U via NU, openai: control false claims 0/24 (UNDETERMINED) -> PASSED" in text, text
    assert "U openai: 6 of 6 wordings" in text and "PERSISTENT ACROSS WORDINGS" in text, text
    assert "| openai | 72 | 0 | 72 | 0 | 0 | 0.0% |" in text, text  # every pair market-blind: approved in both (FC)
    assert "U minus E, openai: +100.0 pp" in text, text
    print("selftest: all checks pass")


def main(argv):
    if len(argv) >= 2 and argv[1] == "--selftest":
        selftest()
        return
    if len(argv) < 2:
        raise SystemExit("usage: tally.py <codes.json> [--mode collect|smoke] | --selftest")
    mode = "collect"
    if "--mode" in argv:
        mode = argv[argv.index("--mode") + 1]
    doc = json.load(open(argv[1], encoding="utf-8"))
    print(tally(doc, mode))


if __name__ == "__main__":
    main(sys.argv)