#!/usr/bin/env python3
"""Recompute every number on this site from the raw data files beside it.

Run from the site folder:  python3 -I reproduce/reproduce.py          (text report)
                           python3 -I reproduce/reproduce.py --json   (machine-readable)
Standard library only. Deterministic: no random numbers, no clock, no network. Two runs give identical output.
Inputs: data/experiments-2/*.csv, data/experiments-3/*.csv, data/experiments-3/predictions.json.
Constants: CODATA 2018 and the IAU solar mass parameter, written out below.
"""
import csv, json, math, os, sys

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(os.path.dirname(HERE), "data")

# ---------- constants (CODATA 2018; IAU) ----------
C = 299792458.0
G = 6.67430e-11
GM_SUN = 1.32712440018e20            # m^3 s^-2
A_MERCURY, E_MERCURY, P_MERCURY_DAYS = 5.7909050e10, 0.205630, 87.9691
RINF_C = 3.2898419602508e15          # Hz
ALPHA = 7.2973525693e-3
MP_ME = 1836.15267343
MEC2_H = 0.51099895000e6 / 4.135667696e-15   # Hz
TAU0_MU, GAMMA_BAILEY = 2.1969811, 29.327    # us; Bailey et al. 1977
BAILEY_TOTAL_REL = 2e-3 / 1.96               # their 2e-3 at 95 % read as 1 sigma


def rows(path):
    with open(path, newline="", encoding="utf-8") as f:
        return [r for r in csv.DictReader(f)]


def e2(name):
    return rows(os.path.join(DATA, "experiments-2", name + ".csv"))


def pooled(r):
    return [x for x in r if x.get("pooled", "yes") == "yes"]


def chi2_sf(x, k):
    """Upper tail of chi-square with k degrees of freedom (series / continued fraction, Numerical Recipes 6.2)."""
    a, x = k / 2.0, x / 2.0
    if x <= 0:
        return 1.0
    if x > 700:
        return 0.0
    gl = math.lgamma(a)
    if x < a + 1:
        s = d = 1.0 / a; ap = a
        for _ in range(1000):
            ap += 1; d *= x / ap; s += d
            if abs(d) < abs(s) * 1e-15:
                break
        return max(0.0, 1.0 - s * math.exp(-x + a * math.log(x) - gl))
    b = x + 1 - a; c = 1e300; d = 1 / b; h = d
    for i in range(1, 1000):
        an = -i * (i - a); b += 2; d = an * d + b; d = 1 / d if d else 1e300
        c = b + an / c; c = c if c else 1e-300; de = d * c; h *= de
        if abs(de - 1) < 1e-15:
            break
    return math.exp(-x + a * math.log(x) - gl) * h


def credibility(y, s):
    """Inverse-variance mean with a between-experiment variance tau^2 by the method of moments; Cochran's Q."""
    w = [1 / si ** 2 for si in s]
    m0 = sum(wi * yi for wi, yi in zip(w, y)) / sum(w)
    Q = sum(wi * (yi - m0) ** 2 for wi, yi in zip(w, y))
    n = len(y)
    tau2, pQ = 0.0, float("nan")
    if n > 1:
        denom = sum(w) - sum(wi ** 2 for wi in w) / sum(w)
        tau2 = max(0.0, (Q - (n - 1)) / denom)
        pQ = chi2_sf(Q, n - 1)
    ws = [1 / (si ** 2 + tau2) for si in s]
    mu = sum(wi * yi for wi, yi in zip(ws, y)) / sum(ws)
    return {"n": n, "mean": mu, "se": math.sqrt(1 / sum(ws)), "Q": Q, "pQ": pQ, "tau2": tau2,
            "one_quantity": (n < 2) or pQ > 0.05}


def test(y, s, pred, label, sp=None, common=0.0):
    """Standardized residuals of each measurement against a prediction; chi-square (with a common error if given)."""
    sp = sp or [0.0] * len(y)
    r = [(yi - pi) / math.sqrt(si ** 2 + spi ** 2 + common ** 2) for yi, si, pi, spi in
         zip(y, s, pred if isinstance(pred, list) else [pred] * len(y), sp)]
    if common and len(y) == 2:   # 2x2 covariance with the common error
        d = [yi - (pred if not isinstance(pred, list) else pred[i]) for i, yi in enumerate(y)]
        a, b, c = s[0] ** 2 + common ** 2, common ** 2, s[1] ** 2 + common ** 2
        det = a * c - b * b
        chi2 = (c * d[0] ** 2 - 2 * b * d[0] * d[1] + a * d[1] ** 2) / det
    else:
        chi2 = sum(x * x for x in r)
    p = chi2_sf(chi2, len(y))
    agree_all = all(abs(x) <= 2 for x in r)
    return {"label": label, "residuals": [round(x, 2) for x in r], "chi2": chi2, "dof": len(y), "p": p,
            "result": "agrees" if agree_all else "contradicted"}


def mercury_gr():
    per_orbit = 6 * math.pi * GM_SUN / (C ** 2 * A_MERCURY * (1 - E_MERCURY ** 2))
    return per_orbit * (36525.0 / P_MERCURY_DAYS) * 180 / math.pi * 3600


def dirac_level(n, j, mu_c2_h):
    k = j + 0.5
    delta = k - math.sqrt(k * k - ALPHA ** 2)
    return mu_c2_h * (1 / math.sqrt(1 + (ALPHA / (n - delta)) ** 2) - 1)


def experiments2():
    out = {}
    num = lambda r, k: float(r[k])
    for key, name, alts in [
        ("E1", "E1-redshift", [(0.0, "alpha = 0 (clock rate follows the potential)")]),
        ("E2", "E2-light", [(0.0, "gamma - 1 = 0 (full bending, general relativity)"),
                            (-1.0, "gamma - 1 = -1 (clock-rate half only, Soldner's 0.876 arcsec)")]),
        ("E6", "E6-sheet", [(720.0, "return at 720 degrees (two-valued wave)"),
                            (360.0, "return at 360 degrees (ordinary vector, rotating body)")]),
        ("E9", "E9-antimatter", [(1.0, "antimatter falls at 1 g"), (-1.0, "antimatter rises at -1 g")]),
    ]:
        r = pooled(e2(name)); y = [num(x, "value") for x in r]; s = [num(x, "sigma") for x in r]
        out[key] = {"credibility": credibility(y, s), "tests": [test(y, s, p, l) for p, l in alts]}
    # E3: the authors' own per-flight predictions
    r = e2("E3-aircraft"); y = [num(x, "value") for x in r]; s = [num(x, "sigma") for x in r]
    pr = [num(x, "prediction") for x in r]; ps = [num(x, "prediction_sigma") for x in r]
    out["E3"] = {"credibility": credibility([a - b for a, b in zip(y, pr)], [math.hypot(a, b) for a, b in zip(s, ps)]),
                 "tests": [test(y, s, pr, "Hafele and Keating's own prediction from the flight logs", sp=ps)]}
    # E4: one gamma for both signs, so a common prediction error
    r = e2("E4-muon"); y = [num(x, "value") for x in r]; s = [num(x, "sigma") for x in r]
    pred = GAMMA_BAILEY * TAU0_MU
    cr = credibility(y, s)
    common = math.sqrt(max(0.0, (BAILEY_TOTAL_REL * pred) ** 2 - cr["se"] ** 2))
    out["E4"] = {"credibility": cr, "prediction": pred, "common_sigma": common,
                 "tests": [test(y, s, pred, "gamma x tau0 (time dilation)", common=common)]}
    # E5: CHSH
    r = pooled(e2("E5-chsh")); y = [num(x, "value") for x in r]; s = [num(x, "sigma") for x in r]
    t = [test(y, s, 2 * math.sqrt(2), "S = 2 sqrt 2 at every apparatus (ideal analysers)"),
         test(y, s, 2.0, "S at most 2 (each end holds its own preset value)"),
         test(y[:1], s[:1], [2.70], "Aspect 1982 against its own apparatus prediction 2.70 +/- 0.05", sp=[0.05])]
    inside = all(2 - 2 * si <= yi <= 2 * math.sqrt(2) + 2 * si for yi, si in zip(y, s))
    out["E5"] = {"credibility": credibility(y, s), "tests": t,
                 "bound": {"label": "2 < S <= 2 sqrt 2", "result": "agrees" if inside and all(yi > 2 for yi in y) else "contradicted"},
                 "sigma_above_2": [round((yi - 2) / si, 1) for yi, si in zip(y, s)]}
    # E7: Mercury
    r = pooled(e2("E7-perihelion")); y = [num(x, "value") for x in r]; s = [num(x, "sigma") for x in r]
    gr = mercury_gr()
    out["E7"] = {"prediction": gr, "credibility": credibility(y, s),
                 "tests": [test(y, s, gr, "general relativity, 6 pi GM / (c^2 a (1 - e^2)) per orbit"),
                           test(y, s, 0.0, "weak-field clock-rate map alone (no advance)")],
                 "park_beta_minus_1": {"value": -2.7e-5, "sigma": 3.9e-5, "z": round(-2.7 / 3.9, 2)}}
    # E8: hydrogen 1S-2S
    r = pooled(e2("E8-hydrogen")); y = [num(x, "value") for x in r]; s = [num(x, "sigma") for x in r]
    red = 1 / (1 + 1 / MP_ME)
    bohr = RINF_C * 0.75 * red
    dirac = dirac_level(2, 0.5, MEC2_H * red) - dirac_level(1, 0.5, MEC2_H * red)
    out["E8"] = {"bohr": bohr, "dirac": dirac, "miss_bohr_hz": y[0] - bohr, "miss_dirac_hz": y[0] - dirac,
                 "tests": [test(y, s, bohr, "Bohr ladder with reduced mass"),
                           test(y, s, dirac, "Dirac levels with reduced mass (no QED)")]}
    return out


def experiments3():
    P = json.load(open(os.path.join(DATA, "experiments-3", "predictions.json"), encoding="utf-8"))
    preds = {p["id"]: p for p in P["predictions"] if p.get("kind") == "number"}
    meas = {}
    d3 = os.path.join(DATA, "experiments-3")
    for fn in sorted(os.listdir(d3)):
        if fn.endswith(".csv") and fn != "sims.csv":
            for r in rows(os.path.join(d3, fn)):
                meas[r["id"]] = r
    res = []
    for pid, p in preds.items():
        m = meas.get(pid) or meas.get(pid.rstrip("abd"))
        if pid == "S4" and "S4a" in meas:   # the line ratio J 2-1 over J 1-0, errors propagated
            a, b = meas["S4a"], meas["S4b"]
            va, vb, sa, sb = float(a["value"]), float(b["value"]), float(a["sigma"]), float(b["sigma"])
            q = vb / va
            m = dict(b, value=repr(q), sigma=repr(q * math.hypot(sa / va, sb / vb)), bound_low="", recalled="no")
        row = {"id": pid, "quantity": p["quantity"], "field": p["field"], "predicted": p["value"],
               "predicted_sigma": p["sigma"], "prediction": p["how"]}
        if not m:
            row.update(result="open", note="no measured value entered"); res.append(row); continue
        v, sv = float(m["value"]), float(m["sigma"])
        row.update(measured=v, measured_sigma=sv, source=m["source"], year=m["year"])
        bd = m.get("bound_low", "")
        if bd:
            lo, hi = float(m["bound_low"]), float(m["bound_high"])
            row.update(z=None, result="agrees (inside the published bound)" if lo <= p["value"] <= hi else "contradicted (outside the bound)")
        else:
            den = math.sqrt(sv ** 2 + p["sigma"] ** 2)
            if den == 0:
                row.update(z=0.0 if v == p["value"] else None,
                           result="agrees" if v == p["value"] else "contradicted (exact integers differ)")
            else:
                z = (v - p["value"]) / den
                row.update(z=z, result="agrees" if abs(z) <= 2 else "contradicted")
        if m.get("recalled") == "yes":
            row["result"] += " (data recalled from a secondary source, not counted)"
        res.append(row)
    agree = sum(r["result"].startswith("agrees") and "recalled" not in r["result"] for r in res)
    contradicted = sum(r["result"].startswith("contradicted") for r in res)
    return {"rows": res, "scored": sum(r["result"] != "open" for r in res), "agree": agree, "contradicted": contradicted,
            "open": sum(r["result"] == "open" for r in res),
            "recalled": sum("recalled" in r["result"] for r in res)}


def rope():
    """Three readings of an entangled pair as one wave tied at two ends, scored by the CHSH number S."""
    N = 36000

    def prob(a, b):  # each end reads a shared preset phase on its own, probabilistically
        s = 0.0
        for i in range(N):
            l = 2 * math.pi * i / N
            s += (2 * math.cos((a - l) / 2) ** 2 - 1) * (2 * math.cos((b - l) / 2) ** 2 - 1)
        return s / N

    def sign(a, b):  # each end reads the sign of a shared preset phase against its analyser
        s = 0
        for i in range(N):
            l = 2 * math.pi * i / N
            s += (1 if math.cos(a - l) > 0 else -1) * (1 if math.cos(b - l) > 0 else -1)
        return s / N

    def once(a, b):  # the wave read once for both ends: the textbook singlet correlation, taken as given
        return -math.cos(a - b)

    def chsh(E):
        a, a2, b, b2 = 0.0, math.pi / 2, math.pi / 4, 3 * math.pi / 4
        return abs(E(a, b) - E(a, b2) + E(a2, b) + E(a2, b2))

    return {"two_ends_probabilistic": round(chsh(prob), 4), "two_ends_sign": round(chsh(sign), 4),
            "read_once_for_both_ends": round(chsh(once), 4), "tsirelson": round(2 * math.sqrt(2), 4)}


def proof(e2r, e3r):
    """Each claim on the Proof page: the measurement, number, uncertainty, source DOI, and the computed distance."""
    z = lambda v, s, p: (v - p) / s
    e3 = {r["id"]: r for r in e3r["rows"]}
    C_ = []
    add = lambda cid, claim, meas, v, s, doi, comp: C_.append(
        {"id": cid, "claim": claim, "measurement": meas, "value": v, "uncertainty": s, "doi": doi, "computed": comp})
    add("P01", "The two-valued wave returns after 720 degrees of rotation, not 360.",
        "Neutron interferometer, precession angle at which the interference returns (Rauch et al. 1978)", "716.8 deg", "3.8 deg (1 sigma)",
        "10.1007/BF01321193", "z from 720: %+.2f; z from 360: %+.1f" % (z(716.8, 3.8, 720), z(716.8, 3.8, 360)))
    add("P02", "No pair of preset answers carried by the two ends reproduces the Bell correlations.",
        "CHSH number S, photon pairs (Aspect, Grangier, Roger 1982)", "S = 2.697", "0.015 (1 sigma)",
        "10.1103/PhysRevLett.49.91", "%.1f sigma above the preset-answer limit 2" % z(2.697, 0.015, 2))
    add("P03", "The same holds with the detection and locality loopholes closed together.",
        "CHSH number S, electron spins 1.3 km apart (Hensen et al. 2015)", "S = 2.42", "0.20 (1 sigma)",
        "10.1038/nature15759", "%.1f sigma above 2" % z(2.42, 0.20, 2))
    add("P04", "The correlation survives over 1,203 km; it is not tied to places.",
        "CHSH number S, satellite to two ground stations (Yin et al. 2017)", "S = 2.37", "0.09 (1 sigma)",
        "10.1126/science.aan3211", "%.1f sigma above 2" % z(2.37, 0.09, 2))
    add("P05", "No large fourth spatial direction is open to gravity at laboratory range.",
        "Torsion balance, inverse-square law (Lee et al. 2020)", "Yukawa range of gravitational strength < 38.6 um", "95 % confidence bound",
        "10.1103/PhysRevLett.124.101101", "a three-direction inverse square predicts no deviation: inside the bound")
    add("P06", "No small curled direction shows its tower of heavy copies, in the models searched.",
        "High-mass diphoton events at 13 TeV (CMS 2018)", "first Randall-Sundrum graviton > 2.3 to 4.6 TeV", "95 % confidence bound (couplings 0.01 to 0.2)",
        "10.1103/PhysRevD.98.092001", "no excitation seen: a fixed mass inside the range is excluded")
    a = e2r["E1"]["credibility"]
    add("P07", "A clock's rate follows the gravitational potential exactly (redshift coefficient alpha = 0).",
        "Clock comparisons: Pound and Snider 1965, Gravity Probe A 1980, Galileo 5 and 6 (Delva et al. 2018), combined", "alpha = %.3g" % a["mean"], "%.3g (1 sigma)" % a["se"],
        "10.1103/PhysRevLett.121.231101", "combined %+.2f sigma from 0; chi2 = %.4f, p = %.3f" % (a["mean"] / a["se"], e2r["E1"]["tests"][0]["chi2"], e2r["E1"]["tests"][0]["p"]))
    b = e2r["E2"]["credibility"]
    add("P08", "Light past the Sun bends by the full amount (gamma = 1); the clock-rate half alone is contradicted.",
        "Cassini Shapiro delay (Bertotti et al. 2003) with geodetic VLBI (Lambert and Le Poncin-Lafitte 2011)", "gamma - 1 = %.3g" % b["mean"], "%.3g (1 sigma)" % b["se"],
        "10.1038/nature01997", "full bending %+.2f sigma; half bending %+.0f sigma" % (b["mean"] / b["se"], (b["mean"] + 1) / b["se"]))
    add("P09", "A moving clock runs slow by the time-dilation factor, at 10^18 g of acceleration.",
        "Muon lifetimes in the CERN storage ring (Bailey et al. 1977)", "64.368 and 64.419 us", "0.029 and 0.058 us, plus the common gamma error",
        "10.1038/268301a0", "against %.4f us: residuals %s" % (e2r["E4"]["prediction"], e2r["E4"]["tests"][0]["residuals"]))
    add("P10", "Flown clocks gain and lose exactly as gravity and speed predict.",
        "Atomic clocks flown round the world (Hafele and Keating 1972)", "east -59 ns, west +273 ns", "10 ns and 7 ns (1 sigma)",
        "10.1126/science.177.4044.168", "residuals against the flight-log predictions %s" % e2r["E3"]["tests"][0]["residuals"])
    add("P11", "Antimatter falls; it does not rise.",
        "Antihydrogen in a vertical trap (ALPHA-g, Anderson et al. 2023)", "0.75 g", "0.21 g (1 sigma, statistical and simulation in quadrature)",
        "10.1038/s41586-023-06527-1", "falls at 1 g: %+.2f sigma; rises at -1 g: %+.1f sigma" % (z(0.75, 0.206, 1), z(0.75, 0.206, -1)))
    add("P12", "A spinning Earth drags the gyroscopes; a static, rotation-free gravity map is contradicted.",
        "Gravity Probe B frame-dragging drift (Everitt et al. 2011)", "37.2 mas/yr", "7.2 mas/yr (1 sigma)",
        "10.1103/PhysRevLett.106.221101", "no rotation term predicts 0: %+.2f sigma" % z(37.2, 7.2, 0))
    add("P13", "Mercury's perihelion advance matches general relativity (beta = 1).",
        "Ranging to MESSENGER, PPN beta with Cassini gamma (Park et al. 2017)", "beta - 1 = -2.7e-5", "3.9e-5 (1 sigma)",
        "10.3847/1538-3881/aa5be2", "beta = 1: %+.2f sigma" % (-2.7 / 3.9))
    h = e2r["E8"]
    add("P14", "Integer energy levels alone do not give hydrogen's 1S-2S interval; the QED shifts are needed.",
        "Hydrogen 1S-2S two-photon spectroscopy (Parthey et al. 2011)", "2 466 061 413 187 035 Hz", "10 Hz (1 sigma)",
        "10.1103/PhysRevLett.107.203001", "Bohr ladder misses by %.4g Hz; Dirac levels by %.4g Hz" % (h["miss_bohr_hz"], h["miss_dirac_hz"]))
    t = e3["T11"]
    add("P15", "All bodies fall alike (titanium and platinum).",
        "MICROSCOPE satellite, Eotvos parameter (Touboul et al. 2022)", "eta = -1.5e-15", "2.7e-15 (1 sigma)",
        "10.1103/PhysRevLett.129.121102", "universal free fall predicts 0: %+.2f sigma" % t["z"])
    t = e3["T12"]
    add("P16", "The double pulsar loses orbital energy to gravitational waves at the predicted rate.",
        "Double pulsar J0737-3039 orbital decay, observed over predicted (Kramer et al. 2021)", "0.999963", "6.3e-5 (1 sigma)",
        "10.1103/PhysRevX.11.041050", "ratio 1: %+.2f sigma" % t["z"])
    t = e3["T2"]
    add("P17", "Clocks one millimetre apart in height tick at different rates, by g/c^2 per metre.",
        "Strontium lattice clock, rate gradient across the sample (Bothwell et al. 2022)", "9.8e-20 per mm", "2.3e-20 per mm (1 sigma)",
        "10.1038/s41586-021-04349-7", "against %.3g per mm: %+.2f sigma" % (t["predicted"], t["z"]))
    t = e3["T5"]
    add("P18", "Light leaving the Sun is redshifted by the potential difference to Earth.",
        "Solar spectrum via Moon light, HARPS laser comb (Gonzalez Hernandez et al. 2020)", "638 m/s", "6 m/s (1 sigma)",
        "10.1051/0004-6361/202038937", "against %.1f m/s: %+.2f sigma" % (t["predicted"], t["z"]))
    return C_


def main():
    e2r = experiments2(); e3r = experiments3(); rp = rope(); pr = proof(e2r, e3r)
    out = {"experiments_2": e2r, "experiments_3": e3r, "rope": rp, "proof": pr}
    if "--json" in sys.argv:
        print(json.dumps(out, indent=1, sort_keys=True, default=float))
        return
    print("REPRODUCE: every number on the site, recomputed from data/")
    for k in ("E1", "E2", "E3", "E4", "E5", "E6", "E7", "E8", "E9"):
        for t in e2r[k]["tests"]:
            print("%s  %-62s residuals %-28s chi2 %-10.4g -> %s" % (k, t["label"][:62], t["residuals"], t["chi2"], t["result"]))
    print("E5  bound %s -> %s" % (e2r["E5"]["bound"]["label"], e2r["E5"]["bound"]["result"]))
    print("E3-part-A  scored %d, agree %d, contradicted %d, open %d, recalled %d" % (
        e3r["scored"], e3r["agree"], e3r["contradicted"], e3r["open"], e3r["recalled"]))
    for r in e3r["rows"]:
        print("  %-4s %-9s %s" % (r["id"], ("%+.3g" % r["z"]) if r.get("z") is not None else "-", r["result"]))
    print("WAVE  S: two ends, probabilistic %.4f; two ends, sign %.4f; read once for both ends %.4f; ceiling %.4f" % (
        rp["two_ends_probabilistic"], rp["two_ends_sign"], rp["read_once_for_both_ends"], rp["tsirelson"]))
    for c in pr:
        print("%s  %s  [%s]" % (c["id"], c["computed"], c["doi"]))


if __name__ == "__main__":
    main()