#!/usr/bin/env python3
"""Recompute a published Daishi study's reading of its metric from its two data files.

    python3 recompute_study.py https://daishi.ai/u/<handle>/studies/<study id> [--json]
    python3 recompute_study.py <study page> --trials trials.csv --seats seats.csv

Downloads a published series' or paired experiment's data files (one row per
trial of each arm, one row per seat of each trial) and the figures its page
prints, from the profile's public JSON. From the files it recomputes the
page's reading of the study's own metric: for each comparison of an arm
against the control, the mean paired difference, the interval at the
declared level, the interval that counts each map once, the exact sign-flip
test and the verdict (under the anytime rule it also reads the design's
tuning constant from the study's recipe); for each arm, its mean over its
trials and their 95% interval; for each model, its mean Fitness Index over
its seats in whole trials, whatever the metric, and its interval counting
each trial once. It prints each
beside the page's figure. Where the metric is a seat score (the Fitness
Index, one of its dimensions, rank, survival or map coverage), it also
rebuilds each trial's value from the seats file and checks the two files
against each other. It follows the rules at /docs/benchmark/reproducing,
"Recomputing a study yourself", and needs Python 3.9 or later and nothing
else: no packages, no account, no key.

Taken on trust, and printed at the end of every run: each seat's scores, as
the match's report published them (recompute.py recomputes any one match
from its public record); which trials played whole; which pairs each
comparison read (the files' pair and in_comparison columns are the page's
own account); why a verdict is withheld, which is the page's reason; and a
trial's value on a metric that is not a seat score. Not recomputed:
secondary outcomes, the detectable gap, the debrief counts, the
output-per-turn table and spend.

A full match means the page's figures on its metric, and its per-model
Fitness Index means, follow from the two files by the platform's own
arithmetic, which this script ports. It does not
show that the files are a true record of play, that a score measures what
its name says, or that this analysis suits the question.

Exit status: 0 when every recomputed figure matches the page and the files
agree with each other and with the page's own counts, 1 when one does not,
2 when the study cannot be recomputed here.
"""
import argparse
import csv
import decimal
import io
import json
import math
import re
import sys
import time
import urllib.error
import urllib.parse
import urllib.request

STUDY_RE = re.compile(r"^(https?://[^/?#]+)?/?u/([^/?#]+)/studies/([^/?#]+)/?(?:[?#].*)?$")
TRIAL_NEEDS = ("study", "trial", "arm", "arm_label", "seed", "metric", "value", "pair", "in_comparison", "stopping_rule", "comparison_level")
SEAT_NEEDS = ("trial", "arm", "seat", "model", "focal", "whole_trial", "fitness_index", "rank", "status")

# Where a metric's values can fall, as the platform's reports bound an
# interval on them; a paired difference on a metric bounded on both sides
# lies within plus or minus its span. A metric missing here is read unbounded.
METRIC_RANGE = dict(
    fitness_index=(0, 100), index_survival=(0, 100), index_economy=(0, 100), index_social=(0, 100),
    index_adaptation=(0, 100), index_competitiveness=(0, 100),
    survived=(0, 1), said_evaluated=(0, 1), map_coverage=(0, 1), compliance_rate=(0, 1), injection_resistance=(0, 1), ring_share=(0, 1),
    rank=(1, None), output_tokens_per_turn=(0, None), compliance_latency=(0, None),
    aggression_per_100_rounds=(0, None), trades_per_100_rounds=(0, None), circular_trades_per_100_rounds=(0, None),
    messages_per_100_rounds=(0, None), aggression_between_seats_per_100_rounds=(0, None),
    trades_between_seats_per_100_rounds=(0, None), betrayals_per_100_rounds=(0, None), witness_sensitivity=(None, None),
)
# A trial's value is the mean of its seats' readings on these, so the seats file recomputes it.
SEAT_READING = dict(
    fitness_index="fitness_index", rank="rank", index_survival="survival", index_economy="economy",
    index_social="social", index_adaptation="adaptation", index_competitiveness="competitiveness", map_coverage="map_coverage",
)


class Unavailable(Exception):
    pass


# ---- Arithmetic exactly as the platform's reports do it ----------------------

def keep_trial_value(x):
    # A trial's value as the record keeps it, to three places, and as the
    # stopping rule reads it. Its rounding belongs to the analysis version:
    # like the platform, it rounds the binary value with halves toward
    # +infinity, as JavaScript's Math.round does, not by the print rule below.
    y = x * 1000
    r = math.floor(y)
    if y - r >= 0.5:
        r += 1
    return r / 1000


def keep_figure(x, places):
    # A figure as the record keeps it, by the rule it prints by: half away
    # from zero on its own digits, read to fifteen significant digits.
    # Python's round() rounds halves to even, and rounding the binary value
    # lets float residue pick a tie's side: 30.4 and 34.3 average
    # 32.349999999999994, which the record keeps as 32.4.
    if not math.isfinite(x):
        return x
    d = decimal.Decimal(format(x, ".15g"))
    return float(d.quantize(decimal.Decimal(1).scaleb(-places), rounding=decimal.ROUND_HALF_UP))


def keep_diff(x, places):
    # A difference as the record keeps it: to its places, except that one that
    # is not zero keeps its side of zero, with up to three more places. It
    # rounds as a figure prints, half away from zero on the difference's own
    # digits to fifteen significant digits, so a tie at the kept place goes
    # the same way on either side of zero.
    if not math.isfinite(x):
        return x
    for p in range(places, places + 4):
        k = keep_figure(x, p)
        if k != 0:
            return k
    return 0.0


def keep_mean_diff(diffs, places, linear):
    # A comparison's mean difference as the record keeps it, from the exact
    # decimal its pairs give: on a linear scale each value is kept to three
    # places, so each difference is a whole number of thousandths, and the
    # mean is kept half away from zero from that fraction, not from a float
    # sum whose residue can decide a tie. On the log scale, the float mean.
    if not linear or places < 3 or not all(math.isfinite(d) for d in diffs):
        return keep_diff(mean(diffs), places)
    s = sum(int(round(d * 1000)) for d in diffs)
    n = len(diffs)
    for p in range(places, places + 4):
        q = (2 * abs(s) * 10 ** (p - 3) + n) // (2 * n)
        if q:
            return math.copysign(q / 10 ** p, s)
    return 0.0


def mean(xs):
    # Summed left to right, as the reports sum: from Python 3.12 sum() compensates
    # float sums and can land a last bit away from them.
    s = 0.0
    for x in xs:
        s += x
    return s / len(xs) if xs else math.nan


def stddev(xs):
    if len(xs) < 2:
        return math.nan
    m = mean(xs)
    s = 0.0
    for x in xs:
        s += (x - m) * (x - m)
    return math.sqrt(s / (len(xs) - 1))


def percentile(ordered, p):
    idx = (len(ordered) - 1) * p
    lo, hi = math.floor(idx), math.ceil(idx)
    return ordered[lo] + (ordered[hi] - ordered[lo]) * (idx - lo)


def mulberry32(seed):
    state = [seed & 0xFFFFFFFF]

    def draw():
        a = (state[0] + 0x6D2B79F5) & 0xFFFFFFFF
        state[0] = a
        t = ((a ^ (a >> 15)) * (1 | a)) & 0xFFFFFFFF
        t = ((t + (((t ^ (t >> 7)) * (61 | t)) & 0xFFFFFFFF)) & 0xFFFFFFFF) ^ t
        return ((t ^ (t >> 14)) & 0xFFFFFFFF) / 4294967296

    return draw


def bootstrap(xs, level, resamples=1000, seed=1337):
    """The seeded percentile bootstrap of the mean: the second gate a fixed-rule verdict clears."""
    draw = mulberry32(seed)
    n = len(xs)
    stats = sorted(mean([xs[math.floor(draw() * n)] for _ in range(n)]) for _ in range(resamples))
    tail = (1 - level) / 2
    return percentile(stats, tail), percentile(stats, 1 - tail)


LANCZOS = (0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313,
           -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)


def ln_gamma(z):
    if z < 0.5:
        return math.log(math.pi / math.sin(math.pi * z)) - ln_gamma(1 - z)
    z -= 1
    x = LANCZOS[0]
    for i in range(1, 9):
        x += LANCZOS[i] / (z + i)
    t = z + 7 + 0.5
    return 0.5 * math.log(2 * math.pi) + (z + 0.5) * math.log(t) - t + math.log(x)


def betacf(x, a, b):
    eps, fpmin = 3e-12, 1e-300
    qab, qap, qam = a + b, a + 1, a - 1
    c, d = 1.0, 1 - qab * x / qap
    if abs(d) < fpmin:
        d = fpmin
    d = 1 / d
    h = d
    for m in range(1, 201):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1 + aa * d
        if abs(d) < fpmin:
            d = fpmin
        c = 1 + aa / c
        if abs(c) < fpmin:
            c = fpmin
        d = 1 / d
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1 + aa * d
        if abs(d) < fpmin:
            d = fpmin
        c = 1 + aa / c
        if abs(c) < fpmin:
            c = fpmin
        d = 1 / d
        delta = d * c
        h *= delta
        if abs(delta - 1) < eps:
            break
    return h


def incomplete_beta(x, a, b):
    if x <= 0:
        return 0.0
    if x >= 1:
        return 1.0
    ln_beta = ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b)
    front = math.exp(a * math.log(x) + b * math.log(1 - x) - ln_beta)
    if x < (a + 1) / (a + b + 2):
        return front * betacf(x, a, b) / a
    return 1 - math.exp(b * math.log(1 - x) + a * math.log(x) - ln_beta) * betacf(1 - x, b, a) / b


def t_p_value(t, df):
    """Two-sided p of a t statistic on df degrees of freedom."""
    if not math.isfinite(t) or not math.isfinite(df) or df <= 0:
        return math.nan
    return incomplete_beta(df / (df + t * t), df / 2, 0.5)


def t_critical(df, level):
    """The two-sided Student t critical value, by bisection on the exact p."""
    if df < 1 or not 0 < level < 1:
        return math.inf
    alpha = 1 - level
    lo, hi = 0.0, 1.0
    while t_p_value(hi, df) > alpha:
        hi *= 2
    for _ in range(100):
        mid = (lo + hi) / 2
        if t_p_value(mid, df) > alpha:
            lo = mid
        else:
            hi = mid
    return hi


def bounded(m, half, bounds):
    lo, hi = m - half, m + half
    if bounds[0] is not None:
        lo = max(bounds[0], lo)
    if bounds[1] is not None:
        hi = min(bounds[1], hi)
    return lo, hi


def t_interval(xs, level, bounds):
    if len(xs) < 2:
        return None
    return bounded(mean(xs), t_critical(len(xs) - 1, level) * (stddev(xs) / math.sqrt(len(xs))), bounds)


def clustered_interval(xs, clusters, level, bounds):
    """The mean's interval counting each cluster once: the cluster-robust (CR1) standard error, on clusters minus 1 degrees of freedom."""
    m = mean(xs)
    residual = {}
    for x, c in zip(xs, clusters):
        residual[c] = residual.get(c, 0) + x - m
    g = len(residual)
    if g < 2:
        return None
    sum_sq = 0.0
    for e in residual.values():
        sum_sq += e * e
    se = math.sqrt((g / (g - 1)) * sum_sq) / len(xs)
    return bounded(m, t_critical(g - 1, level) * se, bounds)


def cluster_totals(xs, clusters):
    totals = {}
    for x, c in zip(xs, clusters):
        totals[c] = totals.get(c, 0) + x
    return list(totals.values())


def sign_flip(xs):
    """The exact two-sided sign-flip test of a zero mean; None with no unit or past 20 non-zero units."""
    if not xs:
        return None
    v = [abs(x) for x in xs if x != 0]
    k = len(v)
    if k == 0:
        return 1.0
    if k > 20:
        return None
    total = 0.0
    for a in v:
        total += a
    observed = 0.0
    for x in xs:
        observed += x
    floor = abs(observed) - 1e-9 * total
    sign = [1] * k
    at_least = 1
    # Gray code order: each step flips one sign.
    for i in range(1, 2 ** k):
        bit = (i & -i).bit_length() - 1
        total -= 2 * sign[bit] * v[bit]
        sign[bit] = -sign[bit]
        if abs(total) >= floor:
            at_least += 1
    return at_least / 2 ** k


def sequence_half_width(sd, t, alpha, t_star):
    """The anytime rule's asymptotic confidence sequence after t pairs."""
    lam = -2 * math.log(alpha)
    rho2 = (lam + math.log(lam + 1)) / t_star
    return sd * math.sqrt(((2 * (t * rho2 + 1)) / (t * t * rho2)) * math.log(math.sqrt(t * rho2 + 1) / alpha))


# ---- Reading the study --------------------------------------------------------

def fetch(base, path, accept):
    req = urllib.request.Request(base + path, headers={"User-Agent": "daishi-recompute-study/1", "Accept": accept})
    for attempt in range(5):
        try:
            with urllib.request.urlopen(req, timeout=60) as r:
                return r.read().decode("utf-8-sig")
        except urllib.error.HTTPError as e:
            if e.code in (404, 409):
                said = e.read().decode("utf-8", "replace").strip()
                try:
                    body = json.loads(said)
                except ValueError:
                    body = None
                if isinstance(body, dict):
                    said = body.get("message") or body.get("error") or said
                raise Unavailable(f"{base}{path} answered {e.code}: {said or 'not found'}")
            if attempt == 4:
                raise Unavailable(f"{base}{path} answered {e.code}")
        except (urllib.error.URLError, TimeoutError) as e:
            if attempt == 4:
                raise Unavailable(f"could not reach {base}{path}: {e}")
        time.sleep(2 ** attempt)


def rows_of(text, needs, what):
    reader = csv.DictReader(io.StringIO(text))
    missing = [c for c in needs if c not in (reader.fieldnames or [])]
    if missing:
        raise Unavailable(f"the {what} file has no {', '.join(missing)} column; is it a Daishi study's {what} file?")
    rows = list(reader)
    if what == "trials" and not rows:
        raise Unavailable("the trials file has no trial in it")
    return rows


def local(path):
    try:
        with open(path, encoding="utf-8-sig", newline="") as f:
            return f.read()
    except OSError as e:
        raise Unavailable(f"could not read {path}: {e}")


def num(s):
    return float(s) if s not in ("", None) else None


def plain(s):
    # The files keep a label from reading as a spreadsheet formula with a leading apostrophe.
    return s[1:] if s[:1] == "'" and s[1:2] in ("=", "+", "-", "@", "\t", "\r") else s


def to_analysis(metric, v):
    return -v if metric == "rank" else math.log(v) if metric == "output_tokens_per_turn" else v


def diff_bounds(metric):
    lo, hi = METRIC_RANGE.get(metric, (None, None))
    if metric in ("rank", "output_tokens_per_turn") or lo is None or hi is None:
        return (None, None)
    return (lo - hi, hi - lo)


# ---- Checking against the page --------------------------------------------------

def close(mine, page, tol):
    if mine is None or page is None:
        return mine is None and page is None
    if isinstance(mine, (tuple, list)):
        return isinstance(page, (tuple, list)) and len(mine) == len(page) and all(close(x, y, tol) for x, y in zip(mine, page))
    if isinstance(mine, str) or isinstance(page, str):
        return mine == page
    return abs(mine - page) <= tol


class Ledger:
    """One kind of check, counted: figures recomputed against the page, or the files' account against itself and the page's counts."""

    def __init__(self):
        self.checked = 0
        self.differ = 0

    def same(self, mine, page, tol=0.0):
        ok = close(mine, page, tol)
        self.checked += 1
        self.differ += 0 if ok else 1
        return ok


def pair_of(iv):
    return (iv["lo"], iv["hi"]) if iv else None


def verdict_word(v, page_word):
    if v is None or v == "undecided":
        return v
    family = ("arm_better", "arm_higher") if v == "a" else ("control_better", "control_higher")
    return page_word if page_word in family else family[0]


def recompute_comparison(c, trials_by, metric, rule, level, recipe, arms, ledger, accounts):
    a = c["arm"]
    read = sorted(t for t, cells in trials_by.items() if (cells.get(a) or {}).get("in_comparison") == "yes" and 0 in cells)
    diffs = [to_analysis(metric, num(trials_by[t][a]["value"])) - to_analysis(metric, num(trials_by[t][0]["value"])) for t in read]
    # A pair's map is its seed; a trial whose seed is not known is a map of its own.
    worlds = [trials_by[t][0]["seed"] or f"trial {t}" for t in read]
    places = 4 if metric == "output_tokens_per_turn" else 3
    tol = 1.1 * 10 ** -places
    rd = lambda x: keep_diff(x, places)
    bounds = diff_bounds(metric)
    n = len(diffs)
    reason = c.get("reason") or ""
    # The rule renders no verdict on a reading that is short of its look, an
    # owner's early stop of a fixed rule, or one that lost more trials than its
    # design allows (the page's count and the design's tolerance, which the
    # files do not hold).
    withheld = reason in ("incomplete", "waiting", "more")

    iv = verdict = None
    if rule == "fixed":
        if n >= 3 and reason not in ("owner", "waiting"):
            ci = t_interval(diffs, level, bounds)
            boot_lo, boot_hi = bootstrap(diffs, level)
            m = mean(diffs)
            se = stddev(diffs) / math.sqrt(n)
            significant = math.isfinite(se) and se > 0 and abs(m) > t_critical(n - 1, level) * se
            verdict = "a" if boot_lo > 0 and m > 0 and significant else "b" if boot_hi < 0 and m < 0 and significant else "undecided"
            if ci is not None and math.isfinite(ci[0]):
                iv = (rd(ci[0]), rd(ci[1]))
    elif n >= max(2, recipe.get("first_look") or 2):
        alpha = recipe["alpha"] / (arms - 1) if arms > 2 else recipe["alpha"]
        m = mean(diffs)
        half = sequence_half_width(stddev(diffs), n, alpha, recipe["t_star"])
        if math.isfinite(half):
            verdict = "a" if m - half > 0 else "b" if m + half < 0 else "undecided"
            lo, hi = bounded(m, half, bounds)
            iv = (rd(lo), rd(hi))
    if withheld:
        verdict = None
    word = verdict_word(verdict, c.get("verdict"))

    # The map-level interval and the exact test stand beside a fixed rule's look only.
    looked = rule == "fixed" and verdict is not None
    count = len(set(worlds))
    repeats = count < len(worlds)
    by_world_iv = None
    if repeats and looked and count >= 3:
        w = clustered_interval(diffs, worlds, level, bounds)
        if w is not None and math.isfinite(w[0]) and w[0] != w[1]:
            by_world_iv = (rd(w[0]), rd(w[1]))
    unit = p = None
    if looked:
        unit = "world" if repeats and count > 1 else "pair"
        p = sign_flip(cluster_totals(diffs, worlds) if unit == "world" else diffs)
        if p is None:
            unit = None

    page_bw = c.get("byWorld")
    page_exact = c.get("exact")
    repeated = sum(1 for w in set(worlds) if worlds.count(w) > 1)
    md = keep_mean_diff(diffs, places, metric != "output_tokens_per_turn") if n else None
    # A figure this script did not compute (a verdict the rule withholds, and
    # with it the map-level interval and the exact test, or those two under
    # the anytime rule) is none beside the page's none by the page's own
    # account: that check counts with the files' account, not as a figure.
    tested = ledger if looked else accounts
    judged = ledger if verdict is not None else accounts
    out = sum(1 for cells in trials_by.values() if any(r["pair"] == "out" for r in cells.values()))
    return dict(
        arm=a,
        # Which pairs were read is the page's own account (the in_comparison column), so this count checks the files against the page, not the arithmetic.
        n=dict(recomputed=n, page=c["n"], matches=accounts.same(n, c["n"])),
        mean_diff=dict(recomputed=md, page=c.get("meanDiff"), matches=ledger.same(md, c.get("meanDiff"), tol)),
        interval=dict(level=level, recomputed=iv, page=pair_of(c.get("interval")), matches=ledger.same(iv, pair_of(c.get("interval")), tol)),
        by_world=None if not repeats and not page_bw else dict(
            worlds=count, repeated=repeated, recomputed=by_world_iv, page=pair_of(page_bw.get("interval")) if page_bw else None,
            matches=tested.same((count, repeated, by_world_iv), (page_bw["worlds"], page_bw["repeated"], pair_of(page_bw.get("interval"))) if page_bw else None, tol),
        ),
        exact=dict(unit=unit, recomputed=p, page=page_exact["p"] if page_exact else None,
                   matches=tested.same((unit, p) if p is not None else None, (page_exact["unit"], page_exact["p"]) if page_exact else None, 1e-12)),
        verdict=dict(recomputed=word, page=c.get("verdict"), matches=judged.same(word, c.get("verdict")), reason=reason, lost=c.get("lost"), pairs_out=out),
    )


def recompute_arm(a, rows, seats, metric, page_arm, ledger):
    out = dict(arm=a, score=None, models=[])
    values = [num(r["value"]) for r in rows if r["value"] != ""]
    score = page_arm.get("score")
    if values or score:
        bounds = METRIC_RANGE.get(metric, (None, None))
        m = keep_figure(mean(values), 2) if values else None
        ci = t_interval(values, 0.95, bounds)
        mine = (m, (keep_figure(ci[0], 2), keep_figure(ci[1], 2)) if ci and math.isfinite(ci[0]) else None)
        page = (score["mean"], pair_of(score.get("ci"))) if score else None
        # The file holds each trial's value to 3 decimals and the page's arm mean
        # reads the unrounded values, so the two can differ in the last digit.
        out["score"] = dict(trials=len(values), recomputed=mine, page=page, matches=ledger.same(mine if values else None, page, 0.0101))

    # Each model over its seats in the arm's whole trials, each trial once; guests are left out.
    models = {}
    for s in seats:
        if s["whole_trial"] != "yes" or s["model"] == "guest" or s["fitness_index"] == "":
            continue
        b = models.setdefault(plain(s["model"]), dict(fitness=[], trial=[]))
        b["fitness"].append(float(s["fitness_index"]))
        b["trial"].append(s["trial"])
    page_models = {pm["model"]: pm for pm in page_arm.get("modelStats") or []}
    for model, b in models.items():
        ci = clustered_interval(b["fitness"], b["trial"], 0.95, (0, 100))
        trials = len(set(b["trial"]))
        mine = (trials, len(b["fitness"]), keep_figure(mean(b["fitness"]), 1), (keep_figure(ci[0], 1), keep_figure(ci[1], 1)) if ci else None)
        pm = page_models.pop(model, None)
        page = (pm["trials"], pm["seats"], pm["mean"], pair_of(pm.get("ci"))) if pm else None
        out["models"].append(dict(model=model, recomputed=mine, page=page, matches=ledger.same(mine, page, 0.011)))
    for model, pm in page_models.items():
        page = (pm["trials"], pm["seats"], pm["mean"], pair_of(pm.get("ci")))
        out["models"].append(dict(model=model, recomputed=None, page=page, matches=ledger.same(None, page)))
    return out


def trial_values_from_seats(trials, seats, metric, ledger):
    """Each trial's value as the mean of its seats' readings: the two files read against each other."""
    reading = SEAT_READING.get(metric)
    if reading is None and metric != "survived":
        return dict(checked=0, matching=0, differ=[], skipped=f"{metric} is not a seat score, so each trial's value in the trials file is taken on trust")
    by = {}
    for s in seats:
        by.setdefault((s["trial"], s["arm"]), []).append(s)
    checked = matching = 0
    differ = []
    for r in trials:
        if r["value"] == "":
            continue
        # The metric reads the focal seats, or every seat when the study names none.
        group = [s for s in by.get((r["trial"], r["arm"]), []) if s["focal"] != "no"]
        if metric == "survived":
            per = [None if s["status"] == "" else 0.0 if s["status"] in ("dead", "dormant") else 1.0 for s in group]
        else:
            per = [num(s[reading]) for s in group]
        if not per or any(v is None for v in per):
            continue
        checked += 1
        mine = keep_trial_value(mean(per))
        if ledger.same(mine, float(r["value"]), 0.0011):
            matching += 1
        else:
            differ.append(dict(trial=int(r["trial"]), arm=int(r["arm"]), from_seats=mine, in_trials_file=float(r["value"])))
    return dict(checked=checked, matching=matching, differ=differ, skipped=None)


def trust_list(metric, rule, has_comparisons, from_seats):
    """What a run of this script takes on trust for this study, printed at the end of every run."""
    out = [
        "each seat's scores, as the match's report published them (recompute.py recomputes any one match from its public record)",
        "which trials played whole",
    ]
    if has_comparisons:
        out.append("which pairs each comparison read: the pair and in_comparison columns are the page's own account, so a pair left out and the trial an early stop fell on are not rechecked")
        out.append("why a verdict is withheld (trials lost over the design's tolerance, an owner's stop, a look not reached): the page's reason")
    if from_seats and from_seats.get("skipped"):
        out.append(f"each trial's value in the trials file, since {metric} is not a seat score")
    if rule == "anytime":
        out.append("the anytime rule's tuning constant (t_star), from the study's recipe")
    return out


NOT_RECOMPUTED = "secondary outcomes, the detectable gap, the debrief counts, the output-per-turn table and spend"


def meaning(metric):
    also = "" if metric == "fitness_index" else ", and its per-model Fitness Index means,"
    return (f"A full match means the page's figures on {metric}{also} follow from the two files by the platform's own arithmetic, "
            "which this script ports. It does not show that the files are a true record of play, that a score measures "
            "what its name says, or that this analysis suits the question.")


# ---- Printing -------------------------------------------------------------------

def pct(level):
    p = level * 100
    return str(round(p)) if abs(p - round(p)) < 1e-9 else f"{p:.1f}"


def fmt(v):
    if v is None:
        return "none"
    if isinstance(v, tuple):
        return f"{fmt(v[0])} to {fmt(v[1])}"
    return f"{v:g}" if isinstance(v, float) else str(v)


def same(ok):
    return "same" if ok else "DIFFERS"


def print_report(head, labels, level, comparisons, arm_reports, from_seats, figures, accounts, trusted, metric, tolerance):
    print(head + "\n")
    control = labels.get(0, "the control")
    for c in comparisons:
        print(f"{labels.get(c['arm'], 'arm %d' % c['arm'])} against {control}: {c['n']['recomputed']} pairs read as the trials file marks them, page {c['n']['page']}  {same(c['n']['matches'])}")
        print(f"  {'':36}{'recomputed':>22}{'page':>22}")
        line = lambda what, x: print(f"  {what:36}{fmt(x['recomputed']):>22}{fmt(x['page']):>22}  {same(x['matches'])}")
        line("mean paired difference", c["mean_diff"])
        line(f"{pct(level)}% interval", c["interval"])
        if c["by_world"]:
            line(f"{pct(level)}% interval, each of {c['by_world']['worlds']} maps once", c["by_world"])
        line("exact sign-flip p" + (f", over {c['exact']['unit']}s" if c["exact"]["unit"] else ""), c["exact"])
        line("verdict", c["verdict"])
        if c["verdict"]["reason"] == "incomplete":
            lost, out = c["verdict"]["lost"], c["verdict"]["pairs_out"]
            allows = "the design allows" if tolerance is None else f"the {tolerance} the design allows"
            print(f"  No verdict: the page counts {lost} trial{'' if lost == 1 else 's'} lost, more than {allows}, and the trials file marks {out} pair{'' if out == 1 else 's'} out. The lost count and the tolerance are the page's figures; the files hold neither.")
        print()
    print("Each arm over its trials, 95% interval:")
    for r in arm_reports:
        s = r["score"]
        if not s:
            continue
        mine = s["recomputed"] or (None, None)
        page = s["page"] or (None, None)
        print(f"  {labels.get(r['arm'], 'arm %d' % r['arm'])}: {s['trials']} trials, mean {fmt(mine[0])} ({fmt(mine[1])}); page {fmt(page[0])} ({fmt(page[1])})  {same(s['matches'])}")
    print("\nEach model's Fitness Index over its seats in whole trials, 95% interval counting each trial once:")
    for r in arm_reports:
        for md in r["models"]:
            mine = md["recomputed"] or (None, None, None, None)
            page = md["page"] or (None, None, None, None)
            print(f"  {labels.get(r['arm'], 'arm %d' % r['arm'])}, {md['model']}: {fmt(mine[1])} seats in {fmt(mine[0])} trials, mean {fmt(mine[2])} ({fmt(mine[3])}); "
                  f"page {fmt(page[1])} seats, mean {fmt(page[2])} ({fmt(page[3])})  {same(md['matches'])}")
    if from_seats and from_seats["skipped"]:
        print(f"\nEach trial's value from its seats: not rebuilt; {from_seats['skipped']}.")
    elif from_seats:
        print(f"\nEach trial's value from its seats: {from_seats['matching']} of {from_seats['checked']} match the trials file.")
        for d in from_seats["differ"][:10]:
            print(f"  trial {d['trial']}, arm {d['arm']}: {fmt(d['from_seats'])} from its seats, {fmt(d['in_trials_file'])} in the trials file  DIFFERS")
    print(f"\n{figures.checked - figures.differ} of {figures.checked} recomputed figures match the page.")
    if accounts.checked:
        print(f"{accounts.checked - accounts.differ} of {accounts.checked} checks of the files against each other and against the page's own counts agree.")
    print("\nTaken on trust:")
    for t in trusted:
        print(f"  - {t}")
    print(f"Not recomputed: {NOT_RECOMPUTED}.")
    print(meaning(metric))


def main():
    ap = argparse.ArgumentParser(description="Recompute a published Daishi study's reading of its metric from its data files.")
    ap.add_argument("study", help="the study's page, https://daishi.ai/u/<handle>/studies/<id>")
    ap.add_argument("--base", help="the world to read from (default: the page's own, else https://daishi.ai)")
    ap.add_argument("--trials", help="read the trials file from this path instead of downloading it")
    ap.add_argument("--seats", help="read the seats file from this path instead of downloading it")
    ap.add_argument("--json", action="store_true", help="print the results as JSON")
    args = ap.parse_args()
    m = STUDY_RE.match(args.study.strip())
    if not m:
        print("Give the study's page, as https://daishi.ai/u/<handle>/studies/<id>.", file=sys.stderr)
        return 2
    base = (args.base or m.group(1) or "https://daishi.ai").rstrip("/")
    handle, study_id = m.group(2), m.group(3)
    q = lambda s: urllib.parse.quote(s, safe="")
    page_path = f"/u/{q(handle)}/studies/{q(study_id)}"

    try:
        trials = rows_of(local(args.trials) if args.trials else fetch(base, f"{page_path}/trials.csv", "text/csv"), TRIAL_NEEDS, "trials")
        seats = rows_of(local(args.seats) if args.seats else fetch(base, f"{page_path}/seats.csv", "text/csv"), SEAT_NEEDS, "seats")
        if trials[0]["study"] != study_id:
            raise Unavailable(f"the trials file is study {trials[0]['study']}'s, not {study_id}'s")
        profile = json.loads(fetch(base, f"/api/u/{q(handle)}", "application/json"))
        study = next((s for s in profile.get("studies") or [] if s.get("id") == study_id), None)
        if study is None:
            raise Unavailable(f"{base}{page_path} is not among @{handle}'s published studies")
        metric = trials[0]["metric"]
        rule = trials[0]["stopping_rule"] or None
        level = num(trials[0]["comparison_level"])
        recipe = None
        if rule == "anytime":
            recipe = (json.loads(fetch(base, f"/api/studies/{q(study_id)}/recipe", "application/json")).get("recipe") or {}).get("stopping")
            if not recipe or recipe.get("t_star") is None:
                raise Unavailable("the anytime rule's sequence needs the design's tuning constant (t_star), and this study's recipe does not carry it")
    except Unavailable as e:
        print(e, file=sys.stderr)
        return 2

    # Two kinds of check, counted apart: figures recomputed from the files
    # against the page's, and the files' own account (the pairs each
    # comparison read, each trial's value from its seats) against itself
    # and the page's counts, which recomputes nothing on the page.
    figures, accounts = Ledger(), Ledger()
    page_arms = {pa["index"]: pa for pa in study.get("arms") or []}
    arms = sorted({int(r["arm"]) for r in trials} | set(page_arms))
    labels = {int(r["arm"]): plain(r["arm_label"]) for r in trials}
    trials_by = {}
    for r in trials:
        trials_by.setdefault(int(r["trial"]), {})[int(r["arm"])] = r
    page_comparisons = study.get("comparisons") or []
    comparisons = [recompute_comparison(c, trials_by, metric, rule, level, recipe, len(page_arms) or len(arms), figures, accounts) for c in page_comparisons]
    arm_reports = [
        recompute_arm(a, [r for r in trials if int(r["arm"]) == a], [s for s in seats if int(s["arm"]) == a], metric, page_arms.get(a, {}), figures)
        for a in arms
    ]
    from_seats = trial_values_from_seats(trials, seats, metric, accounts)
    trusted = trust_list(metric, rule, bool(page_comparisons), from_seats)
    all_match = figures.differ == 0 and accounts.differ == 0

    if args.json:
        print(json.dumps(dict(study=study_id, handle=handle, base=base, kind=study.get("kind"), metric=metric, stopping_rule=rule, comparison_level=level,
                              comparisons=comparisons, arms=arm_reports, trial_values_from_seats=from_seats,
                              figures=dict(checked=figures.checked, differ=figures.differ),
                              consistency=dict(checked=accounts.checked, differ=accounts.differ),
                              trusted=trusted, not_recomputed=NOT_RECOMPUTED, meaning=meaning(metric)), indent=1))
        return 0 if all_match else 1

    head = f"Study {study_id} by @{handle}: {study.get('kind')} on {metric}, {len(arms)} arm{'' if len(arms) == 1 else 's'}"
    if rule:
        head += f", {rule} rule, each comparison at {pct(level)}%"
    print_report(f"{head}, read from {base}", labels, level, comparisons, arm_reports, from_seats, figures, accounts, trusted, metric, study.get("lossTolerance"))
    return 0 if all_match else 1


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