#!/usr/bin/env python3
"""Standalone, dependency-free re-verification of a Review Console
verification bundle (GET /batches/:id/verification-bundle).

The pitch: a buyer's own engineer runs THIS script, on THEIR machine, against
the JSON bundle we handed them, and it recomputes every number in the
bundle's `metrics` block from the raw `judgments`/`defect_records` it also
contains — percent agreement, Cohen's kappa, Krippendorff's alpha, Gwet's
AC1, gold accuracy, defect rate, rubric pass rate, and every bootstrap /
Wilson confidence interval — matching our numbers to within 1e-9, or telling
them exactly where it doesn't. No import beyond the Python standard library
(hashlib, math, json, struct, sys) — nothing to install, nothing to trust
but arithmetic anyone can read.

Every algorithm here is a straight port of:
  server/src/lib/agreement.ts    (point-estimate statistics)
  server/src/lib/statistics.ts   (bootstrap CI + Wilson score interval)
  server/src/lib/sampling.ts     (the seeded mulberry32 RNG bootstrapCI reuses)
The bundle's own `recompute` section documents the same algorithms in prose;
this script is the executable version of that documentation.

Usage:
    python3 verify-bundle.py <bundle.json>
    cat bundle.json | python3 verify-bundle.py -

Prints a JSON report to stdout ({"ok": bool, "checks": [...]}) and exits 0
if every check passed, 1 otherwise.
"""
from __future__ import annotations

import hashlib
import json
import math
import struct
import sys
from typing import Callable, Optional

TOLERANCE = 1e-9
DEFAULT_RESAMPLES = 2000
DEFAULT_LEVEL = 0.95
Z_95 = 1.959963984540054  # standard-normal 97.5th percentile — the only confidence level this build emits (statistics.ts's zScore table)

# ---------------------------------------------------------------------------
# Seeded RNG — mulberry32, seeded from the first 4 bytes (little-endian) of
# SHA-256(seed_string). Every arithmetic step below is kept in unsigned
# 32-bit space (mask with 0xFFFFFFFF) so it matches JavaScript's bitwise
# operator semantics bit-for-bit: Math.imul is an exact multiply reduced to
# its low 32 bits; `>>>` is *always* a logical (zero-fill) shift regardless
# of sign, which is naturally what Python's `>>` does on a non-negative int.
# ---------------------------------------------------------------------------
MASK32 = 0xFFFFFFFF


def _u32(x: int) -> int:
    return x & MASK32


def _imul32(a: int, b: int) -> int:
    """Math.imul(a, b): exact multiply, keep the low 32 bits. Equivalent mod 2^32 whether operands are treated as signed or unsigned (see script header note in docs/API.md's Verification Bundle section for the proof sketch)."""
    return _u32(a * b)


def make_rng(seed: str) -> Callable[[], float]:
    digest = hashlib.sha256(seed.encode("utf-8")).digest()
    state = struct.unpack("<I", digest[:4])[0]  # readUInt32LE equivalent

    box = {"state": state}

    def rng() -> float:
        box["state"] = _u32(box["state"] + 0x6D2B79F5)
        t = box["state"]
        t = _imul32(t ^ _u32(t >> 15), t | 1)
        t = t ^ _u32(t + _imul32(t ^ _u32(t >> 7), t | 61))
        return _u32(t ^ _u32(t >> 14)) / 4294967296.0

    return rng


# ---------------------------------------------------------------------------
# agreement.ts port — point-estimate statistics. Dict/list insertion order is
# kept identical to the TS Map-building order throughout, since floating-
# point summation order can (in principle) affect the last bit or two.
# ---------------------------------------------------------------------------


def percent_agreement(units: list[list[str]]) -> Optional[float]:
    ratable = [u for u in units if len(u) > 0]
    if not ratable:
        return None
    unanimous = sum(1 for u in ratable if all(r == u[0] for r in u))
    return unanimous / len(ratable)


def cohens_kappa(a: list[str], b: list[str]) -> Optional[float]:
    if len(a) != len(b):
        raise ValueError("cohens_kappa: a and b must be the same length.")
    n = len(a)
    if n == 0:
        return None
    agree = 0
    margin_a: dict[str, int] = {}
    margin_b: dict[str, int] = {}
    for av, bv in zip(a, b):
        if av == bv:
            agree += 1
        margin_a[av] = margin_a.get(av, 0) + 1
        margin_b[bv] = margin_b.get(bv, 0) + 1
    po = agree / n
    categories = list(margin_a.keys()) + [c for c in margin_b.keys() if c not in margin_a]
    pe = sum((margin_a.get(c, 0) / n) * (margin_b.get(c, 0) / n) for c in categories)
    if pe == 1:
        return 1.0
    return (po - pe) / (1 - pe)


def krippendorff_alpha(units: list[list[Optional[str]]], metric: str = "nominal") -> Optional[float]:
    """B4: `metric="ordinal"` ports agreement.ts's ordinal branch — see that
    function's doc comment for the Krippendorff (2011) citation."""
    coincidence: dict[str, dict[str, float]] = {}

    def bump(c: str, k: str, w: float) -> None:
        coincidence.setdefault(c, {})
        coincidence[c][k] = coincidence[c].get(k, 0.0) + w

    for unit in units:
        valid = [v for v in unit if v is not None]
        m = len(valid)
        if m < 2:
            continue
        weight = 1.0 / (m - 1)
        for i in range(m):
            for j in range(m):
                if i == j:
                    continue
                bump(valid[i], valid[j], weight)

    categories: list[str] = []
    seen = set()
    for c, row in coincidence.items():
        if c not in seen:
            categories.append(c)
            seen.add(c)
        for k in row.keys():
            if k not in seen:
                categories.append(k)
                seen.add(k)
    if not categories:
        return None

    n = 0.0
    marginals: dict[str, float] = {}
    for c in categories:
        s = sum(coincidence.get(c, {}).get(k, 0.0) for k in categories)
        marginals[c] = s
        n += s
    if n <= 1:
        return None

    if metric == "nominal":
        sum_diag = sum(coincidence.get(c, {}).get(c, 0.0) for c in categories)
        sum_marginals_sq = sum(marginals[c] ** 2 for c in categories)
        denominator = n * n - sum_marginals_sq
        if denominator == 0:
            return 1.0
        return 1 - ((n - 1) * (n - sum_diag)) / denominator

    # Ordinal: rank categories by numeric value ascending; weigh every
    # off-diagonal (c,k) pair by the cumulative-marginal delta2 instead of
    # nominal's flat 1.
    ranked = sorted(categories, key=lambda c: float(c))
    prefix = [0.0]
    for c in ranked:
        prefix.append(prefix[-1] + marginals.get(c, 0.0))

    def cumulative_between(lo_idx: int, hi_idx: int) -> float:
        return prefix[hi_idx + 1] - prefix[lo_idx]

    do_sum = 0.0
    de_sum = 0.0
    for i, ci in enumerate(ranked):
        for j, cj in enumerate(ranked):
            if i == j:
                continue
            lo, hi = min(i, j), max(i, j)
            nc = marginals.get(ci, 0.0)
            nk = marginals.get(cj, 0.0)
            d = cumulative_between(lo, hi) - (nc + nk) / 2
            delta2 = d * d
            do_sum += coincidence.get(ci, {}).get(cj, 0.0) * delta2
            de_sum += nc * nk * delta2
    if de_sum == 0:
        return 1.0
    return 1 - ((n - 1) * do_sum) / de_sum


def weighted_kappa(a: list[float], b: list[float], weights: str, categories: list[float]) -> Optional[float]:
    """Port of agreement.ts's weightedKappa — see that function's doc
    comment for the scikit-learn `cohen_kappa_score` source citation."""
    if len(a) != len(b):
        raise ValueError("weighted_kappa: a and b must be the same length.")
    n = len(a)
    if n == 0:
        return None
    ranked = sorted(set(categories))
    q = len(ranked)
    if q < 2:
        return 1.0
    rank_of = {c: i for i, c in enumerate(ranked)}

    observed = [[0.0] * q for _ in range(q)]
    row_sum = [0.0] * q
    col_sum = [0.0] * q
    for av, bv in zip(a, b):
        ri, rj = rank_of[av], rank_of[bv]
        observed[ri][rj] += 1
        row_sum[ri] += 1
        col_sum[rj] += 1

    numerator = 0.0
    denominator = 0.0
    for i in range(q):
        for j in range(q):
            w = abs(i - j) if weights == "linear" else (i - j) ** 2
            if w == 0:
                continue
            expected = row_sum[i] * col_sum[j] / n
            numerator += w * observed[i][j]
            denominator += w * expected
    if denominator == 0:
        return 1.0
    return 1 - numerator / denominator


def gwet_ac2(a: list[float], b: list[float], weights: str, categories: list[float]) -> Optional[float]:
    """Port of agreement.ts's gwetAC2 — see that function's doc comment for
    the irrCAC (Gwet's own reference implementation) source citation."""
    if len(a) != len(b):
        raise ValueError("gwet_ac2: a and b must be the same length.")
    n = len(a)
    if n == 0:
        return None
    ranked = sorted(set(categories))
    q = len(ranked)
    if q < 2:
        return 1.0
    xmin, xmax = ranked[0], ranked[-1]
    range2 = (xmax - xmin) ** 2
    rng = xmax - xmin

    def w(ci: float, cj: float) -> float:
        return 1 - (ci - cj) ** 2 / range2 if weights == "quadratic" else 1 - abs(ci - cj) / rng

    sum_w = sum(w(ci, cj) for ci in ranked for cj in ranked)

    pa_sum = 0.0
    count_a: dict[float, int] = {}
    count_b: dict[float, int] = {}
    for av, bv in zip(a, b):
        pa_sum += w(av, bv)
        count_a[av] = count_a.get(av, 0) + 1
        count_b[bv] = count_b.get(bv, 0) + 1
    pa = pa_sum / n

    pe_sum = 0.0
    for c in ranked:
        pi = (count_a.get(c, 0) + count_b.get(c, 0)) / (2 * n)
        pe_sum += pi * (1 - pi)
    pe = (sum_w * pe_sum) / (q * (q - 1))

    if 1 - pe == 0:
        return 1.0
    return (pa - pe) / (1 - pe)


def gwet_ac1(a: list[str], b: list[str]) -> Optional[float]:
    if len(a) != len(b):
        raise ValueError("gwet_ac1: a and b must be the same length.")
    n = len(a)
    if n == 0:
        return None
    agree = 0
    pooled: dict[str, int] = {}
    for av, bv in zip(a, b):
        if av == bv:
            agree += 1
        pooled[av] = pooled.get(av, 0) + 1
        pooled[bv] = pooled.get(bv, 0) + 1
    po = agree / n
    categories = list(pooled.keys())
    q = len(categories)
    if q < 2:
        return 1.0
    pe_sum = 0.0
    for c in categories:
        pi = pooled[c] / (2 * n)
        pe_sum += pi * (1 - pi)
    pe_ac1 = pe_sum / (q - 1)
    return (po - pe_ac1) / (1 - pe_ac1)


def gold_accuracy(results: list[Optional[bool]]) -> Optional[float]:
    computable = [r for r in results if r is not None]
    if not computable:
        return None
    return (sum(1 for r in computable if r) / len(computable)) * 100


def defect_rate(item_defect_severities: list[list[str]]) -> Optional[float]:
    if not item_defect_severities:
        return None
    defective = sum(1 for sevs in item_defect_severities if any(s in ("MAJOR", "CRITICAL") for s in sevs))
    return (defective / len(item_defect_severities)) * 100


def rubric_pass_rate(overall_scores: list[float], threshold: float) -> Optional[float]:
    if not overall_scores:
        return None
    passed = sum(1 for s in overall_scores if s >= threshold)
    return (passed / len(overall_scores)) * 100


# ---------------------------------------------------------------------------
# validatorAgreement.ts port — C1_C2_TRUST_BRIEF.md §2.2 (V-P0-5). Validator-
# validator reliability, computed from `bundle["qa_responses"]` — one entry
# per (item, pre-adjudication QA response) pair, already excluding the
# adjudication response itself (see lib/validatorAgreement.ts's module doc
# comment for why). `ci_seed` is the EXACT seed recorded on
# `metrics["validator_agreement"]["interval"]["seed"]` — read back and reused
# verbatim, same discipline as every other interval in this script, never
# reconstructed from inputs_hash.
# ---------------------------------------------------------------------------


def compute_validator_agreement(qa_responses: list[dict], ci_seed: str) -> Optional[dict]:
    by_item: dict[str, list[dict]] = {}
    for r in qa_responses:
        by_item.setdefault(r["item_id"], []).append(r)

    qualifying: list[list[dict]] = []
    for responses in by_item.values():
        usable = [r for r in responses if r["label"] is not None]
        if len(usable) >= 2:
            qualifying.append(usable)
    if not qualifying:
        return None

    distinct_responders = {r["responder_ref"] for responses in qualifying for r in responses}
    is_fixed_pair = len(distinct_responders) == 2 and all(
        len(responses) == 2 and len({r["responder_ref"] for r in responses}) == 2 for responses in qualifying
    )

    if is_fixed_pair:
        id_a, id_b = sorted(distinct_responders)
        a: list[str] = []
        b: list[str] = []
        for responses in qualifying:
            by_responder = {r["responder_ref"]: r["label"] for r in responses}
            a.append(by_responder[id_a])
            b.append(by_responder[id_b])
        value = cohens_kappa(a, b)
        interval = bootstrap_ci(len(a), value, lambda idx: cohens_kappa([a[i] for i in idx], [b[i] for i in idx]), ci_seed)
        return {"value": value, "method": "COHENS_KAPPA", "n_items": len(qualifying), "n_validators": 2, "interval": interval}

    units = [[r["label"] for r in responses] for responses in qualifying]
    value = krippendorff_alpha(units)
    interval = bootstrap_ci(len(units), value, lambda idx: krippendorff_alpha([units[i] for i in idx]), ci_seed)
    return {"value": value, "method": "KRIPPENDORFF_ALPHA", "n_items": len(qualifying), "n_validators": len(distinct_responders), "interval": interval}


# ---------------------------------------------------------------------------
# statistics.ts port — bootstrap percentile CI + Wilson score interval.
# ---------------------------------------------------------------------------


def _percentile(sorted_values: list[float], p: float) -> float:
    idx = p * (len(sorted_values) - 1)
    lo = math.floor(idx)
    hi = math.ceil(idx)
    if lo == hi:
        return sorted_values[lo]
    frac = idx - lo
    return sorted_values[lo] * (1 - frac) + sorted_values[hi] * frac


def bootstrap_ci(
    n: int,
    point: Optional[float],
    stat_on_resample: Callable[[list[int]], Optional[float]],
    seed: str,
    resamples: int = DEFAULT_RESAMPLES,
    level: float = DEFAULT_LEVEL,
) -> Optional[dict]:
    if point is None or n < 2:
        return None
    rng = make_rng(seed)
    values: list[float] = []
    for _ in range(resamples):
        indices = [int(rng() * n) for _ in range(n)]  # int() truncates toward 0; rng()*n is always >= 0, so this equals Math.floor
        v = stat_on_resample(indices)
        if v is not None:
            values.append(v)
    if not values:
        return None
    values.sort()
    alpha = 1 - level
    return {"point": point, "lo": _percentile(values, alpha / 2), "hi": _percentile(values, 1 - alpha / 2)}


def sample_random(ids: list[str], n: int, seed: str) -> list[str]:
    """Port of sampling.ts's sampleRandom — partial Fisher-Yates over a seeded mulberry32 stream."""
    rng = make_rng(seed)
    pool = list(ids)
    take = max(0, min(n, len(pool)))
    for i in range(take):
        j = i + int(rng() * (len(pool) - i))  # int() truncates toward 0; rng()*(...) is always >= 0, so this equals Math.floor
        pool[i], pool[j] = pool[j], pool[i]
    return pool[:take]


def allocate_with_caps(entries: list[dict], total: int) -> dict[str, int]:
    """Port of sampling.ts's allocateWithCaps: deterministic weighted largest-remainder (Hamilton) allocation with per-entry caps."""
    fixed: dict[str, int] = {}
    pool = list(entries)
    remaining = total
    for _ in range(len(entries) + 1):
        if not pool:
            break
        total_weight = sum(e["weight"] for e in pool)
        exact = [{"key": e["key"], "exact": (remaining * e["weight"] / total_weight) if total_weight > 0 else 0.0} for e in pool]
        alloc = {e["key"]: math.floor(e["exact"]) for e in exact}
        leftover = remaining - sum(alloc.values())
        by_remainder_desc = sorted(exact, key=lambda e: (-(e["exact"] - math.floor(e["exact"])), e["key"]))
        for e in by_remainder_desc:
            if leftover <= 0:
                break
            alloc[e["key"]] = alloc.get(e["key"], 0) + 1
            leftover -= 1
        cap_by_key = {e["key"]: e["cap"] for e in pool}
        over_cap = [e for e in pool if alloc.get(e["key"], 0) > cap_by_key[e["key"]]]
        if not over_cap:
            for e in pool:
                fixed[e["key"]] = alloc.get(e["key"], 0)
            return fixed
        for e in over_cap:
            fixed[e["key"]] = cap_by_key[e["key"]]
            remaining -= cap_by_key[e["key"]]
        over_cap_keys = {e["key"] for e in over_cap}
        pool = [e for e in pool if e["key"] not in over_cap_keys]
    return fixed


def replay_risk_weighted_draw(
    population_item_ids: list[str],
    contributor_keys: dict[str, Optional[str]],
    n: int,
    seed: str,
    params: dict,
    recorded_strata: list[dict],
) -> dict:
    """Port of sampling.ts's sampleRiskWeightedByContributor (allocation + draw), given the
    bundle's own recorded inputs. `contributor_keys` covers EVERY population item (sampled or
    not — `bundle["sampling"]["population_contributor_keys"]`), so N_k (stratum population size)
    is independently recomputed here, not trusted from the record. `items_completed`/
    `defects_involved` per stratum (`recorded_strata`, keyed by `key`) come from
    contributor_stats — a live DB table this script has no access to — so smoothed_rate/
    pooled_rate/multiplier are recomputed FROM those recorded counts (a consistency check on
    the arithmetic) rather than from scratch; the allocation and the per-stratum draw are then
    fully independent replays. Returns {sampled_ids: set[str], strata: {key: {...recomputed}}}."""
    a, b, floor_pct = params["a"], params["b"], params["floor_pct"]
    clamp_lo, clamp_hi = params["clamp"]
    unknown_multiplier = params["unknown_multiplier"]

    groups: dict[str, list[str]] = {}
    for item_id in population_item_ids:
        key = contributor_keys.get(item_id) or "(none)"
        groups.setdefault(key, []).append(item_id)
    keys = sorted(groups.keys())
    target = max(0, min(n, len(population_item_ids)))

    recorded_by_key = {s["key"]: s for s in recorded_strata}
    sum_known_d = 0
    sum_known_m = 0
    for key in keys:
        rec = recorded_by_key.get(key, {})
        if rec.get("known") and rec.get("items_completed", 0) > 0:
            sum_known_d += rec["defects_involved"]
            sum_known_m += rec["items_completed"]
    pooled_rate = (sum_known_d + a) / (sum_known_m + a + b)

    strata = []
    for key in keys:
        rec = recorded_by_key.get(key, {})
        N_k = len(groups[key])
        known = bool(rec.get("known"))
        M_k = rec.get("items_completed", 0)
        D_k = rec.get("defects_involved", 0)
        smoothed_rate = (D_k + a) / (M_k + a + b)
        multiplier = min(max(smoothed_rate / pooled_rate, clamp_lo), clamp_hi) if known else unknown_multiplier
        min_k = max(1, math.ceil(floor_pct * N_k))
        strata.append({"key": key, "N_k": N_k, "smoothed_rate": smoothed_rate, "multiplier": multiplier, "min_k": min_k, "cap": N_k - min_k})

    total_min = sum(s["min_k"] for s in strata)
    if target < total_min:
        raise ValueError(f"floor_infeasible: requested {target} below the sum of per-stratum floors ({total_min})")
    residual = target - total_min
    alloc_by_key = allocate_with_caps([{"key": s["key"], "weight": s["N_k"] * s["multiplier"], "cap": s["cap"]} for s in strata], residual)

    sampled_ids: set[str] = set()
    result_strata: dict[str, dict] = {}
    for s in strata:
        n_k = s["min_k"] + alloc_by_key.get(s["key"], 0)
        drawn = sample_random(groups[s["key"]], n_k, f"{seed}:{s['key']}")
        sampled_ids.update(drawn)
        result_strata[s["key"]] = {"population_size": s["N_k"], "sampled_size": len(drawn), "smoothed_rate": s["smoothed_rate"], "multiplier": s["multiplier"]}
    return {"sampled_ids": sampled_ids, "strata": result_strata, "pooled_rate": pooled_rate}


def wilson_interval(successes: int, n: int, level: float = DEFAULT_LEVEL) -> Optional[dict]:
    if n == 0:
        return None
    if level != DEFAULT_LEVEL:
        raise ValueError(f"wilson_interval: unsupported confidence level {level}")
    z = Z_95
    p = successes / n
    z2 = z * z
    denom = 1 + z2 / n
    center = (p + z2 / (2 * n)) / denom
    margin = (z * math.sqrt(p * (1 - p) / n + z2 / (4 * n * n))) / denom
    return {"point": p, "lo": max(0.0, center - margin), "hi": min(1.0, center + margin)}


def _scale_pct(interval: Optional[dict]) -> Optional[dict]:
    if interval is None:
        return None
    return {"point": interval["point"] * 100, "lo": interval["lo"] * 100, "hi": interval["hi"] * 100}


# ---------------------------------------------------------------------------
# §E3 A2 — htEstimator.ts port: Horvitz-Thompson / post-stratified defect
# rate. Within-stratum SRS makes inclusion probability constant per stratum,
# reducing HT to sum_k (N_k/N_eff) * r_k with N_eff = sum of N_k over strata
# that have at least one judged sampled item.
# ---------------------------------------------------------------------------


def horvitz_thompson_defect_rate(strata: list[dict]) -> dict:
    """`strata`: [{"key", "population_size", "defect_flags": [0|1, ...]}]. Returns
    {"defect_rate_pct", "n_eff", "strata": [...], "dropped_strata": [...]}."""
    results = []
    dropped: list[str] = []
    n_eff = 0
    weighted_sum = 0.0
    for s in strata:
        n_judged = len(s["defect_flags"])
        defects = sum(s["defect_flags"])
        rate = (defects / n_judged) if n_judged > 0 else None
        results.append({"key": s["key"], "population_size": s["population_size"], "n_judged": n_judged, "defects": defects, "defect_rate_pct": (rate * 100) if rate is not None else None})
        if n_judged == 0:
            dropped.append(s["key"])
            continue
        n_eff += s["population_size"]
        weighted_sum += s["population_size"] * rate
    return {
        "defect_rate_pct": (weighted_sum / n_eff * 100) if n_eff > 0 else None,
        "n_eff": n_eff,
        "strata": results,
        "dropped_strata": dropped,
    }


def stratified_bootstrap_ci(strata: list[dict], point: Optional[float], seed: str, resamples: int = DEFAULT_RESAMPLES, level: float = DEFAULT_LEVEL) -> Optional[dict]:
    """Port of statistics.ts's stratifiedBootstrapCI. `strata`: [{"key", "weight", "values": [0|1,...]}]."""
    if point is None:
        return None
    sorted_strata = sorted(strata, key=lambda s: s["key"])
    degenerate = sum(1 for s in sorted_strata if len(s["values"]) == 1)
    if all(len(s["values"]) < 2 for s in sorted_strata):
        return None

    rngs = {s["key"]: make_rng(f"{seed}:stratum:{s['key']}") for s in sorted_strata}
    draws: list[float] = []
    for _ in range(resamples):
        stat = 0.0
        for s in sorted_strata:
            n = len(s["values"])
            if n == 0:
                continue
            rng = rngs[s["key"]]
            total = sum(s["values"][int(rng() * n)] for _ in range(n))
            stat += s["weight"] * (total / n)
        draws.append(stat)
    draws.sort()

    alpha = 1 - level
    return {
        "point": point,
        "lo": _percentile(draws, alpha / 2),
        "hi": _percentile(draws, 1 - alpha / 2),
        "method": "STRATIFIED_BOOTSTRAP_PERCENTILE",
        "level": level,
        "resamples": resamples,
        "seed": seed,
        "degenerate_strata": degenerate,
    }


# ---------------------------------------------------------------------------
# Bundle-driven recomputation.
# ---------------------------------------------------------------------------


def _sub(d: dict, key: str) -> dict:
    """d.get(key, {}) has a classic gotcha: if `key` is PRESENT with value None (as every unset interval in this bundle is — json.dumps(None) round-trips to a real `null`, not an absent key), .get's default is never used and the None passes through, crashing the next .get() in a chain. This normalizes both "absent" and "present but null" to {}."""
    return d.get(key) or {}


def load_bundle(path: str) -> dict:
    text = sys.stdin.read() if path == "-" else open(path, "r", encoding="utf-8").read()
    return json.loads(text)


def recompute(bundle: dict) -> list[dict]:
    checks: list[dict] = []

    def check(name: str, expected, actual, tolerance: float = TOLERANCE) -> None:
        if expected is None or actual is None:
            ok = expected == actual
            diff = None
        else:
            diff = abs(expected - actual)
            ok = diff <= tolerance
        checks.append({"name": name, "expected": expected, "actual": actual, "diff": diff, "ok": ok})

    def check_interval(name: str, expected: Optional[dict], actual: Optional[dict], tolerance: float = TOLERANCE) -> None:
        if expected is None or actual is None:
            check(f"{name}.presence", expected is None, actual is None, 0)
            return
        check(f"{name}.point", expected.get("point"), actual["point"], tolerance)
        check(f"{name}.lo", expected.get("lo"), actual["lo"], tolerance)
        check(f"{name}.hi", expected.get("hi"), actual["hi"], tolerance)

    def check_eq(name: str, expected, actual) -> None:
        """Exact-equality check for non-numeric values (sets, lists, strings) — check()'s abs(expected-actual) would crash on these."""
        checks.append({"name": name, "expected": expected, "actual": actual, "diff": None, "ok": expected == actual})

    judgments = bundle["judgments"]
    defect_records = bundle["defect_records"]
    metrics = bundle["metrics"]
    intervals = metrics.get("intervals") or {}
    rubric = bundle.get("rubric")

    # --- reconstruct metricsEngine.ts's per-statistic input arrays, in the
    # exact order `judgments` was shipped (which mirrors work_items.position,
    # per lib/metricsEngine.ts's loadJudgedItemContext ORDER BY) ----------
    paired_units = [(j["original_label"], j["final_label"]) for j in judgments if j["kind"] != "GOLD" and j["final_label"] is not None]
    complete_pairs = [(o, f) for (o, f) in paired_units if o is not None]
    gold_results = [j["gold_correct"] for j in judgments if j["kind"] == "GOLD"]
    defects_by_item: dict[str, list[str]] = {}
    for d in defect_records:
        defects_by_item.setdefault(d["item_id"], []).append(d["severity"])
    item_defect_severities = [defects_by_item.get(j["item_id"], []) for j in judgments]
    rubric_overall_scores = [j["rubric_overall_score"] for j in judgments if j["kind"] != "GOLD" and j["rubric_overall_score"] is not None]

    originals = [o for (o, _f) in paired_units]
    finals = [f for (_o, f) in paired_units]
    complete_originals = [o for (o, _f) in complete_pairs]
    complete_finals = [f for (_o, f) in complete_pairs]

    # B4 -- ordinal-labelKind batches report gwet_ac2 instead of gwet_ac1
    # (mutually exclusive, lib/metricsEngine.ts) -- that key's presence is
    # the bundle's own signal for "this batch's item type is ordinal",
    # independent of which agreement_method the pairing actually selected
    # (a batch with zero paired units still falls back to PERCENT_AGREEMENT
    # regardless of labelKind).
    supplementary = _sub(_sub(metrics, "breakdowns"), "supplementary_agreement")
    is_ordinal = "gwet_ac2" in supplementary
    # Only ordinal-labelKind batches carry numeric labels ("4", "5", ...) --
    # a nominal batch's labels ("A", "B", a rubric's "PASS"/"FAIL", a
    # RANKING's candidate id, ...) are never float-convertible and are never
    # read through these two below (every nominal dispatch branch uses the
    # string originals/finals directly).
    if is_ordinal:
        numeric_originals = [float(o) for o in complete_originals]
        numeric_finals = [float(f) for f in complete_finals]
        ordinal_categories = sorted(set(numeric_originals) | set(numeric_finals))
    else:
        numeric_originals = []
        numeric_finals = []
        ordinal_categories = []

    # --- agreement: whichever method the console selected -----------------
    method = metrics.get("agreement_method")
    if method == "PERCENT_AGREEMENT":
        agreement_value = percent_agreement([[o, f] if o is not None else [f] for (o, f) in paired_units])
    elif method == "COHENS_KAPPA":
        agreement_value = cohens_kappa(originals, finals)
    elif method == "KRIPPENDORFF_ALPHA":
        agreement_value = krippendorff_alpha([[o, f] for (o, f) in paired_units])
    elif method == "WEIGHTED_KAPPA_QUADRATIC":
        agreement_value = weighted_kappa(numeric_originals, numeric_finals, "quadratic", ordinal_categories)
    elif method == "KRIPPENDORFF_ALPHA_ORDINAL":
        agreement_value = krippendorff_alpha([[o, f] for (o, f) in paired_units], metric="ordinal")
    else:
        agreement_value = None
    check("agreement_value", metrics.get("agreement_value"), agreement_value)

    agreement_interval = _sub(intervals, "agreement")
    if method == "COHENS_KAPPA" and paired_units:
        ci = bootstrap_ci(len(originals), agreement_value, lambda idx: cohens_kappa([originals[i] for i in idx], [finals[i] for i in idx]), agreement_interval.get("seed", ""))
    elif method == "KRIPPENDORFF_ALPHA" and paired_units:
        units_list = [[o, f] for (o, f) in paired_units]
        ci = bootstrap_ci(len(units_list), agreement_value, lambda idx: krippendorff_alpha([units_list[i] for i in idx]), agreement_interval.get("seed", ""))
    elif method == "WEIGHTED_KAPPA_QUADRATIC" and paired_units:
        ci = bootstrap_ci(
            len(numeric_originals),
            agreement_value,
            lambda idx: weighted_kappa([numeric_originals[i] for i in idx], [numeric_finals[i] for i in idx], "quadratic", ordinal_categories),
            agreement_interval.get("seed", ""),
        )
    elif method == "KRIPPENDORFF_ALPHA_ORDINAL" and paired_units:
        units_list = [[o, f] for (o, f) in paired_units]
        ci = bootstrap_ci(
            len(units_list),
            agreement_value,
            lambda idx: krippendorff_alpha([units_list[i] for i in idx], metric="ordinal"),
            agreement_interval.get("seed", ""),
        )
    else:
        ci = None
    check_interval("intervals.agreement", intervals.get("agreement"), ci)

    # --- gwet's AC1 (nominal) / AC2 (B4, ordinal) -- supplementary, mutually exclusive ---
    if not is_ordinal:
        ac1 = gwet_ac1(complete_originals, complete_finals) if complete_pairs else None
        check("breakdowns.supplementary_agreement.gwet_ac1", supplementary.get("gwet_ac1"), ac1)
        ac1_ci = None
        if complete_pairs:
            ac1_seed = _sub(intervals, "gwet_ac1").get("seed", "")
            ac1_ci = bootstrap_ci(len(complete_originals), ac1, lambda idx: gwet_ac1([complete_originals[i] for i in idx], [complete_finals[i] for i in idx]), ac1_seed)
        check_interval("intervals.gwet_ac1", intervals.get("gwet_ac1"), ac1_ci)
    else:
        ac2 = gwet_ac2(numeric_originals, numeric_finals, "quadratic", ordinal_categories) if complete_pairs else None
        check("breakdowns.supplementary_agreement.gwet_ac2", supplementary.get("gwet_ac2"), ac2)
        ac2_ci = None
        if complete_pairs:
            ac2_seed = _sub(intervals, "gwet_ac2").get("seed", "")
            ac2_ci = bootstrap_ci(
                len(numeric_originals),
                ac2,
                lambda idx: gwet_ac2([numeric_originals[i] for i in idx], [numeric_finals[i] for i in idx], "quadratic", ordinal_categories),
                ac2_seed,
            )
        check_interval("intervals.gwet_ac2", intervals.get("gwet_ac2"), ac2_ci)

    # --- gold accuracy -------------------------------------------------------
    gold_interval = _sub(intervals, "gold_accuracy")
    gold_acc = gold_accuracy(gold_results)
    check("gold_accuracy_pct", metrics.get("gold_accuracy_pct"), gold_acc)
    gold_ci = bootstrap_ci(len(gold_results), gold_acc, lambda idx: gold_accuracy([gold_results[i] for i in idx]), _sub(gold_interval, "bootstrap").get("seed", ""))
    check_interval("intervals.gold_accuracy.bootstrap", gold_interval.get("bootstrap"), gold_ci)
    computable_gold = [r for r in gold_results if r is not None]
    gold_wilson = _scale_pct(wilson_interval(sum(1 for r in computable_gold if r), len(computable_gold))) if computable_gold else None
    check_interval("intervals.gold_accuracy.wilson", gold_interval.get("wilson"), gold_wilson)

    # --- defect rate ---------------------------------------------------------
    # Naive (unweighted sample) rate + its Wilson interval: computed the same
    # way regardless of sampling method — `.wilson` was NEVER on any weighted
    # estimate, before or after risk weighting existed (§E3 A2 §4.3).
    naive_defect_rate_value = defect_rate(item_defect_severities)
    defective_count = sum(1 for sevs in item_defect_severities if any(s in ("MAJOR", "CRITICAL") for s in sevs))
    defect_wilson = _scale_pct(wilson_interval(defective_count, len(item_defect_severities))) if item_defect_severities else None
    defect_interval = _sub(intervals, "defect_rate")
    check_interval("intervals.defect_rate.wilson", defect_interval.get("wilson"), defect_wilson)

    sampling = bundle.get("sampling")
    is_risk_weighted = bool(sampling) and sampling.get("method") == "RISK_WEIGHTED_BY_CONTRIBUTOR"

    if not is_risk_weighted:
        # Every other method: byte-identical to the pre-A2 path — the
        # canonical column is the naive rate, exactly as always.
        check("defect_rate_pct", metrics.get("defect_rate_pct"), naive_defect_rate_value)
    else:
        # §E3 A2 §4.2/§4.4 — canonical defect_rate_pct becomes the
        # Horvitz-Thompson / post-stratified estimate; the naive rate moves
        # to breakdowns.defect_rate_naive_pct (still independently checked
        # above via the shared .wilson interval, which stays naive).
        flags_by_key: dict[str, list[int]] = {}
        for j in judgments:
            key = j.get("original_contributor_key") or "(none)"
            has_defect = 1 if any(s in ("MAJOR", "CRITICAL") for s in defects_by_item.get(j["item_id"], [])) else 0
            flags_by_key.setdefault(key, []).append(has_defect)
        ht_strata_input = [{"key": s["key"], "population_size": s["populationSize"], "defect_flags": flags_by_key.get(s["key"], [])} for s in sampling["strata"]]
        ht = horvitz_thompson_defect_rate(ht_strata_input)
        check("defect_rate_pct", metrics.get("defect_rate_pct"), ht["defect_rate_pct"])

        breakdowns = _sub(metrics, "breakdowns")
        check("breakdowns.defect_rate_naive_pct", breakdowns.get("defect_rate_naive_pct"), naive_defect_rate_value)
        check_eq("breakdowns.biased_under_risk_weighting", True, breakdowns.get("biased_under_risk_weighting"))
        check_eq("breakdowns.dropped_strata", sorted(ht["dropped_strata"]), sorted(breakdowns.get("dropped_strata") or []))
        recorded_slices = {s["key"]: s.get("defect_rate_pct") for s in _sub(breakdowns, "slices").get("by_contributor") or []}
        for s in ht["strata"]:
            check(f"breakdowns.slices.by_contributor.{s['key']}.defect_rate_pct", recorded_slices.get(s["key"]), s["defect_rate_pct"])

        check_interval("intervals.defect_rate.wilson_naive", defect_interval.get("wilson_naive"), defect_wilson)

        n_eff = ht["n_eff"]
        bootstrap_strata = [
            {"key": s["key"], "weight": (s["population_size"] / n_eff) if n_eff > 0 else 0.0, "values": flags_by_key.get(s["key"], [])}
            for s in ht["strata"]
            if s["n_judged"] > 0
        ]
        recorded_stratified = defect_interval.get("stratified_bootstrap")
        stratified_seed = recorded_stratified.get("seed", "") if recorded_stratified else f"bootstrap:{metrics.get('inputs_hash')}:defect_rate_ht"
        # stratified_bootstrap_ci resamples raw 0|1 values (0-1 scale) — the
        # point estimate must be on the same scale; rescale the result back
        # to *_pct afterward (same convention lib/metricsEngine.ts uses).
        ht_proportion = (ht["defect_rate_pct"] / 100) if ht["defect_rate_pct"] is not None else None
        stratified_raw = stratified_bootstrap_ci(bootstrap_strata, ht_proportion, stratified_seed)
        stratified = _scale_pct(stratified_raw) if stratified_raw else None
        if stratified is not None:
            stratified["method"] = stratified_raw["method"]
            stratified["level"] = stratified_raw["level"]
            stratified["resamples"] = stratified_raw["resamples"]
            stratified["seed"] = stratified_raw["seed"]
            stratified["degenerate_strata"] = stratified_raw["degenerate_strata"]
        check_interval("intervals.defect_rate.stratified_bootstrap", recorded_stratified, stratified)

        # --- replay the draw itself: assert the sampled set matches judgments ---
        replay = replay_risk_weighted_draw(
            sampling["population_item_ids"],
            sampling["population_contributor_keys"],
            sampling["n"],
            sampling["seed"],
            sampling["params"],
            sampling["strata"],
        )
        judged_ids = {j["item_id"] for j in judgments}
        check_eq("risk_sampling.sampled_set", sorted(judged_ids), sorted(replay["sampled_ids"]))
        check("risk_sampling.pooled_rate", sampling.get("pooled_rate"), replay["pooled_rate"])
        for s in sampling["strata"]:
            recomputed = replay["strata"].get(s["key"], {})
            check(f"risk_sampling.strata.{s['key']}.sampledSize", s["sampledSize"], recomputed.get("sampled_size"), 0)
            check(f"risk_sampling.strata.{s['key']}.multiplier", s["multiplier"], recomputed.get("multiplier"))

    # --- rubric pass rate (RUBRIC_SCORED batches only) ------------------------
    if rubric is not None and rubric.get("overall_pass_threshold") is not None:
        rubric_interval = _sub(intervals, "rubric_pass_rate")
        threshold = rubric["overall_pass_threshold"]
        pass_rate = rubric_pass_rate(rubric_overall_scores, threshold)
        check("rubric_pass_rate_pct", metrics.get("rubric_pass_rate_pct"), pass_rate)
        rubric_ci = bootstrap_ci(len(rubric_overall_scores), pass_rate, lambda idx: rubric_pass_rate([rubric_overall_scores[i] for i in idx], threshold), _sub(rubric_interval, "bootstrap").get("seed", ""))
        check_interval("intervals.rubric_pass_rate.bootstrap", rubric_interval.get("bootstrap"), rubric_ci)
        passed_count = sum(1 for s in rubric_overall_scores if s >= threshold)
        rubric_wilson = _scale_pct(wilson_interval(passed_count, len(rubric_overall_scores))) if rubric_overall_scores else None
        check_interval("intervals.rubric_pass_rate.wilson", rubric_interval.get("wilson"), rubric_wilson)

    # --- validator-validator reliability (C1_C2_TRUST_BRIEF.md §2.2, V-P0-5) ---
    # QA P1-3(a): this was previously dead code — compute_validator_agreement
    # was defined but never called from here. Wired in now: recompute from
    # `bundle["qa_responses"]` and diff against the recorded
    # `metrics["validator_agreement"]`.
    recorded_va = metrics.get("validator_agreement")
    va_seed = ((recorded_va or {}).get("interval") or {}).get("seed", "")
    recomputed_va = compute_validator_agreement(bundle.get("qa_responses", []), va_seed)
    if recorded_va is None or recomputed_va is None:
        check_eq("validator_agreement.presence", recorded_va is not None, recomputed_va is not None)
    else:
        check("validator_agreement.value", recorded_va.get("value"), recomputed_va["value"])
        check_eq("validator_agreement.method", recorded_va.get("method"), recomputed_va["method"])
        check_eq("validator_agreement.n_items", recorded_va.get("n_items"), recomputed_va["n_items"])
        check_eq("validator_agreement.n_validators", recorded_va.get("n_validators"), recomputed_va["n_validators"])
        check_interval("validator_agreement.interval", recorded_va.get("interval"), recomputed_va["interval"])

    return checks


def main() -> int:
    if len(sys.argv) != 2:
        print("usage: verify-bundle.py <bundle.json | ->", file=sys.stderr)
        return 2
    bundle = load_bundle(sys.argv[1])
    checks = recompute(bundle)
    ok = all(c["ok"] for c in checks)
    print(json.dumps({"ok": ok, "checks": checks}, indent=2))
    return 0 if ok else 1


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