"""The real statistic, run only after controls.py has passed.

    python analyse.py > results.txt

Everything it computes is defined in rule.py, which was sealed at the 1f916
registry on 2026-09-22 at 13:10:48Z, before the catalogue was fetched.
Writes results.txt here and curve.csv (events per degree of syzygy distance).
"""

import csv
import json
import pathlib

import numpy as np

import controls
import rule

HERE = pathlib.Path(__file__).parent
DATA = pathlib.Path.home() / "data" / "quake-moon"
rng = np.random.default_rng(rule.SEED + 1)


def report(name, theta, f_exp, distance=rule.syzygy_distance, bootstrap=True,
           cluster_id=None):
    e = rule.excess(theta, f_exp, distance=distance)
    n = len(theta)
    inside = int((distance(theta) <= rule.HALF).sum())
    if cluster_id is not None:
        lo, hi = rule.bootstrap_excess_by_cluster(theta, cluster_id, f_exp, rng,
                                                  draws=2000, distance=distance)
    elif bootstrap:
        lo, hi = rule.bootstrap_excess(theta, f_exp, rng, distance=distance)
    else:
        _, lo, hi = controls.ci_normal((distance(theta) <= rule.HALF).astype(float),
                                       f_exp)
    v = rule.verdict(lo, hi)
    print(f"{name}")
    print(f"  n = {n}   inside the syzygy half = {inside} "
          f"({100*inside/n:.2f} %)   expected {100*f_exp:.2f} %")
    print(f"  E = {100*e:+.2f} %   95 % [{100*lo:+.2f}, {100*hi:+.2f}]"
          f"   bar +{100*rule.BAR:.0f} %   -> {v}")
    return {"n": n, "inside": inside, "e": e, "lo": lo, "hi": hi, "verdict": v}


def main():
    times, lat, lon, dep, mag = controls.load()
    f_exp, f_full = np.load(DATA / "f_exp.npy")
    theta_all = rule.phases(times)
    mains, cluster = rule.decluster(times, lat, lon, mag)
    th = theta_all[mains]
    out = {}

    print("PRIMARY -- declustered M>=6.0, syzygy half, mainshock bootstrap\n")
    out["primary"] = report("primary", th, f_exp)

    # the shortcut used in the floor, checked here against the real bootstrap
    inside = (rule.syzygy_distance(th) <= rule.HALF).astype(float)
    ne, nlo, nhi = controls.ci_normal(inside, f_exp)
    print(f"  normal approximation used for the floor: "
          f"[{100*nlo:+.2f}, {100*nhi:+.2f}]   "
          f"max gap to the percentile bootstrap "
          f"{100*max(abs(nlo-out['primary']['lo']), abs(nhi-out['primary']['hi'])):.3f} pts")

    r, p = rule.schuster(th)
    print(f"\n  Schuster on 2*theta: R = {r:.2f}, n = {len(th)}, "
          f"p = {p:.4f}   (independence assumed; the verdict is not carried "
          f"by this)")
    out["schuster"] = {"R": r, "p": p, "n": len(th)}

    print("\nROBUSTNESS -- the full catalogue, whole clusters resampled\n")
    out["full"] = report("all events, cluster bootstrap", theta_all, f_exp,
                         cluster_id=cluster)

    print("\nSECONDARIES\n")
    shallow = mains[np.nan_to_num(dep[mains], nan=1e9) <= 70.0]
    out["S1_shallow"] = report("S1  shallow, depth <= 70 km",
                               theta_all[shallow], f_exp)
    m7 = mains[mag[mains] >= 7.0]
    out["S2_m7"] = report("S2  M >= 7.0  (floor said void: unfailable)",
                          theta_all[m7], f_exp)
    out["S3_full"] = report("S3  full moon alone", th, f_full,
                            distance=rule.full_distance)
    mains_s, _ = rule.decluster(times, lat, lon, mag, symmetric=True)
    out["S4_symmetric"] = report("S4  symmetric declustering window",
                                 theta_all[mains_s], f_exp)

    # the shape, for the chart: events per 5 degrees of syzygy distance,
    # against the time spent there
    th_hour = controls.grid_phases(3600, "hour")
    edges = np.arange(0, 95, 5)
    d_ev = rule.syzygy_distance(th)
    d_tm = rule.syzygy_distance(th_hour)
    obs, _ = np.histogram(d_ev, bins=edges)
    tm, _ = np.histogram(d_tm, bins=edges)
    exp = tm / tm.sum() * obs.sum()
    with (HERE / "curve.csv").open("w", newline="") as f:
        w = csv.writer(f)
        w.writerow(["d_lo", "d_hi", "observed", "expected", "ratio"])
        for i in range(len(obs)):
            w.writerow([edges[i], edges[i + 1], int(obs[i]), round(float(exp[i]), 2),
                        round(float(obs[i] / exp[i]), 4)])
    print("\nshape by syzygy distance (5-degree bins), observed / expected:")
    for i in range(len(obs)):
        print(f"  {edges[i]:>2}-{edges[i+1]:>2} deg  {obs[i]:>5}  "
              f"{exp[i]:>8.1f}  {obs[i]/exp[i]:.3f}")

    out["f_exp"] = float(f_exp)
    out["f_full"] = float(f_full)
    (HERE / "results.json").write_text(json.dumps(out, indent=1) + "\n")


if __name__ == "__main__":
    main()