#!/usr/bin/env python3
"""Run every engine on the corpus, hold each answer to its reference within its tolerance, and
say what passed, what failed and by how much.

    compare.py --tier ci|full [--out RESULTS.json] [--jobs N] [--inject FAULT] [--only ID,...]

Exit status 0 when every check passes, 1 when one fails, 2 when the run itself could not be
made. `--inject` corrupts one engine's answers after it gives them, or makes the engine under
test refuse, crash, hang or accept what it must not, the way a wrong engine would, so that a run
can be seen to fail; FAULTS lists what it can do.

Every tolerance is derived, never fitted, and tolerance.py says from what. A check another engine
cannot make (an engine without the gate, a method that does not take the circuit) is recorded as
`not run`, with the engine's own words, and is not a pass. The engine under test is held to more:
a crash, a time-out, or a refusal that corpus/expected-refusals.json does not list in those words
fails the run, and so does an answer to a circuit that file says it must refuse.
"""
import argparse
import hashlib
import itertools
import json
import math
import os
import platform
import shutil
import signal
import subprocess
import sys
import tempfile
import threading
import time
from concurrent.futures import ThreadPoolExecutor
from fractions import Fraction

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
import corpus  # noqa: E402
import decoding  # noqa: E402
import exact  # noqa: E402
import tolerance as tol  # noqa: E402

ROOT = os.path.dirname(HERE)
WAI = os.path.dirname(ROOT)

STATE_ENGINES = ["wai.exact", "wai.statevector", "aer.statevector", "aer.qasm_import", "quimb.dense", "quimb.tensor_network", "ddsim.statevector", "qsim.statevector"]
COUNT_ENGINES = ["wai.auto", "wai.statevector", "aer.statevector", "ddsim.statevector", "qsim.statevector"]
CLIFFORD = {"id", "x", "y", "z", "h", "s", "sdg", "sx", "sxdg", "cx", "cy", "cz", "swap"}
# Every sampled test of a run (counts, error-correction rates, stabilizer outcomes, and the
# estimates of the noisy, mitigated and sampled methods) shares this: a correct engine fails a run
# on sampling with probability at most this, a quarter to each kind.
FAMILY_WISE_ALPHA = 1e-6
SHARES = 4
# The closed-form checks of the integer reference, in units of 2^-256: far above its own rounding
# (a few units a gate), far below anything an engine is held to.
REFERENCE_UNITS = 2**40
# The most amplitude updates (the op-log's length times 2^n) the referee reruns bit for bit: its
# pure-Python state vector takes about 0.2 µs an update on the build host, so this is a few minutes.
REFEREE_WORK = 2**30
# mqt.ddsim's angle tolerance (its parameter_tolerance setting): an angle below it may be read as 0.
DDSIM_PARAMETER_TOLERANCE = 1e-13

FAULTS = {
    "amplitude": "aer.statevector: one amplitude of every state moved by 1e-9",
    "bit_order": "wai.statevector: every state read with its qubits in reverse order",
    "counts": "stim.tableau: 5% of the shots of the most frequent outcome moved to the least frequent outcome that can occur",
    "counts_spread": "aer.statevector: every outcome's count moved by a tenth of itself, up and down in turn, the total kept (an L1 distance of about 0.1)",
    "dem": "wai.qec: the most likely mechanism's probability raised by one part in 10^6",
    "qec_rate": "wai.qec: the sampler's most frequent detector count raised by a quarter",
    "stabilizer": "stim.tableau: the sign of the last stabilizer flipped",
    "fixed_point": "wai.exact: the lowest bit of one fixed-point amplitude flipped",
    "drop_gate": "wai.exact: the op-log's last phase dropped, and the integers and hash made again from what is left",
    "drop_phase": "wai.statevector: the state of the circuit with its small-angle gate left out",
    "refuse": "wai.mps: the engine refuses every job",
    "crash": "wai.auto: the engine exits with status 101 and no answer",
    "hang": "wai.qec: the engine's sampler never answers",
    "accept": "wai.exact: every angle moved to the nearest m/2^32 of a turn before the engine reads it",
    "qec_words": "wai.qec: the engine refuses each model with a line number, in words no list holds",
    "qec_crash": "wai.qec: the engine crashes on each model",
    "exact_words": "wai.exact: the engine refuses each state, in words other than its list's",
    "qec_short": "wai.qec: the sampler's answer lists no observable",
    "exact_turn": "wai.exact: every state turned by a global phase of 0.5 rad",
    "phase_flip": "wai.statevector: every state turned by a half turn, as a phase applied with the wrong sign, or not at all",
    "mps_phase": "aer.mps: every amplitude turned by 0.5 rad",
    "counts_collapse": "stim.tableau: every shot given to one outcome",
    "pauli_bound": "wai.pauli: every expectation value, at every threshold of its ladder, moved by twice its truncation bound and 10^-9",
    "stabrank_amp": "wai.stabrank: every amplitude of an exact sum moved by 10^-9",
    "cutrank_amp": "wai.cutrank: one amplitude moved by 10^-8",
    "trajectory_mean": "wai.trajectories: every expectation value moved by 0.05",
    "density_value": "wai.density_matrix: every expectation value of a noisy run moved by 10^-8",
    "frame_sign": "wai.stabilizer: every noisy expectation value negated",
    "zne_raw": "wai.mitigate: each zero-noise extrapolation replaced by the value at scale 1",
    "readout_raw": "wai.mitigate: each readout-mitigated quasi-probability replaced by the raw frequency",
    "vd_purity": "wai.mitigate: the stated purity raised by 10^-6",
    "pulse_p1": "wai.pulse: P(1) moved by ten times its stated bound and 10^-12",
    "decode_flip": "wai.matching: every prediction's observables flipped",
    "uf_flip": "wai.union_find: every prediction's observables flipped",
    "member_bit": "wai.fabric: the lowest bit of one amplitude of every GPU member's state flipped",
}
# The faults the engine's adapter carries out (engines.py), rather than this file.
DRIVER_FAULTS = {"refuse": "wai.mps", "crash": "wai.auto", "hang": "wai.qec", "accept": "wai.exact", "qec_words": "wai.qec", "qec_crash": "wai.qec", "exact_words": "wai.exact"}


def is_clifford(e):
    for s in e["steps"]:
        for x in s.get("then", [s]):
            if "g" in x and x["g"] not in CLIFFORD:
                return False
    return True


def has_dynamic(e):
    k = len(e["steps"])
    while k and "measure" in e["steps"][k - 1]:
        k -= 1
    return any("g" not in s for s in e["steps"][:k]) or k < len(e["steps"]) and e["expect"] != "distribution"


def group_rss_bytes(pgid):
    """The resident memory of every process in a process group, as the system counts each."""
    try:
        out = subprocess.run(["ps", "-A", "-o", "pgid=,rss="], capture_output=True, text=True, timeout=30).stdout
    except (OSError, subprocess.SubprocessError):
        return 0
    total = 0
    for line in out.splitlines():
        f = line.split()
        if len(f) == 2 and f[0] == str(pgid):
            total += int(f[1]) * 1024
    return total


class Runner:
    def __init__(self, work, python, timeout, inject, memory_gb=None):
        self.work = work
        self.python = python
        self.timeout = timeout
        self.inject = inject
        # The memory guard: a job whose processes together hold more than this is stopped, with
        # every process it started, and recorded like a time-out.
        self.memory = int(memory_gb * 2**30) if memory_gb else None
        self.counter = itertools.count(1)
        self.lock = threading.Lock()

    def run(self, engine, job):
        """One engine on one job in a process group of its own: its answer and wall time. On a
        time-out the whole group is killed, the engine's driver with its adapter."""
        with self.lock:
            k = next(self.counter)
        tag = "%s-%s-%d" % (job.get("id", "job"), engine, k)
        jp = os.path.join(self.work, tag + ".job.json")
        ap = os.path.join(self.work, tag + ".answer.json")
        if job["kind"] == "state" and "state_out" not in job:
            job = dict(job, state_out=os.path.join(self.work, tag + ".state"))
        if self.inject in DRIVER_FAULTS and DRIVER_FAULTS[self.inject] == engine:
            job = dict(job, inject=self.inject)
        with open(jp, "w") as f:
            json.dump(job, f)
        t0 = time.perf_counter()
        ep = os.path.join(self.work, tag + ".stderr")
        with open(ep, "wb") as ef:
            p = subprocess.Popen([self.python, "-I", os.path.join(HERE, "engines.py"), engine, jp, ap], stdout=subprocess.DEVNULL, stderr=ef, start_new_session=True)
            stopped = None
            while p.poll() is None:
                try:
                    p.wait(timeout=1.0)
                except subprocess.TimeoutExpired:
                    pass
                if p.poll() is not None:
                    break
                if time.perf_counter() - t0 > self.timeout:
                    stopped = {"ok": False, "timeout": self.timeout}
                elif self.memory and group_rss_bytes(p.pid) > self.memory:
                    stopped = {"ok": False, "guard": "stopped at the memory guard: its processes held more than %.0f GiB" % (self.memory / 2**30)}
                if stopped:
                    try:
                        os.killpg(p.pid, signal.SIGKILL)
                    except ProcessLookupError:
                        pass
                    p.wait()
                    return stopped
        err = open(ep, "rb").read()
        wall = time.perf_counter() - t0
        try:
            ans = json.load(open(ap))
        except Exception:  # noqa: BLE001
            return {"ok": False, "error": "no answer; exit %d: %s" % (p.returncode, err.decode(errors="replace")[-1500:])}
        ans["wall_seconds"] = wall
        ans["job"] = job
        return self.corrupt(engine, ans)

    def corrupt(self, engine, ans):
        """The injected fault, applied as the engine's own answer."""
        f = self.inject
        if not f or not ans.get("ok"):
            return ans
        import numpy as np

        job = ans["job"]
        if f == "amplitude" and engine == "aer.statevector" and job["kind"] == "state":
            v = np.fromfile(job["state_out"], dtype=np.complex128)
            v[len(v) // 2] += 1e-9
            v.tofile(job["state_out"])
        elif f == "bit_order" and engine == "wai.statevector" and job["kind"] == "state" and job["n"] > 1:
            import engines

            v = np.fromfile(job["state_out"], dtype=np.complex128)
            engines.reverse_qubits(v, job["n"]).tofile(job["state_out"])
        elif f == "counts" and engine == "stim.tableau" and job["kind"] == "counts":
            c = ans["counts"]
            possible = [k for k, p in exact.distribution(job["n"], job["bits"], job["steps"]).items() if p > 0]
            top = max(c, key=c.get)
            dest = min((k for k in possible if k != top), key=lambda k: c.get(k, 0), default=None)
            if dest is not None:
                moved = c[top] // 20
                c[top] -= moved
                c[dest] = c.get(dest, 0) + moved
        elif f == "counts_spread" and engine == "aer.statevector" and job["kind"] == "counts":
            c = ans["counts"]
            keys = sorted(c)
            total = sum(c.values())
            for i, k in enumerate(keys):
                c[k] += (c[k] // 10) * (1 if i % 2 == 0 else -1)
            c[keys[0]] += total - sum(c.values())
        elif f == "dem" and engine == "wai.qec" and job["kind"] == "qec_model" and ans["mechanisms"]:
            m = max(ans["mechanisms"], key=lambda x: x[0])
            m[0] *= 1 + 1e-6
        elif f == "qec_rate" and engine == "wai.qec" and job["kind"] == "qec_sample" and ans.get("detectors"):
            d = ans["detectors"]
            i = max(range(len(d)), key=lambda j: d[j])
            d[i] += d[i] // 4
        elif f == "stabilizer" and engine == "stim.tableau" and job["kind"] == "stabilizers":
            last = ans["stabilizers"][-1]
            ans["stabilizers"][-1] = ("-" if last[0] == "+" else "+") + last[1:]
        elif f == "fixed_point" and engine == "wai.exact" and job["kind"] == "state":
            fx = ans["settings"]["fixed"]
            a = np.load(fx)
            a[0, 0] ^= 1
            np.save(fx, a)
        elif f == "drop_gate" and engine == "wai.exact" and job["kind"] == "state":
            ops = ans["ops"]
            last = max((i for i, op in enumerate(ops) if op[0] == 9), default=None)
            if last is not None:
                del ops[last]
                ref = referee()
                amps = ref.simulate(job["n"], [ref.gate(op, t, tuple(c), p) for op, c, t, p in ops], ref.derive_v2()[0])
                ints = np.array([[int(a), int(b)] for a, b in amps], dtype=np.int64)
                np.save(ans["settings"]["fixed"], ints)
                digest = ref.statevector_hash(job["n"], amps)
                ans["state_hash"] = digest.hex() if isinstance(digest, (bytes, bytearray)) else digest
                ((ints[:, 0] / float(1 << 30)) + 1j * (ints[:, 1] / float(1 << 30))).tofile(job["state_out"])
        elif f == "qec_short" and engine == "wai.qec" and job["kind"] == "qec_sample":
            ans["observables"] = []
        elif f == "exact_turn" and engine == "wai.exact" and job["kind"] == "state":
            fx = ans["settings"]["fixed"]
            a = np.load(fx)
            z = (a[:, 0] + 1j * a[:, 1]) * complex(math.cos(0.5), math.sin(0.5))
            a[:, 0], a[:, 1] = np.round(z.real).astype(np.int64), np.round(z.imag).astype(np.int64)
            np.save(fx, a)
            ((a[:, 0] / float(1 << 30)) + 1j * (a[:, 1] / float(1 << 30))).tofile(job["state_out"])
        elif f == "phase_flip" and engine == "wai.statevector" and job["kind"] == "state":
            v = np.fromfile(job["state_out"], dtype=np.complex128)
            (-v).tofile(job["state_out"])
        elif f == "mps_phase" and engine == "aer.mps" and job["kind"] == "amplitudes":
            t = complex(math.cos(0.5), math.sin(0.5))
            ans["amplitudes"] = {b: [(complex(*a) * t).real, (complex(*a) * t).imag] for b, a in ans["amplitudes"].items()}
        elif f == "counts_collapse" and engine == "stim.tableau" and job["kind"] == "counts":
            c = ans["counts"]
            ans["counts"] = {sorted(c)[0]: sum(c.values())}
        elif f == "pauli_bound" and engine == "wai.pauli" and job["kind"] == "expectations":
            for x in ans["result"]["expectations"]:
                for r in [x] + x.get("ladder", []):
                    r["value"] += 2 * r["bound"] + 1e-9
        elif f == "stabrank_amp" and engine == "wai.stabrank" and job["kind"] == "amplitudes" and (ans.get("accuracy") or {}).get("class") == "exact":
            ans["amplitudes"] = {b: [a[0] + 1e-9, a[1]] for b, a in ans["amplitudes"].items()}
        elif f == "cutrank_amp" and engine == "wai.cutrank" and job["kind"] == "amplitudes":
            b = sorted(ans["amplitudes"])[0]
            ans["amplitudes"][b] = [ans["amplitudes"][b][0] + 1e-8, ans["amplitudes"][b][1]]
        elif f == "trajectory_mean" and engine == "wai.trajectories" and job["kind"] == "expectations":
            for x in ans["result"].get("expectations", []):
                x["value"] += 0.05
        elif f == "density_value" and engine == "wai.density_matrix" and job["kind"] == "expectations":
            for x in ans["result"].get("expectations", []):
                x["value"] += 1e-8
        elif f == "frame_sign" and engine == "wai.stabilizer" and job["kind"] == "expectations":
            for x in ans["result"].get("expectations", []):
                x["value"] = -x["value"]
        elif f == "zne_raw" and engine == "wai.mitigate" and job["kind"] == "mitigated" and ans["result"]["mitigation"]["method"] == "zne":
            for x in ans["result"]["mitigation"]["observables"]:
                x["mitigated"] = x["raw"]
        elif f == "readout_raw" and engine == "wai.mitigate" and job["kind"] == "mitigated" and ans["result"]["mitigation"]["method"] == "readout":
            for x in ans["result"]["mitigation"]["distribution"]:
                x["mitigated"] = x["raw"]
        elif f == "vd_purity" and engine == "wai.mitigate" and job["kind"] == "mitigated" and ans["result"]["mitigation"]["method"] == "vd":
            ans["result"]["mitigation"]["details"]["purity"] += 1e-6
        elif f == "pulse_p1" and engine == "wai.pulse" and job["kind"] == "pulse":
            r = ans["result"]
            r["pulse"]["populations"]["1"] += 10 * r["accuracy"]["bound"] + 1e-12
        elif f in ("decode_flip", "uf_flip") and engine == {"decode_flip": "wai.matching", "uf_flip": "wai.union_find"}[f] and job["kind"] == "decode" and job["syndromes"]:
            ans["predicted"] = [[0] if not p else [] for p in ans["predicted"]]
        elif f == "member_bit" and engine == "wai.fabric" and job["kind"] == "fabric":
            import numpy as np

            for m in ans["members"]:
                if m["family"] != "cpu" and m.get("state_out"):
                    a = np.fromfile(m["state_out"], dtype="<i8")
                    a[0] ^= 1
                    a.tofile(m["state_out"])
        elif f == "drop_phase" and engine == "wai.statevector" and job["kind"] == "state" and job.get("probe") is not None:
            steps = [s for i, s in enumerate(job["steps"]) if i != job["probe"]]
            w = exact.run_unitary(job["n"], steps)
            np.array([complex(float(Fraction(r, exact.ONE)), float(Fraction(i, exact.ONE))) for r, i in zip(w.re, w.im)]).tofile(job["state_out"])
        return ans


_REFEREE = None


def referee():
    """The conformance corpus's independent integer reader of the bit-exact method (section 3.0)."""
    global _REFEREE
    if _REFEREE is None:
        import importlib.util

        spec = importlib.util.spec_from_file_location("referee", os.path.join(WAI, "quantum-conformance", "tools", "referee.py"))
        _REFEREE = importlib.util.module_from_spec(spec)
        spec.loader.exec_module(_REFEREE)
    return _REFEREE


class Report:
    def __init__(self):
        self.checks = []

    def add(self, entry, engine, quantity, measured, tolerance, ok, note=None, status=None):
        c = {"entry": entry, "engine": engine, "quantity": quantity, "measured": measured, "tolerance": tolerance, "status": status or ("pass" if ok else "fail")}
        if note:
            c["note"] = note
        self.checks.append(c)
        return c

    def not_run(self, entry, engine, quantity, why):
        return self.add(entry, engine, quantity, None, None, True, note=why, status="not run")


# ------------------------------------------------------------------ the engine under test, held to its answers


class Expected:
    """corpus/expected-refusals.json: the refusals by the engine under test a run accepts."""

    def __init__(self):
        self.rows = json.load(open(os.path.join(ROOT, "corpus", "expected-refusals.json")))["refusals"]
        self.by = {(r["entry"], r["engine"]): r for r in self.rows}
        self.answered = set()

    def get(self, entry, engine):
        return self.by.get((entry, engine))


EXPECTED = None


def under_test(engine):
    return engine.startswith("wai.")


def words(ans):
    return ans.get("refused") or ans.get("unsupported") or ans.get("error") or ans.get("guard") or ("did not finish within %s s" % ans.get("timeout"))


def no_answer(rep, entry, eng, ans, quantity, fail_on_error=False):
    """An engine gave no answer. From the engine under test, a refusal that expected-refusals.json
    lists for this entry, in its listed words, is recorded as `refused`; anything else (a crash, a
    time-out, a refusal not listed, or in other words) fails the run. From another engine it is
    `not run`, with its words, or a failure where it crashed on a check it must make."""
    why = str(words(ans))
    if under_test(eng):
        x = EXPECTED.get(entry, eng)
        if x and ans.get("refused") and x["text"] in str(ans["refused"]):
            return rep.add(entry, eng, "refuses, as listed: " + x["reason"], ans["refused"], x["text"], True, status="refused")
        kind = "it crashed" if ans.get("error") else "it did not finish" if ans.get("timeout") or ans.get("guard") else "a refusal no list holds" if not x else "a refusal in other words than its list's"
        return rep.add(entry, eng, "%s: the engine under test gave no answer (%s)" % (quantity, kind), why[:1500], None, False)
    if fail_on_error and ans.get("error"):
        return rep.add(entry, eng, quantity, None, None, False, note=why)
    return rep.not_run(entry, eng, quantity, why)


def answered(rep, entry, eng):
    """The engine under test answered. Where it must refuse, that fails the run."""
    if under_test(eng):
        x = EXPECTED.get(entry, eng)
        if x:
            EXPECTED.answered.add((entry, eng))
            if x["must"] == "refuse":
                rep.add(entry, eng, "accepts a circuit it must refuse: " + x["reason"], True, False, False)
                return False
    return True


# ------------------------------------------------------------------ states


def load_state(path):
    import numpy as np

    return np.fromfile(path, dtype=np.complex128)


def exact_distance(ref, vec):
    """(strict, modulo a global phase, the phase) between the integer reference and a state."""
    re = exact.ints_from_floats([float(x) for x in vec.real])
    im = exact.ints_from_floats([float(x) for x in vec.imag])
    d, dm, phase, _, _ = exact.distance(ref.re, ref.im, re, im)
    return d, dm, phase


def float_distance(a, b):
    """(strict, modulo a global phase) between two double states, the sums taken exactly enough:
    math.fsum over the squares."""
    import numpy as np

    def norm(v):
        return math.sqrt(math.fsum((v.real**2).tolist()) + math.fsum((v.imag**2).tolist()))

    d = norm(a - b)
    # Turn b by the phase that best aligns it with a, then measure the difference itself:
    # sqrt(|a|^2 + |b|^2 - 2|<b|a>|) would lose half its digits to cancellation.
    ip = np.vdot(b, a)
    turn = ip / abs(ip) if abs(ip) > 0 else 1.0
    dm = norm(a - turn * b)
    return d, dm


def ref_vector(ref):
    import numpy as np

    return np.array([complex(float(Fraction(r, exact.ONE)), float(Fraction(i, exact.ONE))) for r, i in zip(ref.re, ref.im)])


def state_tolerance(engine, ans, e):
    gates = tol.gate_count(e["steps"])
    n = e["n"]
    if engine == "wai.exact":
        return tol.fixed_point_bound(n, ans["ops"]), "the extension's published bound (quantum-sim.md section 5), summed over the op-log the engine ran"
    if engine in ("wai.statevector",):
        return tol.double_state(gates, n, width=3), "double precision, fusion width 3"
    if engine in ("aer.statevector", "aer.qasm_import"):
        g = max(gates, ans.get("settings", {}).get("gates", gates))
        width = 5 if n >= 14 else 2
        return tol.double_state(g, n, width=width), "double precision, fusion up to %d qubits from 14 qubits" % 5
    if engine == "quimb.dense":
        return tol.double_state(max(gates, ans.get("gates", gates)), n, width=3), "double precision, gate by gate"
    if engine == "quimb.tensor_network":
        w = ans.get("settings", {}).get("contraction_width") or n
        return tol.contraction(ans.get("gates", gates) + n, min(float(w), 40.0)), "double precision, a contraction of width %.1f" % float(w)
    if engine == "ddsim.statevector":
        return tol.decision_diagram(gates, n), "decision-diagram edge weights merged within 2^-42"
    if engine == "qsim.statevector":
        return tol.single_state(gates, n, width=2), "single precision, fusion width 2"
    raise KeyError(engine)


def run_oplog(n, ops):
    """The bit-exact method's op-log run in double precision: its ideal state, to a double's
    rounding. Base gates are quantum-sim.md's (section 3.0): I X Y Z H S S† T T† P = 0..9, P(k)
    the phase e^(2πi/2^k); qubit 0 is the lowest bit of the index."""
    import numpy as np

    r = 1 / math.sqrt(2)
    fixed = {0: ((1, 0), (0, 1)), 1: ((0, 1), (1, 0)), 2: ((0, -1j), (1j, 0)), 3: ((1, 0), (0, -1)), 4: ((r, r), (r, -r)), 5: ((1, 0), (0, 1j)), 6: ((1, 0), (0, -1j)), 7: ((1, 0), (0, complex(r, r))), 8: ((1, 0), (0, complex(r, -r)))}
    v = np.zeros(1 << n, dtype=np.complex128)
    v[0] = 1
    idx = np.arange(1 << n)
    pairs = {}
    for base, controls, target, param in ops:
        m = ((1, 0), (0, np.exp(2j * math.pi / (1 << param)) if param > 2 else (-1 if param == 1 else 1j))) if base == 9 else fixed[base]
        key = (tuple(controls), target)
        if key not in pairs:
            mask = 0
            for c in controls:
                mask |= 1 << c
            lo = idx[((idx & mask) == mask) & ((idx >> target) & 1 == 0)]
            pairs[key] = (lo, lo | (1 << target))
        lo, hi = pairs[key]
        a, b = v[lo], v[hi]
        v[lo], v[hi] = m[0][0] * a + m[0][1] * b, m[1][0] * a + m[1][1] * b
    return v


def oplog_check(rep, e, n, ops, want):
    """The op-log the bit-exact method ran is the circuit: the op-log's ideal state against the
    circuit's, modulo the global phase the lowering drops. Within the double rounding of that run,
    the distance the reader's rule moves angles (corpus.reader_class: a phase moved by δ turns
    changes a state by at most 2πδ), and 2^-50 for the reference's conversion to double. The
    referee then ties the integers to the op-log, so this ties them to the circuit."""
    cls, moved = corpus.reader_class(e)
    _, dm = float_distance(want, run_oplog(n, ops))
    bound = tol.double_state(len(ops), n, width=1) + 2 * math.pi * float(moved) + 2.0**-50
    note = "the reader moved its angles by %.3g of a turn in all, as its rule has it" % float(moved) if cls == "moved" else None
    rep.add(e["id"], "wai.exact", "the op-log is the circuit: the op-log's ideal state against the circuit's, modulo a global phase", dm, bound, dm <= bound, note=note)


def referee_hash(n, ops, fixed, engine_hash):
    """The referee's integer simulator run on the engine's op-log must give the integers the
    engine wrote, and the referee's BLAKE3 of them the hash the engine reported."""
    ref = referee()
    table = ref.derive_v2()[0]
    oplist = [ref.gate(op, t, tuple(c), p) for op, c, t, p in ops]
    amps = ref.simulate(n, oplist, table)
    mine = [(int(a), int(b)) for a, b in fixed]
    same = [tuple(x) for x in amps] == mine
    digest = ref.statevector_hash(n, amps)
    digest = digest.hex() if isinstance(digest, (bytes, bytearray)) else digest
    if not same:
        return {"match": False, "note": "the referee's state differs from the engine's"}
    if digest != engine_hash:
        return {"match": False, "note": "the referee's hash %s, the engine's %s" % (digest, engine_hash)}
    return {"match": True}


def exact_answer(rep, e, n, ans, want, quantity="state"):
    """The bit-exact method's checks beside its distance: the op-log is the circuit, and the
    referee reruns the op-log to the engine's integers and hash (to 14 qubits)."""
    import numpy as np

    fixed = np.load(ans["settings"]["fixed"])
    oplog_check(rep, e, n, ans["ops"], want)
    quantity = "state hash, recomputed by the referee from the op-log and from the integers written"
    if len(ans["ops"]) << n <= REFEREE_WORK:
        h = referee_hash(n, ans["ops"], fixed, ans["state_hash"])
        rep.add(e["id"], "wai.exact", quantity, h["match"], True, h["match"], note=h.get("note"))
    else:
        rep.not_run(e["id"], "wai.exact", quantity, "not rerun: %d ops on 2^%d amplitudes is past the referee's 2^30 updates" % (len(ans["ops"]), n))
    return fixed


def reference_checks(e, ref, rep):
    """A closed form, where there is one, holds the integer reference itself."""
    n = e["n"]
    unit = "units of 2^-256"
    if "qft_x" in e:
        worst = 0
        for y in range(1 << n):
            qr, qi = exact.qft_amplitude(n, e["qft_x"], y)
            worst = max(worst, abs(qr - ref.re[y]) + abs(qi - ref.im[y]))
        rep.add(e["id"], "exact", "integer reference against the closed form, " + unit, worst, REFERENCE_UNITS, worst <= REFERENCE_UNITS)
    if e["family"] == "ghz":
        r = exact.rsqrt2()
        worst = max(abs(ref.re[i] - (r if i in (0, (1 << n) - 1) else 0)) + abs(ref.im[i]) for i in range(1 << n))
        rep.add(e["id"], "exact", "integer reference against the closed form, " + unit, worst, REFERENCE_UNITS, worst <= REFERENCE_UNITS)
    if "grover" in e:
        m, k = e["grover"]
        p = ref.re[m] ** 2 + ref.im[m] ** 2
        err = abs(p - exact.grover_law(n, k)) >> exact.PREC
        rep.add(e["id"], "exact", "integer reference: P(marked) against sin²((2k+1)θ), computed in integers, " + unit, err, REFERENCE_UNITS, err <= REFERENCE_UNITS)
    if "qpe" in e:
        m, bits = e["qpe"]
        i = m | (1 << bits)
        err = abs(ref.re[i] ** 2 + ref.im[i] ** 2 - (1 << (2 * exact.PREC))) >> exact.PREC
        rep.add(e["id"], "exact", "integer reference: P(m) against 1, the closed form, " + unit, err, REFERENCE_UNITS, err <= REFERENCE_UNITS)
    if "qpe_third" in e:
        bits = e["qpe_third"]
        worst = 0
        for y in range(1 << bits):
            i = y | (1 << bits)
            worst = max(worst, abs(ref.re[i] ** 2 + ref.im[i] ** 2 - exact.qpe_third_law(bits, y)) >> exact.PREC)
        rep.add(e["id"], "exact", "integer reference: every outcome's probability against sin²(π2^tδ)/(2^2t sin²(πδ)), " + unit, worst, REFERENCE_UNITS, worst <= REFERENCE_UNITS)


def probe_states(e):
    """The circuit with its probed gate left out, and that gate's angle in radians."""
    steps = [s for i, s in enumerate(e["steps"]) if i != e["probe"]]
    return exact.run_unitary(e["n"], steps), abs(exact.float_radians(e["steps"][e["probe"]]["a"][0]))


def phase_held(rep, e, eng, ref, probe, re_ints, im_ints):
    """Whether an engine holds the probed gate or drops it: its distance alone cannot tell them
    apart once the gate's effect is below the engine's bound. Exactly, in integers, modulo a global
    phase: it holds the gate when its answer is within half the gate's effect of the ideal state, and
    drops it when within half of the state without the gate. When it is within half of neither, its
    own rounding is larger than the gate's effect, and the answer says nothing either way: that is
    recorded and not judged. The share of the step it took (exact.held_share) is recorded beside."""
    without, theta = probe
    share = exact.held_share(ref, without, re_ints, im_ints)
    step = exact.distance(ref.re, ref.im, without.re, without.im)[1]
    to_ideal = exact.distance(ref.re, ref.im, re_ints, im_ints)[1]
    to_without = exact.distance(without.re, without.im, re_ints, im_ints)[1]
    verdict = "holds it" if to_ideal <= step / 2 else "drops it" if to_without <= step / 2 else None
    q = "holds the phase of its small-angle gate: the share of the step from the state without that gate to the ideal (1 holds it, 0 drops it)"
    if verdict is None:
        c = rep.add(e["id"], eng, q, share, [0.5, 1.5], True, status="info", note="neither can be told: its distance from the ideal state, %.3g, is more than half the gate's effect, %.3g" % (to_ideal, step))
    elif eng == "ddsim.statevector" and theta < DDSIM_PARAMETER_TOLERANCE:
        c = rep.add(e["id"], eng, q, share, [0.5, 1.5], True, status="info", note="%s: the angle, %.3g rad, is below its parameter tolerance of %g rad, where it may read an angle as 0" % (verdict, theta, DDSIM_PARAMETER_TOLERANCE))
    else:
        c = rep.add(e["id"], eng, q, share, [0.5, 1.5], verdict == "holds it", note=verdict)
    c["step"] = step


def extension_figure(rep, e, fixed):
    """quantum-sim.md section 5: P2(19) applied 2^18 times to |1> ends at modulus 1.0000189 and
    within 10^-5 rad of the half turn. The modulus is stated to seven decimals."""
    re, im = int(fixed[1][0]), int(fixed[1][1])
    modulus = math.hypot(re, im) / (1 << 30)
    off = abs(math.remainder(math.atan2(im, re) - math.pi, 2 * math.pi))
    rep.add(e["id"], "wai.exact", "modulus of the half turn's amplitude against the extension's figure, 1.0000189 (section 5)", modulus, [1.00001885, 1.00001895], 1.00001885 <= modulus <= 1.00001895)
    rep.add(e["id"], "wai.exact", "its angle from the half turn, in radians, which the extension states is within 1e-5", off, 1e-5, off <= 1e-5)


def check_state_entry(e, runner, rep, engines):
    import numpy as np

    n = e["n"]
    ref = exact.run_unitary(n, e["steps"])
    want = ref_vector(ref)
    reference_checks(e, ref, rep)
    strict = bool(e.get("phase_is_result"))
    probe = probe_states(e) if "probe" in e else None
    exact_ops = None
    jobs = {}
    for eng in engines:
        if eng == "aer.qasm_import" and n > 14:
            continue
        if eng == "quimb.tensor_network" and n > 12:
            continue
        jobs[eng] = {"id": e["id"], "kind": "state", "n": n, "steps": e["steps"], "qasm": e["qasm"], "probe": e.get("probe")}
    answers = run_all(runner, jobs)
    quantity = "distance from the exact state" + (", with its phase" if strict else ", modulo a global phase")
    for eng, ans in answers.items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "state", fail_on_error=True)
            continue
        if not answered(rep, e["id"], eng):
            continue
        vec = load_state(ans["job"]["state_out"])
        if len(vec) != 1 << n:
            rep.add(e["id"], eng, "state length", len(vec), 1 << n, False)
            continue
        if eng == "wai.exact":
            fixed = np.load(ans["settings"]["fixed"])
            re = exact.ints_from_fixed([int(x) for x in fixed[:, 0]])
            im = exact.ints_from_fixed([int(x) for x in fixed[:, 1]])
            d, dm, phase, _, _ = exact.distance(ref.re, ref.im, re, im)
            bound, why = state_tolerance(eng, ans, e)
            measured = d if strict else dm
            c = rep.add(e["id"], eng, quantity, measured, bound, measured <= bound, note=why)
            c["strict"] = d
            c["global_phase"] = phase
            # The phase the method's result states: with it, the state is the reference's, which writes
            # each gate as method statevector does, nothing fitted.
            turn = complex(math.cos(ans.get("global_phase", 0.0)), math.sin(ans.get("global_phase", 0.0)))
            stated = exact.distance(ref.re, ref.im, *[exact.ints_from_floats([float(x) for x in part]) for part in ((vec * turn).real, (vec * turn).imag)])[0]
            rep.add(e["id"], eng, "distance from the exact state, with the global phase its result states", stated, bound + 8 * tol.U_DOUBLE, stated <= bound + 8 * tol.U_DOUBLE)
            exact_ops = ans["ops"]
            exact_answer(rep, e, n, ans, want)
            if e.get("extension_figure"):
                extension_figure(rep, e, fixed)
            continue
        bound, why = state_tolerance(eng, ans, e)
        d, dm, phase = exact_distance(ref, vec)
        measured = d if strict else dm
        c = rep.add(e["id"], eng, quantity, measured, bound, measured <= bound, note=why)
        c["strict"] = d
        c["global_phase"] = phase
        c["seconds"] = ans.get("seconds")
        if ans.get("notes"):
            c["adapter"] = ans["notes"]
        if probe:
            phase_held(rep, e, eng, ref, probe, exact.ints_from_floats([float(x) for x in vec.real]), exact.ints_from_floats([float(x) for x in vec.imag]))
    if is_clifford(e) and e["expect"] == "exact_state":
        check_stabilizer_entry(e, runner, rep)
    # Through the engine's whole path, as a user asks: chosen amplitudes, method auto.
    picks = sorted({0, (1 << n) - 1, 1 % (1 << n)} | {max(range(1 << n), key=lambda i: ref.re[i] ** 2 + ref.im[i] ** 2)})
    bits = [format(i, "0%db" % n) for i in picks]
    ans = runner.run("wai.auto", {"id": e["id"], "kind": "amplitudes", "n": n, "steps": e["steps"], "bitstrings": bits})
    if not ans.get("ok"):
        no_answer(rep, e["id"], "wai.auto", ans, "chosen amplitudes through the job path")
        return
    planned = ans.get("planned")
    # The bit-exact method lowers rz to a phase gate, and so on; its result states the phase that
    # drops, and its amplitudes are held to the reference's with that phase, nothing fitted.
    turn = 1
    if planned == "exact":
        stated_phase = ((ans.get("accuracy") or {}).get("global_phase") or {}).get("radians", 0.0)
        turn = complex(math.cos(stated_phase), math.sin(stated_phase))
    worst = 0.0
    for b in bits:
        a = complex(*ans["amplitudes"][b]) * turn
        worst = max(worst, exact.amplitude_error((ref.re[int(b, 2)], ref.im[int(b, 2)]), a.real, a.imag))
    if planned == "exact" and exact_ops is None:
        # The bound is the op-log's; the bit-exact method gives it.
        a2 = runner.run("wai.exact", {"id": e["id"], "kind": "state", "n": n, "steps": e["steps"]})
        exact_ops = a2.get("ops") if a2.get("ok") else None
    if planned == "exact" and exact_ops is not None:
        bound = tol.fixed_point_bound(n, exact_ops) + 4 * tol.U_DOUBLE
    else:
        bound = tol.double_state(tol.gate_count(e["steps"]), n, width=3)
    rep.add(e["id"], "wai.auto", "chosen amplitudes through the job path, worst error" + (", with the global phase the result states" if planned == "exact" else ""), worst, bound, worst <= bound, note="planned %s" % planned)


# ------------------------------------------------------------------ stabilizers


def check_stabilizer_entry(e, runner, rep):
    n = e["n"]
    tb = exact.Tableau(n)
    for s in e["steps"]:
        if "g" in s:
            tb.gate(s["g"], s["q"])
    rows, rank, x_rank = exact.reduced_group(tb.stabilizers(), n)
    texts = [exact.pauli_text(r, n) for r in rows]
    digest = hashlib.sha256("\n".join(texts).encode()).hexdigest()
    rep.add(e["id"], "exact", "the integer tableau's group has full rank", rank, n, rank == n)
    jobs = {"stim.tableau": {"id": e["id"], "kind": "stabilizers", "n": n, "steps": e["steps"]}, "wai.stabilizer": {"id": e["id"], "kind": "stabilizers", "n": n, "steps": e["steps"], "paulis": [t[1:] for t in texts]}}
    if n <= 1000:
        jobs["aer.stabilizer"] = {"id": e["id"], "kind": "stabilizers", "n": n, "steps": e["steps"]}
    answers = run_all(runner, jobs)
    for eng, ans in answers.items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "stabilizer group", fail_on_error=True)
            continue
        if not answered(rep, e["id"], eng):
            continue
        if eng == "wai.stabilizer":
            want = [-1 if t[0] == "-" else 1 for t in texts]
            wrong = sum(1 for a, b in zip(ans["values"], want) if a != b)
            c = rep.add(e["id"], eng, "generators of the reference group whose expectation is not their sign", wrong, 0, wrong == 0 and len(ans["values"]) == n)
            c["seconds"] = ans.get("seconds")
            rep.add(e["id"], eng, "outcome dimension against the rank of the group's X part", ans["outcome_dimension"], x_rank, ans["outcome_dimension"] == x_rank)
            continue
        theirs = [exact.row_from_text(s) for s in ans["stabilizers"]]
        r2, rank2, _ = exact.reduced_group(theirs, n)
        d2 = hashlib.sha256("\n".join(exact.pauli_text(r, n) for r in r2).encode()).hexdigest()
        c = rep.add(e["id"], eng, "reduced stabilizer group equals the reference's", d2 == digest, True, d2 == digest and rank2 == n)
        c["seconds"] = ans.get("seconds")


def fixed_bits(z_masks, n):
    """The bits whose value the affine space fixes: bit i is fixed when e_i lies in the span of
    the Z rows' masks (GF(2) elimination)."""
    basis = {}
    for m in z_masks:
        while m:
            top = m.bit_length() - 1
            if top not in basis:
                basis[top] = m
                break
            m ^= basis[top]
    out = set()
    for i in range(n):
        v = 1 << i
        while v:
            top = v.bit_length() - 1
            if top not in basis:
                break
            v ^= basis[top]
        if v == 0:
            out.add(i)
    return out


def check_stabilizer_counts(e, runner, rep, alpha):
    """Outcomes of a Clifford circuit measured in Z: uniform on an affine space of dimension x_rank.
    Every shot must lie in it; no outcome may come up more often than its probability 2^-x_rank
    allows (by Bernstein, over all 2^x_rank of them); and each bit the space leaves free must read
    1 half the time. A sampler confined to a smaller space that leaves every bit free is not seen at
    this width; the README says so."""
    n = e["n"]
    tb = exact.Tableau(n)
    gates = [s for s in e["steps"] if "g" in s]
    for s in gates:
        tb.gate(s["g"], s["q"])
    rows, _, x_rank = exact.reduced_group(tb.stabilizers(), n)
    z_rows = [r for r in rows if r[0] == 0]
    free = [i for i in range(n) if i not in fixed_bits([r[1] for r in z_rows], n)]
    shots = e["shots"]
    p = 2.0**-x_rank
    jobs = {eng: {"id": e["id"], "kind": "counts", "n": n, "bits": e["bits"], "steps": e["steps"], "shots": shots, "seed": 7} for eng in ("wai.stabilizer", "stim.tableau", "aer.stabilizer")}
    for eng, ans in run_all(runner, jobs).items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "outcomes inside the stabilizer's affine space")
            continue
        if not answered(rep, e["id"], eng):
            continue
        counts = ans["counts"]
        got = sum(counts.values())
        rep.add(e["id"], eng, "shots returned", got, shots, got == shots)
        outside = 0
        for k, v in counts.items():
            x = int(k, 2)
            for zr in z_rows:
                if ((zr[1] & x).bit_count() + zr[2]) % 2:
                    outside += v
                    break
        rep.add(e["id"], eng, "shots outside the affine space the stabilizers fix (dimension %d)" % x_rank, outside, 0, outside == 0)
        a = alpha / 3 / 2
        t = tol.bernstein(shots, p, a / 2.0**x_rank)
        most = max(counts.values(), default=0)
        rep.add(e["id"], eng, "the most frequent outcome's count above shots·2^-%d, in units of its Bernstein bound over all 2^%d outcomes" % (x_rank, x_rank), (most - shots * p) / t, 1.0, most - shots * p <= t)
        worst = 0.0
        bad = []
        for i in free:
            ones = sum(v for k, v in counts.items() if int(k, 2) >> i & 1)
            tb_ = tol.bernstein(shots, 0.5, a / max(len(free), 1))
            worst = max(worst, abs(ones - shots / 2) / tb_)
            if abs(ones - shots / 2) > tb_:
                bad.append("bit %d: %d ones of %d" % (i, ones, shots))
        rep.add(e["id"], eng, "each of the %d bits the space leaves free against 1/2, largest deviation in units of its Bernstein bound" % len(free), worst, 1.0, not bad, note="; ".join(bad[:4]) if bad else None)
        rep.add(e["id"], eng, "distinct outcomes, of %d shots over 2^%d" % (shots, x_rank), len(counts), None, True, status="info")


# ------------------------------------------------------------------ sampled distributions


def check_counts_entry(e, runner, rep, plan):
    dist = exact.distribution(e["n"], e["bits"], e["steps"])
    engines = list(COUNT_ENGINES)
    if is_clifford(e):
        engines += ["wai.stabilizer", "stim.tableau"]
    engines = e.get("count_engines") or engines
    jobs = {eng: {"id": e["id"], "kind": "counts", "n": e["n"], "bits": e["bits"], "steps": e["steps"], "shots": e["shots"], "seed": 11} for eng in engines}
    for eng, ans in run_all(runner, jobs).items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "counts", fail_on_error=True)
            continue
        if not answered(rep, e["id"], eng):
            continue
        plan.append((e, eng, ans, dist))


def judge_counts(plan, rep, alpha):
    """Every outcome's count against its exact probability, by Bernstein's inequality, and the
    whole distribution at once against the exact one, in L1 distance (tolerance.py); each test
    at a share of `alpha`, so that a correct engine fails them with probability at most alpha."""
    tests = sum(len(set(d) | set(a["counts"])) + 1 for _, _, a, d in plan)
    a = alpha / max(tests, 1)
    for e, eng, ans, dist in plan:
        shots = sum(ans["counts"].values())
        if shots != e["shots"]:
            rep.add(e["id"], eng, "shots returned", shots, e["shots"], False)
            continue
        worst = 0.0
        bad = []
        for k in sorted(set(dist) | set(ans["counts"])):
            p = float(dist.get(k, 0))
            x = ans["counts"].get(k, 0)
            if p == 0:
                if x:
                    bad.append("%s: %d shots, probability 0" % (k, x))
                continue
            t = tol.bernstein(shots, p, a)
            dev = abs(x - shots * p)
            worst = max(worst, dev / t)
            if dev > t:
                bad.append("%s: %d shots, expected %.1f ± %.1f" % (k, x, shots * p, t))
        c = rep.add(e["id"], eng, "largest deviation of an outcome's count, in units of its Bernstein bound", worst, 1.0, not bad, note="; ".join(bad[:4]) if bad else None)
        if ans.get("planned"):
            c["planned"] = ans["planned"]
        support = sum(1 for v in dist.values() if v > 0)
        l1 = math.fsum(abs(ans["counts"].get(k, 0) / shots - float(dist.get(k, 0))) for k in set(dist) | set(ans["counts"]))
        eps = tol.l1_distance(shots, support + 1, a)
        rep.add(e["id"], eng, "L1 distance from the exact distribution, Σ|count/shots − p| over every outcome (twice the total variation), against its bound", l1, eps, l1 <= eps)


# ------------------------------------------------------------------ wider circuits


def check_mirror_entry(e, runner, rep):
    import numpy as np

    n = e["n"]
    engines = ["wai.exact", "wai.statevector", "aer.statevector", "quimb.dense", "ddsim.statevector", "qsim.statevector"]
    jobs = {eng: {"id": e["id"], "kind": "state", "n": n, "steps": e["steps"]} for eng in engines}
    for eng, ans in run_all(runner, jobs).items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "distance from |0…0>")
            continue
        if not answered(rep, e["id"], eng):
            continue
        if eng == "wai.exact":
            fixed = np.load(ans["settings"]["fixed"])
            v = fixed[:, 0].astype(np.float64) / (1 << 30) + 1j * fixed[:, 1].astype(np.float64) / (1 << 30)
        else:
            v = load_state(ans["job"]["state_out"])
        target = np.zeros_like(v)
        target[0] = 1
        d, dm = float_distance(v, target)
        bound, why = state_tolerance(eng, ans, e)
        c = rep.add(e["id"], eng, "distance from |0…0>, modulo a global phase", dm, bound, dm <= bound, note=why)
        c["strict"] = d
        c["seconds"] = ans.get("seconds")
        if eng == "wai.exact":
            exact_answer(rep, e, n, ans, target)


def check_pairwise_entry(e, runner, rep):
    import numpy as np

    n = e["n"]
    engines = ["wai.statevector", "aer.statevector", "quimb.dense", "ddsim.statevector", "qsim.statevector"]
    if corpus.reader_class(e)[0] == "exact":
        engines.insert(0, "wai.exact")
    jobs = {eng: {"id": e["id"], "kind": "state", "n": n, "steps": e["steps"]} for eng in engines}
    got = {}
    for eng, ans in run_all(runner, jobs).items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "state, against the other engines")
            continue
        if not answered(rep, e["id"], eng):
            continue
        if eng == "wai.exact":
            fixed = np.load(ans["settings"]["fixed"])
            v = fixed[:, 0].astype(np.float64) / (1 << 30) + 1j * fixed[:, 1].astype(np.float64) / (1 << 30)
        else:
            v = load_state(ans["job"]["state_out"])
        got[eng] = (v, state_tolerance(eng, ans, e)[0], ans)
    names = sorted(got)
    for i, a in enumerate(names):
        for b in names[i + 1:]:
            _, dm = float_distance(got[a][0], got[b][0])
            bound = got[a][1] + got[b][1]
            rep.add(e["id"], "%s ~ %s" % (a, b), "distance between the two states, modulo a global phase", dm, bound, dm <= bound)
    if "qft_x" in e:
        r = corpus.Stream("pick/" + e["id"])
        ys = [r.below(1 << n) for _ in range(64)]
        for eng in names:
            v = got[eng][0]
            worst = max(exact.amplitude_error(exact.qft_amplitude(n, e["qft_x"], y), float(v[y].real), float(v[y].imag)) for y in ys)
            rep.add(e["id"], eng, "64 amplitudes against the closed form, worst error, modulo nothing", worst, got[eng][1], worst <= got[eng][1])


def wide_amplitudes(e):
    n = e["n"]
    r = corpus.Stream("wide/" + e["id"])
    if e.get("wide_kind") == "ghz":
        bits = ["0" * n, "1" * n, "".join(str(r.below(2)) for _ in range(n))]
        ref = {bits[0]: (exact.rsqrt2(), 0), bits[1]: (exact.rsqrt2(), 0), bits[2]: (0, 0)}
        return bits, ref
    if e.get("wide_kind") == "qft":
        ys = [r.below(1 << n) for _ in range(6)]
        bits = [format(y, "0%db" % n) for y in ys]
        return bits, {b: exact.qft_amplitude(n, e["qft_x"], int(b, 2)) for b in bits}
    if e["expect"] == "wide_mirror":
        b = "0" * n
        return [b], {b: (exact.ONE, 0)}
    bits = ["".join(str(r.below(2)) for _ in range(n)) for _ in range(8)]
    return bits, None


def mps_bound(eng, ans, e):
    """tolerance.mps for each engine's matrix product state, charged what that engine reports it
    truncated: the engine under test its summed discarded weight over its truncations, the
    simulator with a log its every truncation's weight and its two-site updates, and the library
    with an absolute cutoff the most that cutoff drops at each update (engines.py)."""
    steps = tol.svd_steps(e["steps"])
    acc = ans.get("accuracy") or {}
    if eng == "wai.mps":
        chi = acc.get("max_bond_used") or acc.get("max_bond") or 64
        trunc = tol.truncation_from(total=acc.get("discarded_weight", 0.0), count=acc.get("truncations", 0))
        why = "%d two-site updates at bond %d; truncation from its reported discarded weight" % (steps, chi)
    elif eng == "aer.mps":
        log = ans["settings"]["mps_log"]
        steps, chi = log["updates"], log["max_bond"]
        trunc = tol.truncation_from(weights=log["discarded"])
        why = "%d two-site updates at bond %d, and %d truncations, from its log" % (steps, chi, len(log["discarded"]))
    else:
        s = ans["settings"]
        chi = s.get("max_bond_used") or 1
        trunc = 2 * steps * math.sqrt(s["cutoff"])
        why = "%d two-site updates at bond %d; its cutoff drops at most %g of weight at each" % (steps, chi, s["cutoff"])
    return tol.mps(steps, chi, trunc), why


def check_wide_entry(e, runner, rep):
    n = e["n"]
    bits, ref = wide_amplitudes(e)
    engines = ["wai.mps", "wai.tensor_network", "aer.mps", "quimb.mps", "quimb.tensor_network"]
    if n <= 26:
        engines += ["wai.statevector", "aer.statevector"]
    engines = e.get("only_engines") or engines
    jobs = {eng: {"id": e["id"], "kind": "amplitudes", "n": n, "steps": e["steps"], "bitstrings": bits, "memory_gb": 24, "max_bond": e.get("max_bond")} for eng in engines}
    got = {}
    gates = tol.gate_count(e["steps"])
    smallest = min((math.hypot(*ref[b]) / exact.ONE for b in bits if ref and any(ref[b])), default=None) if ref else None
    for eng, ans in run_all(runner, jobs).items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "amplitudes")
            continue
        if not answered(rep, e["id"], eng):
            continue
        if eng in ("wai.mps", "aer.mps", "quimb.mps"):
            bound, why = mps_bound(eng, ans, e)
        elif eng in ("wai.tensor_network", "quimb.tensor_network"):
            if "contraction_flops" in ans:
                w = math.log2(max(ans["contraction_flops"], 2.0))
            else:
                w = float(ans.get("settings", {}).get("contraction_width") or n)
            bound = tol.contraction(gates + n, min(w, 40.0))
            why = "contraction of width %.1f" % w
        else:
            bound = tol.double_state(gates, n, width=5)
            why = "double precision state vector"
        got[eng] = (ans["amplitudes"], bound, ans)
        if ref is not None:
            worst = max(exact.amplitude_error(ref[b], *ans["amplitudes"][b]) for b in bits)
            # A bound as large as half the smallest amplitude checked cannot see a wrong phase on it.
            if smallest is not None and bound >= smallest / 2:
                c = rep.add(e["id"], eng, "amplitudes against the closed form, worst error", worst, bound, True, status="info", note=why + "; this bound is not below half the amplitudes it checks (%.3g), so it cannot tell a wrong phase" % smallest)
            else:
                c = rep.add(e["id"], eng, "amplitudes against the closed form, worst error", worst, bound, worst <= bound, note=why)
            c["seconds"] = ans.get("seconds")
            if ans.get("accuracy"):
                c["accuracy"] = ans["accuracy"]
    if ref is None:
        names = sorted(got)
        for i, a in enumerate(names):
            for b in names[i + 1:]:
                worst = max(abs(complex(*got[a][0][x]) - complex(*got[b][0][x])) for x in bits)
                bound = got[a][1] + got[b][1]
                rep.add(e["id"], "%s ~ %s" % (a, b), "chosen amplitudes, worst difference", worst, bound, worst <= bound)


def check_expectation_entry(e, runner, rep):
    kind, p, k = e["closed"]
    p = Fraction(p)
    if kind == "depol1":
        ref = (-1) ** k * (1 - 4 * p / 3) ** k
    else:
        ref = (1 - 16 * p / 15) ** k
    jobs = {eng: {"id": e["id"], "kind": "expectation", "n": e["n"], "steps": e["steps"], "noise": e["noise"], "observable": e["observable"]} for eng in ("wai.density_matrix", "aer.density_matrix")}
    for eng, ans in run_all(runner, jobs).items():
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "expectation under noise")
            continue
        if not answered(rep, e["id"], eng):
            continue
        d = abs(Fraction(ans["value"]) - ref)
        bound = tol.density(k, e["n"])
        rep.add(e["id"], eng, "<%s> against the closed form" % e["observable"], float(d), bound, float(d) <= bound)


# ------------------------------------------------------------------ error correction


def qec_text(q, runner, rep):
    """The committed circuit, which the pinned generator must still write byte for byte."""
    text = open(os.path.join(ROOT, "corpus", "qec", q["id"] + ".circuit")).read()
    ans = runner.run("stim.generate", dict(q, kind="generate"))
    same = ans.get("ok") and ans.get("text") == text
    rep.add(q["id"], "stim.generate", "the pinned generator writes the committed circuit, byte for byte", bool(same), True, bool(same), note=None if same else (ans.get("error") or "the text differs"))
    return text


def merged(mechs):
    out = {}
    for p, ds, ls in mechs:
        key = (tuple(ds), tuple(ls))
        q = out.get(key, 0.0)
        out[key] = q + p - 2 * q * p
    return out


def check_qec_entry(q, text, runner, rep, qec_plan):
    jobs = {eng: {"id": q["id"], "kind": "qec_model", "text": text} for eng in ("wai.qec", "stim.tableau")}
    models = run_all(runner, jobs)
    a, b = models["wai.qec"], models["stim.tableau"]
    if not a.get("ok"):
        no_answer(rep, q["id"], "wai.qec", a, "detector error model")
        if not b.get("ok"):
            rep.not_run(q["id"], "stim.tableau", "detector error model", str(words(b))[:300])
        return
    if not answered(rep, q["id"], "wai.qec"):
        return
    if not b.get("ok"):
        # A channel the other engine gives no exact model of (a general Pauli channel): its
        # sampler still holds the engine under test's model, by the rates below.
        rep.not_run(q["id"], "stim.tableau", "detector error model", "gives no model of this circuit, which the engine under test models: %s" % str(words(b))[:300])
        shots = 200000 if q["tier"] == "ci" else 1000000
        sj = {eng: {"id": q["id"], "kind": "qec_sample", "text": text, "shots": shots, "seed": 5} for eng in ("wai.qec", "stim.tableau")}
        qec_plan.append((q, merged(a["mechanisms"]), a["detectors"], a["observables"], run_all(runner, sj), shots))
        return
    rep.add(q["id"], "wai.qec ~ stim.tableau", "detectors and observables", [a["detectors"], a["observables"]], [b["detectors"], b["observables"]], [a["detectors"], a["observables"]] == [b["detectors"], b["observables"]])
    ma, mb = merged(a["mechanisms"]), merged(b["mechanisms"])
    only_a = sorted(set(ma) - set(mb))
    only_b = sorted(set(mb) - set(ma))
    sites = text_noise_sites(text)
    bound = tol.dem(sites)
    worst = max([abs(ma[k] - mb[k]) for k in set(ma) & set(mb)] or [0.0])
    ok = not only_a and not only_b and worst <= bound
    note = None
    if only_a or only_b:
        note = "symptom sets only in the engine under test: %d, only in the other: %d (first: %s / %s)" % (len(only_a), len(only_b), only_a[:1], only_b[:1])
    c = rep.add(q["id"], "wai.qec ~ stim.tableau", "mechanisms: same symptom sets, probabilities within the bound (largest difference)", worst, bound, ok, note=note)
    c["mechanisms"] = [len(ma), len(mb)]
    c["seconds"] = [a.get("seconds"), b.get("seconds")]
    shots = 200000 if q["tier"] == "ci" else 1000000
    sj = {eng: {"id": q["id"], "kind": "qec_sample", "text": text, "shots": shots, "seed": 5} for eng in ("wai.qec", "stim.tableau")}
    samples = run_all(runner, sj)
    qec_plan.append((q, ma, a["detectors"], a["observables"], samples, shots))


def text_noise_sites(text):
    """How many noise channels a circuit applies, its blocks unrolled: an upper bound on how many
    channels one mechanism can merge."""
    import re

    total = 0
    stack = [1]
    for line in text.splitlines():
        s = line.strip()
        if s.startswith("REPEAT"):
            stack.append(stack[-1] * int(s.split()[1]))
            continue
        if s == "}":
            stack.pop()
            continue
        m = re.match(r"([A-Z_0-9]+)(\[[^\]]*\])?\(([^)]*)\)\s*(.*)", s)
        if m and (m.group(1).endswith("ERROR") or m.group(1).startswith(("DEPOLARIZE", "PAULI_CHANNEL")) or m.group(1) in ("E", "M", "MX", "MR", "MRX", "MY")):
            total += stack[-1] * max(len(m.group(4).split()), 1)
    return total


def judge_qec(qec_plan, rep, alpha):
    """Each detector's and observable's firing count against the rate the model gives,
    by Bernstein's inequality at a share of `alpha`."""
    tests = sum(d + o for _, _, d, o, _, _ in qec_plan) * 2
    a = alpha / max(tests, 1)
    for q, model, nd, no, samples, shots in qec_plan:
        # The rate at which a detector fires: an odd number of its mechanisms, each independent.
        prod_d = [1.0] * nd
        prod_o = [1.0] * no
        for (ds, ls), p in model.items():
            for d in ds:
                prod_d[d] *= 1 - 2 * p
            for o in ls:
                prod_o[o] *= 1 - 2 * p
        rates = [(1 - x) / 2 for x in prod_d] + [(1 - x) / 2 for x in prod_o]
        for eng, ans in samples.items():
            if not ans.get("ok"):
                no_answer(rep, q["id"], eng, ans, "detection rates")
                continue
            shape = [len(ans.get("detectors") or []), len(ans.get("observables") or []), ans.get("shots")]
            if shape != [nd, no, shots]:
                rep.add(q["id"], eng, "detectors, observables and shots returned", shape, [nd, no, shots], False)
                continue
            counts = ans["detectors"] + ans["observables"]
            worst = 0.0
            bad = []
            for i, (x, r) in enumerate(zip(counts, rates)):
                if r == 0:
                    if x:
                        bad.append("%d fired %d times at rate 0" % (i, x))
                    continue
                t = tol.bernstein(shots, r, a)
                dev = abs(x - shots * r)
                worst = max(worst, dev / t)
                if dev > t:
                    bad.append("%s %d: %d, expected %.1f ± %.1f" % ("D" if i < nd else "L", i if i < nd else i - nd, x, shots * r, t))
            c = rep.add(q["id"], eng, "detection and flip rates against the model, largest deviation in units of its Bernstein bound", worst, 1.0, not bad, note="; ".join(bad[:3]) if bad else None)
            c["seconds"] = ans.get("seconds")


# ------------------------------------------------------------------ the engine's other methods


class Estimates:
    """The run's statistical tests of estimates (expectation values from trajectories and Pauli
    frames, mitigated values from shots, the sampled stabilizer-rank sum's stated confidence), each
    judged once the run knows how many there are, at an equal share of the kind's family-wise rate.
    A test is (how many tests it is, a function of each test's alpha that writes its checks)."""

    def __init__(self):
        self.items = []

    def add(self, count, judge):
        self.items.append((count, judge))

    def judge(self, alpha):
        total = sum(c for c, _ in self.items)
        for _, fn in self.items:
            fn(alpha / max(total, 1))


def gates_of(e):
    return [s for s in e["steps"] if "g" in s]


def measured_of(e):
    return [(s["measure"], s["bit"]) for s in e["steps"] if "measure" in s]


def pauli_rotations(steps):
    total = 0
    for s in steps:
        if s["g"] not in corpus.PAULI_ROTATIONS:
            raise ValueError("no rotation count for %s" % s["g"])
        total += corpus.PAULI_ROTATIONS[s["g"]]
    return total


def check_pauli_entry(e, runner, rep):
    """Pauli propagation's expectation values against the integer state's: at each threshold of
    each ladder, |value − exact| at most the truncation bound the method states, plus its rounding
    (tolerance.pauli_rounding)."""
    n = e["n"]
    ref = exact.run_unitary(n, e["steps"])
    want = {o: exact.pure_expectation(ref, o) for o in e["observables"]}
    rotations = pauli_rotations(gates_of(e))
    for th in e["thresholds"]:
        ans = runner.run("wai.pauli", {"id": e["id"], "kind": "expectations", "n": n, "steps": e["steps"], "observables": e["observables"], "options": {"threshold": th, "ladder": True}})
        if not ans.get("ok"):
            no_answer(rep, e["id"], "wai.pauli", ans, "expectation values by Pauli propagation", fail_on_error=True)
            continue
        if not answered(rep, e["id"], "wai.pauli"):
            continue
        rows = ans["result"]["expectations"]
        for o, x in zip(e["observables"], rows):
            for r in x.get("ladder") or [dict(x, threshold=th)]:
                b, terms = r["bound"], r["peak_terms"]
                err = abs(Fraction(r["value"]) - want[o])
                allowed = b + tol.pauli_bound_rounding(b, terms, rotations) + tol.pauli_rounding(rotations, terms, x["coeff"])
                c = rep.add(e["id"], "wai.pauli", "<%s> within its truncation bound and rounding of the integer state's, threshold %g" % (o, r["threshold"]), float(err), allowed, err <= Fraction(allowed), note="bound %.3g, %d terms at the peak, at most %d rotations" % (b, terms, rotations))
                c["seconds"] = ans.get("seconds")
            if th == 0:
                rep.add(e["id"], "wai.pauli", "<%s>: the truncation bound at threshold 0, where only exact zeros drop" % o, x["bound"], 0.0, x["bound"] == 0.0)


def chosen_bits(e, ref, picks):
    """Every basis state, or `picks` of them drawn from the entry's stream, with the most likely."""
    n = e["n"]
    if picks is None:
        return [format(i, "0%db" % n) for i in range(1 << n)]
    r = corpus.Stream("amplitudes/" + e["id"])
    xs = {r.below(1 << n) for _ in range(picks)}
    xs.add(max(range(1 << n), key=lambda i: ref.re[i] ** 2 + ref.im[i] ** 2))
    return [format(i, "0%db" % n) for i in sorted(xs)]


AMPLITUDES_PER_JOB = 64  # the most amplitudes one job may ask for


def run_amplitudes(runner, eng, job, bits):
    """A job's amplitudes asked for in as many jobs as the job service's limit needs: one answer,
    the first that is not ok, or the amplitudes of them all with the first's other fields."""
    first = None
    for k in range(0, len(bits), AMPLITUDES_PER_JOB):
        a = runner.run(eng, dict(job, kind="amplitudes", bitstrings=bits[k:k + AMPLITUDES_PER_JOB]))
        if not a.get("ok"):
            return a
        if first is None:
            first = a
        else:
            first["amplitudes"].update(a["amplitudes"])
            first["seconds"] = (first.get("seconds") or 0) + (a.get("seconds") or 0)
    return first


def amplitudes_distance(ref, amps, bits):
    """(the distance over the chosen amplitudes modulo a global phase, the strict distance), from the
    integer reference: the phase that best aligns them, then the exact differences."""
    ip = 0j
    for b in bits:
        i = int(b, 2)
        a = complex(*amps[b])
        ip += a.conjugate() * complex(float(Fraction(ref.re[i], exact.ONE)), float(Fraction(ref.im[i], exact.ONE)))
    turn = ip / abs(ip) if abs(ip) > 0 else 1.0
    d2 = s2 = 0.0
    for b in bits:
        i = int(b, 2)
        a = complex(*amps[b])
        r = (ref.re[i], ref.im[i])
        d2 += exact.amplitude_error(r, (a * turn).real, (a * turn).imag) ** 2
        s2 += exact.amplitude_error(r, a.real, a.imag) ** 2
    return math.sqrt(d2), math.sqrt(s2)


def check_stabrank_entry(e, runner, rep, est):
    """The stabilizer-rank sum. Exactly: every amplitude, and each expectation value through the
    sandwich of twice the T gates, within the rounding of the sum's one evaluation
    (tolerance.stabrank_exact). Sampled: each draw states a bound that holds with a stated
    confidence; over the seeds the bound may fail no more often than a binomial count at that
    confidence allows (tolerance.binomial_upper)."""
    n = e["n"]
    ref = exact.run_unitary(n, e["steps"])
    bits = chosen_bits(e, ref, None if n <= 8 else 64)
    ans = run_amplitudes(runner, "wai.stabrank", {"id": e["id"], "n": n, "steps": e["steps"]}, bits)
    if not ans.get("ok"):
        no_answer(rep, e["id"], "wai.stabrank", ans, "amplitudes of the exact stabilizer-rank sum", fail_on_error=True)
        return
    if not answered(rep, e["id"], "wai.stabrank"):
        return
    acc = ans.get("accuracy") or {}
    t = acc.get("t_count", 0)
    bound = tol.stabrank_exact(t, n)
    worst = max(exact.amplitude_error((ref.re[int(b, 2)], ref.im[int(b, 2)]), *ans["amplitudes"][b]) for b in bits)
    c = rep.add(e["id"], "wai.stabrank", "amplitudes of the exact sum against the integer state, worst error, with the circuit's phase (%d T gates)" % t, worst, bound, worst <= bound and acc.get("class") == "exact", note="class %s, %s terms" % (acc.get("class"), acc.get("terms")))
    c["seconds"] = ans.get("seconds")
    if e.get("observables"):
        a2 = runner.run("wai.stabrank", {"id": e["id"], "kind": "expectations", "n": n, "steps": e["steps"], "observables": e["observables"]})
        if not a2.get("ok"):
            no_answer(rep, e["id"], "wai.stabrank", a2, "expectation values of the exact stabilizer-rank sum", fail_on_error=True)
        else:
            for o, x in zip(e["observables"], a2["result"]["expectations"]):
                want = exact.pure_expectation(ref, o)
                b = tol.stabrank_exact(2 * t, n, x["coeff"]) + abs(x["value"]) * tol.U_DOUBLE
                err = abs(Fraction(x["value"]) - want)
                rep.add(e["id"], "wai.stabrank", "<%s> of the exact sum (the sandwich of %d T gates) against the integer state's" % (o, 2 * t), float(err), b, err <= Fraction(b))
    sp = e.get("sampled")
    if not sp:
        return
    draws = []
    for seed in range(1, sp["seeds"] + 1):
        a = run_amplitudes(runner, "wai.stabrank", {"id": e["id"], "n": n, "steps": e["steps"], "seed": seed, "options": {"terms": sp["terms"], "confidence": sp["confidence"]}}, bits)
        if not a.get("ok"):
            no_answer(rep, e["id"], "wai.stabrank", a, "amplitudes of the sampled stabilizer-rank sum", fail_on_error=True)
            return
        draws.append(a)

    def judge(alpha):
        delta = 1 - sp["confidence"]
        misses = []
        for seed, a in enumerate(draws, 1):
            eps = a["accuracy"]["bound"]
            # The bound holds for the state; its evaluation adds the sum's own rounding.
            room = eps + tol.stabrank_exact(t, n)
            off = max(exact.amplitude_error((ref.re[int(b, 2)], ref.im[int(b, 2)]), *a["amplitudes"][b]) for b in bits)
            if off > room:
                misses.append("seed %d: %.3g past %.3g" % (seed, off, room))
        most = tol.binomial_upper(len(draws), delta, alpha)
        rep.add(e["id"], "wai.stabrank", "sampled draws whose stated bound failed (%d terms, confidence %g), against the most a binomial count at that confidence allows" % (sp["terms"], sp["confidence"]), len(misses), most - 1, len(misses) < most, note="; ".join(misses[:4]) or None)

    est.add(1, judge)


def check_cutrank_entry(e, runner, rep):
    """The cut-rank sum's amplitudes against the integer state's, over every amplitude or a seeded
    choice, modulo a global phase in the 2-norm, within √m times its per-amplitude rounding
    (tolerance.cutrank); the distance with the phase it states is recorded beside."""
    n = e["n"]
    ref = exact.run_unitary(n, e["steps"])
    bits = chosen_bits(e, ref, e.get("picks"))
    ans = run_amplitudes(runner, "wai.cutrank", {"id": e["id"], "n": n, "steps": e["steps"]}, bits)
    if not ans.get("ok"):
        no_answer(rep, e["id"], "wai.cutrank", ans, "amplitudes by the cut-rank sum", fail_on_error=True)
        return
    if not answered(rep, e["id"], "wai.cutrank"):
        return
    dm, d = amplitudes_distance(ref, ans["amplitudes"], bits)
    bound = math.sqrt(len(bits)) * tol.cutrank(n)
    acc = ans.get("accuracy") or {}
    c = rep.add(e["id"], "wai.cutrank", "%d amplitudes against the integer state, 2-norm modulo a global phase" % len(bits), dm, bound, dm <= bound and acc.get("class") == "exact", note="class %s, widest cut %s" % (acc.get("class"), acc.get("max_cut")))
    c["strict"] = d
    c["seconds"] = ans.get("seconds")


def check_noisy_entry(e, runner, rep, plan, est):
    """A noisy program against the integer density matrix (exact.noisy_density): the density
    matrix's expectation values within its rounding (tolerance.density_steps), the estimates of
    trajectories and Pauli frames within Hoeffding's bound on their samples, and every method's
    counts against the exact law of the bits, readout flips included (judge_counts)."""
    n = e["n"]
    gates = gates_of(e)
    rho = exact.noisy_density(n, gates, e["noise"])
    want = {o: rho.expectation(o) for o in e["observables"]}
    readout = exact.noise_of(e["noise"])["readout"]
    dist = exact.noisy_distribution(rho, measured_of(e), e["bits"], readout)
    steps = tol.noisy_steps([len(g["q"]) for g in gates], e["noise"])
    for eng in e["engines"]:
        if eng == "aer.density_matrix":
            for o in e["observables"]:
                a = runner.run(eng, {"id": e["id"], "kind": "expectation", "n": n, "steps": gates, "noise": e["noise"], "observable": o})
                if not a.get("ok"):
                    no_answer(rep, e["id"], eng, a, "<%s> under noise" % o)
                    continue
                b = tol.density(len(gates), n)
                err = abs(Fraction(a["value"]) - want[o])
                rep.add(e["id"], eng, "<%s> under noise against the integer density matrix's" % o, float(err), b, err <= Fraction(b))
            continue
        ans = runner.run(eng, {"id": e["id"], "kind": "expectations", "n": n, "bits": e["bits"], "steps": e["steps"], "noise": e["noise"], "observables": e["observables"], "shots": e["shots"], "seed": 13})
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "counts and expectation values under noise", fail_on_error=True)
            continue
        if not answered(rep, e["id"], eng):
            continue
        r = ans["result"]
        acc = r.get("accuracy") or {}
        plan.append((e, eng, {"counts": r.get("counts", {}), "planned": ans.get("planned")}, dist))
        rows = r.get("expectations", [])
        if eng == "wai.density_matrix":
            for o, x in zip(e["observables"], rows):
                b = abs(x["coeff"]) * tol.density_steps(n, steps)
                err = abs(Fraction(x["value"]) - want[o])
                rep.add(e["id"], eng, "<%s> under noise against the integer density matrix's" % o, float(err), b, err <= Fraction(b), note="%d superoperators" % steps)
            continue
        samples = acc.get("trajectories") if eng == "wai.trajectories" else acc.get("shots")
        # Each trajectory's state is a double-precision run of the gates and the Paulis its errors
        # insert (at most one for each place a channel acts); each frame's value is exact.
        state_round = 2 * tol.double_state(len(gates) + steps, n, width=5) if eng == "wai.trajectories" else 0.0

        def judge(alpha, eng=eng, rows=rows, samples=samples, state_round=state_round):
            for o, x in zip(e["observables"], rows):
                spread = 2 * abs(x["coeff"])
                b = tol.hoeffding(samples, spread, alpha) + abs(x["coeff"]) * state_round
                err = abs(Fraction(x["value"]) - want[o])
                c = rep.add(e["id"], eng, "<%s> under noise, the mean of %d samples, against the integer density matrix's within Hoeffding's bound" % (o, samples), float(err), b, err <= Fraction(b), note="stated standard error %.3g" % x.get("stderr", float("nan")))
                c["stderr"] = x.get("stderr")

        est.add(len(rows), judge)


def richardson(factors):
    """The coefficients that carry values at `factors` to zero by the polynomial through them."""
    fs = [Fraction(f) for f in factors]
    out = []
    for i, fi in enumerate(fs):
        c = Fraction(1)
        for j, fj in enumerate(fs):
            if j != i:
                c *= (0 - fj) / (fi - fj)
        out.append(c)
    return out


def check_mitigation_entry(e, runner, rep, est):
    """Mitigation, against what each estimate tends to, computed exactly from the integer density
    matrix: zne's value at each scale (the folded program's noisy value) and its extrapolation
    (the same extrapolation of those exact values); vd's purity Tr(ρ²) and its ratio
    Tr(ρ²O)/Tr(ρ²); readout's quasi-probabilities (the law of the bits before their flips). Each
    value from shots is held within Hoeffding's bound on them; the ideal value each result states
    within the double-precision state's rounding."""
    n = e["n"]
    gates = gates_of(e)
    m = e["mitigate"]
    method = m["method"]
    job = {"id": e["id"], "kind": "mitigated", "n": n, "bits": e["bits"], "steps": e["steps"], "noise": e["noise"], "observables": e["observables"], "mitigate": m, "shots": e.get("shots", 0), "seed": 17, "method": "density_matrix"}
    ans = runner.run("wai.mitigate", job)
    if not ans.get("ok"):
        no_answer(rep, e["id"], "wai.mitigate", ans, "mitigation (%s)" % method, fail_on_error=True)
        return
    if not answered(rep, e["id"], "wai.mitigate"):
        return
    res = ans["result"]["mitigation"]
    rows = res.get("observables", [])
    ideal_state = exact.run_unitary(n, gates)
    for o, x in zip(e["observables"], rows):
        want = exact.pure_expectation(ideal_state, o)
        b = 2 * tol.double_state(len(gates), n, width=5)
        err = abs(Fraction(x["ideal"]) - want)
        rep.add(e["id"], "wai.mitigate", "%s: the ideal <%s> it states, against the integer state's" % (method, o), float(err), b, err <= Fraction(b))
    if method == "zne":
        shots = m["shots"]
        coeffs = richardson(m["factors"])
        scales = res["details"]["scales"]
        exact_at = []
        for f in m["factors"]:
            fold = (int(f) - 1) // 2
            rho = exact.noisy_density(n, gates, e["noise"], fold=fold)
            g = len(gates) * (2 * fold + 1)
            exact_at.append(({o: rho.expectation(o) for o in e["observables"]}, tol.density_steps(n, tol.noisy_steps([len(x["q"]) for x in gates] * (2 * fold + 1), e["noise"])), g))

        def judge(alpha):
            t = tol.hoeffding(shots, 2.0, alpha)
            for i, o in enumerate(e["observables"]):
                for k, f in enumerate(m["factors"]):
                    v = scales[k]["values"][i]["value"]
                    x, dm, _ = exact_at[k]
                    err = abs(Fraction(v) - x[o])
                    rep.add(e["id"], "wai.mitigate", "zne: <%s> at scale %g, from %d shots, against the folded program's exact noisy value" % (o, f, shots), float(err), t + dm, err <= Fraction(t + dm))
                want = sum(c * exact_at[k][0][o] for k, c in enumerate(coeffs))
                room = sum(abs(float(c)) * (t + exact_at[k][1]) for k, c in enumerate(coeffs))
                got = rows[i]["mitigated"]
                err = abs(Fraction(got) - want)
                rep.add(e["id"], "wai.mitigate", "zne: <%s> extrapolated, against the Richardson extrapolation of the exact noisy values at its scales" % o, float(err), room, err <= Fraction(room), note="exact extrapolation %.6g, ideal %.6g" % (float(want), rows[i]["ideal"]))

        est.add(len(e["observables"]) * len(m["factors"]), judge)
    elif method == "vd":
        shots = m["shots"]
        rho = exact.noisy_density(n, gates, e["noise"])
        purity = rho.purity()
        eps_f = tol.density_steps(n, tol.noisy_steps([len(x["q"]) for x in gates], e["noise"])) / math.sqrt(2**n)
        d = 1 << n
        pb = 2 * eps_f + 4 * d * tol.U_DOUBLE
        err = abs(Fraction(res["details"]["purity"]) - purity)
        rep.add(e["id"], "wai.mitigate", "vd: the purity Tr(ρ²) it states, against the integer density matrix's", float(err), pb, err <= Fraction(pb))

        # Tr(ρ²P) from the engine's ρ: within 2·ε_F·√d of the exact one (|Tr(AP)| ≤ ‖A‖_F‖P‖_F,
        # ‖P‖_F = √d), with the product's own rounding.
        pn = 2 * eps_f * math.sqrt(d) + 4 * d * tol.U_DOUBLE

        def judge(alpha):
            t = tol.hoeffding(shots, 2.0, alpha)
            for i, o in enumerate(e["observables"]):
                raw = rho.expectation(o)
                err = abs(Fraction(rows[i]["raw"]) - raw)
                rep.add(e["id"], "wai.mitigate", "vd: <%s> raw, from %d shots, against the integer density matrix's" % (o, shots), float(err), t + eps_f * math.sqrt(d), err <= Fraction(t + eps_f * math.sqrt(d)))
                num = exact.trace_rho2_pauli(rho, o)
                want = num / purity
                tn, td = t + pn, t + pb
                room = (tn + abs(float(want)) * td) / (float(purity) - td)
                err = abs(Fraction(rows[i]["mitigated"]) - want)
                rep.add(e["id"], "wai.mitigate", "vd: <%s> as Tr(ρ²O)/Tr(ρ²), each from %d shots, against the exact ratio" % (o, shots), float(err), room, err <= Fraction(room), note="exact ratio %.6g" % float(want))

        est.add(3 * len(e["observables"]), judge)
    elif method == "readout":
        p = exact.noise_of(e["noise"])["readout"]
        nb = e["bits"]
        rho = exact.noisy_density(n, gates, e["noise"])
        before = exact.noisy_distribution(rho, measured_of(e), nb)
        listed = {x["bits"]: x["mitigated"] for x in res["distribution"]}
        total = res["details"]["shots"]
        s = 1 / (1 - 2 * float(p))
        spread = s**nb * (1 - float(p)) ** (nb - 1)
        # The engine's inversion: nb passes, each two products and a sum, of values within
        # spread; entries it leaves out are those of magnitude at most 10^-12.
        rounding = 8 * nb * tol.U_DOUBLE * s**nb + 1e-12

        def judge(alpha):
            t = tol.hoeffding(total, spread, alpha)
            worst, bad = 0.0, []
            for x in range(1 << nb):
                k = format(x, "0%db" % nb)
                q = listed.get(k, 0.0)
                err = abs(Fraction(q) - before.get(k, 0))
                worst = max(worst, float(err) / (t + rounding))
                if err > Fraction(t + rounding):
                    bad.append("%s: %.5f, exactly %.5f" % (k, q, float(before.get(k, 0))))
            rep.add(e["id"], "wai.mitigate", "readout: every quasi-probability against the exact law before the flips, largest deviation in units of Hoeffding's bound on %d shots" % total, worst, 1.0, not bad, note="; ".join(bad[:4]) or None)

        est.add(1 << nb, judge)


def check_pulse_entry(e, runner, rep):
    """A pulse program on one qubit against its closed form: P(1) and P(0) within the bound the
    run states (its series' truncation, proved, and its rounding, estimated)."""
    rabi = Fraction(50 * 10**6)

    def frac(x):
        num, _, den = x.partition("/")
        return Fraction(num) / (Fraction(den) if den else 1)

    if len(e["closed"]) == 1:
        c = e["closed"][0]
        p1 = exact.rabi_p1(rabi, frac(c["amplitude"]), frac(c["seconds"]), frac(c.get("detuning_hz", "0")))
    else:
        # A turn about X, then one about an orthogonal axis: P(1) = s2²c1² + c2²s1².
        (c1, s1), (c2, s2) = [exact.cos_sin_turns(rabi * frac(c["amplitude"]) * frac(c["seconds"]) / 2) for c in e["closed"]]
        p1 = Fraction(s2 * s2 * c1 * c1 + c2 * c2 * s1 * s1, 1 << (4 * exact.PREC))
    req = {"qasm": e["qasm"], "method": "pulse", "shots": 0, "pulse": {"dt_ns": 0.5, "levels": 2, "tolerance": 1e-9}}
    ans = runner.run("wai.pulse", {"id": e["id"], "kind": "pulse", "n": 1, "steps": [], "request": req})
    if not ans.get("ok"):
        no_answer(rep, e["id"], "wai.pulse", ans, "a pulse run", fail_on_error=True)
        return
    if not answered(rep, e["id"], "wai.pulse"):
        return
    r = ans["result"]
    bound = r["accuracy"]["bound"]
    pops = r["pulse"]["populations"]
    for k, want in (("1", p1), ("0", 1 - p1)):
        err = abs(Fraction(pops[k]) - want)
        c = rep.add(e["id"], "wai.pulse", "P(%s) against the closed form, within the bound the run states" % k, float(err), bound, err <= Fraction(bound), note="closed form %.12f; %s" % (float(want), r["pulse"]["integrator"].get("steps")))
        c["seconds"] = ans.get("seconds")


def check_decoder_entry(e, runner, rep):
    """Syndromes of a circuit's model decoded, against the least weight of every correction
    (decoding.Graph): exact matching's prediction must be a least-weight correction's up to its
    weights' rounding (tolerance.matching_rounding); union-find must correct every single edge,
    as it states. Every single fault, pairs of faults, and seeded sets of two to six."""
    first = runner.run("wai.matching", {"id": e["id"], "kind": "decode", "circuit": e["qasm"], "syndromes": []})
    if not first.get("ok"):
        no_answer(rep, e["id"], "wai.matching", first, "the decoders' model", fail_on_error=True)
        return
    g = decoding.Graph(first["model"])
    edges = [k for k, _ in g.edges]
    r = corpus.Stream("decode/" + e["id"])
    sets = [[x] for x in edges]
    if e["pairs"] == "all":
        sets += [[a, b] for i, a in enumerate(edges) for b in edges[i + 1:]]
    else:
        sets += [[edges[r.below(len(edges))], edges[r.below(len(edges))]] for _ in range(e["pairs"])]
    sets += [[edges[r.below(len(edges))] for _ in range(2 + r.below(5))] for _ in range(e["sampled"])]
    syndromes = []
    seen = set()
    for st in sets:
        fired, _ = g.syndrome_of(st)
        if fired and tuple(fired) not in seen:
            seen.add(tuple(fired))
            syndromes.append(fired)
    least = [g.least(f) for f in syndromes]
    singles = {}
    for (u, v), mask in g.single.items():
        fired, _ = g.syndrome_of([(u, v, mask)])
        if fired:
            singles[tuple(fired)] = mask
    for dec in e["decoders"]:
        eng = "wai." + dec
        ans = runner.run(eng, {"id": e["id"], "kind": "decode", "circuit": e["qasm"], "syndromes": syndromes})
        if not ans.get("ok"):
            no_answer(rep, e["id"], eng, ans, "decoded syndromes", fail_on_error=True)
            continue
        if not answered(rep, e["id"], eng):
            continue
        preds = [sum(1 << k for k in p) for p in ans["predicted"]]
        if len(preds) != len(syndromes):
            rep.add(e["id"], eng, "predictions returned", len(preds), len(syndromes), False)
            continue
        worst, off, agree = 0.0, [], 0
        for f, w, pm in zip(syndromes, least, preds):
            room = tol.matching_rounding(len(f), g.detectors)
            excess = w.get(pm, math.inf) - min(w.values())
            if excess <= room:
                agree += 1
            ratio = excess / room if room else math.inf
            worst = max(worst, ratio)
            if excess > room and len(off) < 4:
                off.append("%s: predicted %d, least weights %s" % (f, pm, {k: round(v, 4) for k, v in w.items()}))
        if dec == "matching":
            c = rep.add(e["id"], eng, "predictions not a least-weight correction's, largest excess weight in units of its rounding (%d syndromes)" % len(syndromes), worst, 1.0, worst <= 1.0, note="; ".join(off) or None)
        else:
            c = rep.add(e["id"], eng, "predictions that are a least-weight correction's, of %d syndromes (no claim of least weight)" % len(syndromes), agree, None, True, status="info")
            wrong = [f for f, pm in zip(syndromes, preds) if tuple(f) in singles and pm != singles[tuple(f)]]
            rep.add(e["id"], eng, "single edges of the graph not corrected exactly, of %d" % len(singles), len(wrong), 0, not wrong, note=("first: %s" % wrong[:3]) if wrong else None)
        c["seconds"] = ans.get("seconds")


# ------------------------------------------------------------------ the GPU members (the gpu tier)


def gpu_refused(rep, entry, eng, ans, quantity):
    """In the gpu tier a refusal names a device or member this host lacks: recorded as not run,
    with the engine's words. A crash still fails the run."""
    if ans.get("refused"):
        rep.not_run(entry, eng, quantity, str(ans["refused"])[:300])
    else:
        no_answer(rep, entry, eng, ans, quantity, fail_on_error=True)


def check_fabric_entry(e, runner, rep, seen_absent):
    """Each GPU member's bit-exact state, against the CPU's: the same integers and hash as the
    bit-exact method's on the CPU, which the state checks hold to the integer reference within the
    extension's bound. A member this host lacks, or one that refuses the circuit, is recorded as
    not run, with why."""
    import numpy as np

    n = e["n"]
    cpu = runner.run("wai.exact", {"id": e["id"], "kind": "state", "n": n, "steps": e["steps"]})
    if not cpu.get("ok"):
        no_answer(rep, e["id"], "wai.exact", cpu, "the CPU's bit-exact state")
        return
    mine = np.load(cpu["settings"]["fixed"]).reshape(-1)
    ref = exact.run_unitary(n, e["steps"])
    re = exact.ints_from_fixed([int(x) for x in mine[0::2]])
    im = exact.ints_from_fixed([int(x) for x in mine[1::2]])
    d = exact.distance(ref.re, ref.im, re, im)[1]
    bound = tol.fixed_point_bound(n, cpu["ops"])
    rep.add(e["id"], "wai.exact", "the CPU's bit-exact state against the integer reference, modulo a global phase", d, bound, d <= bound)
    ans = runner.run("wai.fabric", {"id": e["id"], "kind": "fabric", "n": n, "steps": e["steps"], "state_out": os.path.join(runner.work, "fabric-%s" % e["id"])})
    if not ans.get("ok"):
        no_answer(rep, e["id"], "wai.fabric", ans, "the fabric's members", fail_on_error=True)
        return
    for a in ans["absent"]:
        key = (a["family"], a.get("key"), a["why"])
        if key not in seen_absent:
            seen_absent.add(key)
            rep.not_run("fabric", "wai.fabric[%s%s]" % (a["family"], (":" + a["key"]) if a.get("key") else ""), "the member's bit-exact states", a["why"][:300])
    for m in ans["members"]:
        name = "wai.fabric[%s:%s]" % (m["family"], m["key"])
        if m.get("refused"):
            rep.not_run(e["id"], name, "the bit-exact state", m["refused"][:300])
            continue
        theirs = np.fromfile(m["state_out"], dtype="<i8")
        same = theirs.shape == mine.shape and bool((theirs == mine).all())
        c = rep.add(e["id"], name, "the bit-exact state's integers and hash, against the CPU's bit-exact method", same and m["state_hash"] == cpu["state_hash"], True, same and m["state_hash"] == cpu["state_hash"], note=m.get("implementation"))
        c["seconds"] = m.get("seconds")


# A single-precision method this build or host does not have, with its words: recorded once.
NO_GPU_METHOD = {}


def check_single_entry(e, runner, rep):
    """The single-precision methods on the platform's GPU against the integer state: the state
    vector (gate by gate at fusion width 3, u = 2^-24) and the cut-rank table, over every
    amplitude or a seeded choice, modulo a global phase in the 2-norm. A method the build or the
    host lacks (`needs a GPU`) is recorded as not run once, and not asked again."""
    n = e["n"]
    engines = [x for x in ["wai.statevector_single"] + (["wai.cutrank_single"] if e["family"] == "cutrank" else []) if x not in NO_GPU_METHOD]
    if not engines:
        return
    ref = exact.run_unitary(n, gates_of(e))
    bits = chosen_bits(e, ref, None if n <= 12 else 64)
    for eng in engines:
        ans = run_amplitudes(runner, eng, {"id": e["id"], "n": n, "steps": gates_of(e)}, bits)
        if not ans.get("ok") and "needs a GPU" in str(ans.get("refused")):
            NO_GPU_METHOD[eng] = ans["refused"]
            rep.not_run("gpu", eng, "amplitudes on the GPU in single precision", str(ans["refused"])[:300])
            continue
        if not ans.get("ok"):
            gpu_refused(rep, e["id"], eng, ans, "amplitudes on the GPU in single precision")
            continue
        dm, d = amplitudes_distance(ref, ans["amplitudes"], bits)
        if eng == "wai.statevector_single":
            bound = tol.single_state(tol.gate_count(e["steps"]), n, width=3)
        else:
            bound = math.sqrt(len(bits)) * tol.cutrank(n) * tol.U_SINGLE / tol.U_DOUBLE
        c = rep.add(e["id"], eng, "%d amplitudes in single precision against the integer state, 2-norm modulo a global phase" % len(bits), dm, bound, dm <= bound, note=str((ans.get("accuracy") or {}).get("device") or ""))
        c["strict"] = d
        c["seconds"] = ans.get("seconds")


GPU_FAMILIES = ("bit_order", "ghz", "qft", "grover", "qpe", "random_clifford_t", "random_dyadic", "mirror", "gate_conventions", "clifford_large", "random_rotations", "cutrank")


def gpu_entries(es):
    """The gpu tier's circuits: those the bit-exact method runs (by the reader's rule) up to 20
    qubits, for the fabric's members; and the double-precision state entries and cut-rank sums up
    to 20 qubits, for the single-precision methods."""
    fabric, single = [], []
    for e in es:
        if e["family"] not in GPU_FAMILIES or e["n"] > 20 or any("g" not in s for s in e["steps"]):
            continue
        if e["expect"] in ("exact_state", "known_state", "mirror") and corpus.reader_class(e)[0] == "exact" and "wai.exact" in e.get("engines", ["wai.exact"]):
            fabric.append(e)
        if e["expect"] in ("exact_state", "known_state", "cutrank") and len(e["steps"]) < 4096:
            single.append(e)
    return fabric, single


# ------------------------------------------------------------------ the run


def run_all(runner, jobs):
    with ThreadPoolExecutor(max_workers=runner.jobs) as ex:
        futs = {eng: ex.submit(runner.run, eng, job) for eng, job in jobs.items()}
        return {eng: f.result() for eng, f in futs.items()}


def host():
    info = {"system": platform.system(), "machine": platform.machine(), "release": platform.release()}
    try:
        if platform.system() == "Darwin":
            info["cpu"] = subprocess.run(["sysctl", "-n", "machdep.cpu.brand_string"], capture_output=True, text=True).stdout.strip()
            info["memory_bytes"] = int(subprocess.run(["sysctl", "-n", "hw.memsize"], capture_output=True, text=True).stdout)
            info["os"] = "macOS " + platform.mac_ver()[0]
        else:
            for line in open("/proc/cpuinfo"):
                if line.startswith("model name"):
                    info["cpu"] = line.split(":", 1)[1].strip()
                    break
            for line in open("/proc/meminfo"):
                if line.startswith("MemTotal"):
                    info["memory_bytes"] = int(line.split()[1]) * 1024
            try:
                for line in open("/etc/os-release"):
                    if line.startswith("PRETTY_NAME="):
                        info["os"] = line.split("=", 1)[1].strip().strip('"')
            except OSError:
                pass
        info["logical_cpus"] = os.cpu_count()
        info["load_average"] = list(os.getloadavg())
    except Exception:  # noqa: BLE001
        pass
    return info


def run_reproducer(path):
    """A committed reproducer of another engine's disagreement: it prints one JSON object per
    line, each a claim, what it saw, and whether the claim holds."""
    out = subprocess.run([sys.executable, "-I", os.path.join(ROOT, path)], capture_output=True, text=True, timeout=600)
    claims = []
    for line in out.stdout.splitlines():
        line = line.strip()
        if line.startswith("{"):
            claims.append(json.loads(line))
    if out.returncode != 0 and not claims:
        claims.append({"claim": "the reproducer runs", "observed": out.stderr[-500:], "holds": False})
    return claims


def apply_known(rep, entries_run):
    """A disagreement of another engine with the exact reference, found and examined before, is
    listed in corpus/known-disagreements.json with what was found and a reproducer committed beside
    it. A failing check that matches one is reported as `known` rather than as a failure of this
    run, and the reproducer's claims are rerun and recorded; one listed that no longer shows is
    reported too. Only other engines may be listed: a failure of the engine under test is always a
    failure."""
    known = json.load(open(os.path.join(ROOT, "corpus", "known-disagreements.json")))["disagreements"]
    seen = set()
    for i, k in enumerate(known):
        if k["engine"].startswith("wai."):
            rep.add(k["entry"], k["engine"], "a known disagreement listed for the engine under test", True, False, False, note="only other engines may be listed")
            seen.add(i)
        if not os.path.exists(os.path.join(ROOT, k.get("reproducer", ""))) or not k.get("reproducer"):
            rep.add(k["entry"], k["engine"], "a known disagreement has a reproducer committed beside it", k.get("reproducer"), True, False)
            seen.add(i)
    for c in rep.checks:
        if c["status"] != "fail":
            continue
        for i, k in enumerate(known):
            quantities = k["quantity"] if isinstance(k["quantity"], list) else [k["quantity"]]
            if i not in seen and c["entry"] == k["entry"] and c["engine"] == k["engine"] and any(x in c["quantity"] for x in quantities):
                c["status"] = "known"
                c["note"] = ((c.get("note") or "") + " | known: " + k["finding"]).strip(" |")
                c["right"] = k["right"]
                k["_matched"] = True
    done = set()
    for k in known:
        if k["entry"] not in entries_run or k["engine"].startswith("wai."):
            continue
        if not k.get("_matched"):
            rep.add(k["entry"], k["engine"], "a listed disagreement no longer shows", False, True, True, status="info", note=k["finding"])
        if k.get("reproducer") and k["reproducer"] not in done and os.path.exists(os.path.join(ROOT, k["reproducer"])):
            done.add(k["reproducer"])
            for cl in run_reproducer(k["reproducer"]):
                rep.add(k["entry"], k["engine"], "reproducer (%s): %s" % (k["reproducer"], cl["claim"]), cl["observed"], cl.get("expected"), True, status="known" if cl["holds"] else "info", note=None if cl["holds"] else "the claim no longer holds")


def check_expected(rep, entries_run):
    """The expected-refusals file is what corpus.py's rule and list give, and a listed refusal that
    no longer comes is reported (a `may` that the engine now answers)."""
    want = corpus.expected_refusals()
    rep.add("corpus", "corpus.py", "corpus/expected-refusals.json is the reader's rule and corpus.py's list", EXPECTED.rows == want, True, EXPECTED.rows == want)
    for entry, eng in sorted(EXPECTED.answered):
        x = EXPECTED.get(entry, eng)
        if x["must"] == "may" and entry in entries_run:
            rep.add(entry, eng, "a listed refusal no longer comes: the engine answered", True, False, True, status="info", note=x["reason"])


def phase_table_check(rep):
    """quantum-sim.md section 3.2's table: P(k) for k = 1..32, cos and sin of 2π/2^k rounded to the
    nearest 2^-30. The integer reference must give the referee's table, entry for entry."""
    table = referee().derive_v2()[0]
    half = 1 << (exact.PREC - 31)
    wrong = []
    for k in range(1, 33):
        c, s = exact.cos_sin("t:1/%d" % k)
        mine = ((c + half) >> (exact.PREC - 30), (s + half) >> (exact.PREC - 30))
        if mine != tuple(table[k]):
            wrong.append(k)
    rep.add("corpus", "exact", "the extension's phase table P(1..32): the integer reference's cos and sin, rounded to the nearest 2^-30, against the referee's", wrong, [], not wrong)


def manifest_check(es, rep):
    """The generator still gives the bytes the manifest pins."""
    m = json.load(open(os.path.join(ROOT, "corpus", "manifest.json")))
    pinned = {r["id"]: r.get("sha256") for r in m["circuits"]}
    drift = [e["id"] for e in es if pinned.get(e["id"]) != hashlib.sha256(e["qasm"].encode()).hexdigest()]
    for q in corpus.qec_entries():
        text = open(os.path.join(ROOT, "corpus", "qec", q["id"] + ".circuit")).read()
        if pinned.get(q["id"]) != hashlib.sha256(text.encode()).hexdigest():
            drift.append(q["id"])
    rep.add("corpus", "corpus.py", "circuits whose text differs from the manifest", drift, [], not drift)
    for e in es:
        p = os.path.join(ROOT, "corpus", "qasm", e["id"] + ".qasm")
        if os.path.exists(p) and open(p).read() != e["qasm"]:
            rep.add(e["id"], "corpus.py", "the committed program equals the generator's", False, True, False)


CI_HAND = ("qec_inverted_and_tagged", "qec_correlated", "qec_refused_random_detector")


def main():
    global EXPECTED
    ap = argparse.ArgumentParser()
    ap.add_argument("--tier", default="ci", choices=["ci", "full", "gpu"])
    ap.add_argument("--out")
    ap.add_argument("--jobs", type=int, default=2)
    ap.add_argument("--timeout", type=float, default=900)
    ap.add_argument("--memory-gb", type=float, help="stop a job whose processes hold more than this")
    ap.add_argument("--inject", choices=sorted(FAULTS))
    ap.add_argument("--only")
    ap.add_argument("--keep", action="store_true")
    a = ap.parse_args()
    if not os.environ.get("WAI_CROSS_SIMULATOR"):
        print("compare.py: WAI_CROSS_SIMULATOR must name the engine's driver (reproduce.sh builds it)", file=sys.stderr)
        sys.exit(2)
    work = tempfile.mkdtemp(prefix="qcorpus-")
    runner = Runner(work, sys.executable, a.timeout, a.inject, a.memory_gb)
    runner.jobs = a.jobs
    rep = Report()
    EXPECTED = Expected()
    started = time.time()
    es = corpus.build()
    manifest_check(es, rep)
    phase_table_check(rep)
    only = set(a.only.split(",")) if a.only else None
    if a.tier == "gpu":
        fabric, single = gpu_entries(es)
        if only:
            fabric = [e for e in fabric if e["id"] in only]
            single = [e for e in single if e["id"] in only]
        seen_absent = set()
        for e in fabric:
            print("· fabric", e["id"], file=sys.stderr, flush=True)
            check_fabric_entry(e, runner, rep, seen_absent)
        for e in single:
            print("· single", e["id"], file=sys.stderr, flush=True)
            check_single_entry(e, runner, rep)
        # What the host declares it has (QCORPUS_GPU_EXPECT: family:key for a fabric member, or a
        # method's engine name) must have passed a check, or the run fails: a host whose GPU went
        # missing is a failure, not a run with every member not run.
        for want in filter(None, os.environ.get("QCORPUS_GPU_EXPECT", "").split(",")):
            eng = "wai.fabric[%s]" % want if ":" in want else want
            ran = sum(1 for c in rep.checks if c["engine"] == eng and c["status"] == "pass")
            rep.add("fabric", eng, "checks passed by what this host declares it has", ran, ">0", ran > 0)
        return finish(a, rep, work, started)
    if a.tier == "ci":
        es = [e for e in es if e["tier"] == "ci"]
    if only:
        es = [e for e in es if e["id"] in only]
    count_plan = []
    est = Estimates()
    alpha = FAMILY_WISE_ALPHA / SHARES
    stab_counts = sum(1 for e in es if e["expect"] == "stabilizer_counts")
    for e in es:
        print("·", e["id"], file=sys.stderr, flush=True)
        x = e["expect"]
        if x in ("exact_state", "known_state"):
            check_state_entry(e, runner, rep, e.get("engines") or STATE_ENGINES)
        elif x == "distribution":
            check_counts_entry(e, runner, rep, count_plan)
        elif x == "stabilizer":
            check_stabilizer_entry(e, runner, rep)
        elif x == "stabilizer_counts":
            check_stabilizer_counts(e, runner, rep, alpha / max(stab_counts, 1))
        elif x == "mirror":
            check_mirror_entry(e, runner, rep)
        elif x == "pairwise":
            check_pairwise_entry(e, runner, rep)
        elif x in ("amplitudes", "wide_pairwise", "wide_mirror"):
            check_wide_entry(e, runner, rep)
        elif x == "expectation":
            check_expectation_entry(e, runner, rep)
        elif x == "pauli":
            check_pauli_entry(e, runner, rep)
        elif x == "stabrank":
            check_stabrank_entry(e, runner, rep, est)
        elif x == "cutrank":
            check_cutrank_entry(e, runner, rep)
        elif x == "noisy":
            check_noisy_entry(e, runner, rep, count_plan, est)
        elif x == "mitigation":
            check_mitigation_entry(e, runner, rep, est)
        elif x == "pulse":
            check_pulse_entry(e, runner, rep)
        elif x == "decoders":
            check_decoder_entry(e, runner, rep)
    judge_counts(count_plan, rep, alpha)
    est.judge(alpha)
    qplan = []
    qs = [q for q in corpus.qec_entries() if a.tier == "full" or q["tier"] == "ci"]
    hand = json.load(open(os.path.join(ROOT, "corpus", "qec-hand-written.json")))
    if only:
        qs = [q for q in qs if q["id"] in only]
        hand = {k: v for k, v in hand.items() if k in only}
    for q in qs:
        print("·", q["id"], file=sys.stderr, flush=True)
        text = qec_text(q, runner, rep)
        check_qec_entry(q, text, runner, rep, qplan)
        if q["code"] == "surface_code:rotated_memory_z" and os.environ.get("WAI_QUANTUM"):
            mine = subprocess.run([os.environ["WAI_QUANTUM"], "lab", "qec", "--circuit", "-d", str(q["distance"]), "--rounds", str(q["rounds"]), "--p", repr(q["p"])], capture_output=True, text=True).stdout
            rep.add(q["id"], "wai.qec", "the engine's own generator writes the same text, byte for byte", mine == text, True, mine == text)
    for k, v in hand.items():
        if a.tier == "ci" and k not in CI_HAND:
            continue
        print("·", k, file=sys.stderr, flush=True)
        check_qec_entry({"id": k, "tier": a.tier}, v["text"], runner, rep, qplan)
    judge_qec(qplan, rep, alpha)
    entries_run = {e["id"] for e in es} | {q["id"] for q in qs} | {k for k in hand if a.tier == "full" or k in CI_HAND}
    apply_known(rep, entries_run)
    check_expected(rep, entries_run)
    return finish(a, rep, work, started)


def finish(a, rep, work, started):
    """The results file, every failure and refusal printed, and the verdict."""
    vers = json.loads(subprocess.run([sys.executable, "-I", os.path.join(HERE, "engines.py"), "--versions"], capture_output=True, text=True).stdout)
    statuses = {}
    for c in rep.checks:
        statuses[c["status"]] = statuses.get(c["status"], 0) + 1
    out = {
        "about": "One run of the cross-simulator corpus: every check, its measured value and its tolerance. tools/compare.py writes it.",
        "tier": a.tier,
        "inject": a.inject,
        "host": host(),
        "versions": vers,
        "engine_under_test": {"driver": os.path.basename(os.environ["WAI_CROSS_SIMULATOR"]), "build": os.environ.get("WAI_BUILD_INFO", "")},
        "seconds": round(time.time() - started, 1),
        "summary": statuses,
        "checks": rep.checks,
    }
    if a.out:
        with open(a.out, "w") as f:
            json.dump(out, f, indent=1, sort_keys=False, default=str)
            f.write("\n")
    fails = [c for c in rep.checks if c["status"] == "fail"]
    for c in rep.checks:
        if c["status"] in ("fail", "refused", "not run", "known"):
            print("%-8s %-34s %-34s %s | measured %s, tolerance %s%s" % (c["status"].upper(), c["entry"][:34], c["engine"][:34], c["quantity"][:110], fmt(c["measured"]), fmt(c["tolerance"]), (" — " + c["note"][:240]) if c.get("note") else ""))
    print("summary: %s; %.0f s" % (", ".join("%s %d" % kv for kv in sorted(statuses.items())), time.time() - started))
    print("FAIL" if fails else "PASS")
    if not a.keep:
        shutil.rmtree(work, ignore_errors=True)
    sys.exit(1 if fails else 0)


def fmt(x):
    if isinstance(x, float):
        return "%.3g" % x
    return str(x)[:80]


if __name__ == "__main__":
    main()