"""Apply the committed rule to the fetched counts, then describe the rest.

    python analyse.py > results.txt

1. The sealed verdict: rule.py on yearly.csv at minmag 7.0, unchanged.
2. Described, no bar: the same trend at M >= 4.5, 5.0, 6.0.
3. Described, no bar: magnitude of completeness per window by maximum
   curvature (the 0.1-unit bin holding the most events), and the share of each
   window's M >= 4.5 events that sit below the worst window's completeness.
"""

import csv
import pathlib
from collections import defaultdict

from rule import trend, verdict

HERE = pathlib.Path(__file__).parent

yearly = defaultdict(dict)
with (HERE / "yearly.csv").open() as f:
    for row in csv.DictReader(f):
        yearly[float(row["minmag"])][int(row["year"])] = int(row["count"])

print("== 1. Sealed: M >= 7.0, 1973-2025 ==")
years = sorted(yearly[7.0])
assert years == list(range(1973, 2026)), "yearly.csv is incomplete at M7"
t = trend(years, [yearly[7.0][y] for y in years])
print(
    f"mean {t['mean_per_year']:.2f}/yr over {t['n_years']} years; "
    f"change per decade {t['change']:+.1%} (95 % CI {t['lo']:+.1%} to {t['hi']:+.1%}); "
    f"dispersion phi {t['phi']:.2f}"
)
print("VERDICT:", verdict(t))
print("M7 counts by year:", " ".join(f"{y}:{yearly[7.0][y]}" for y in years))

print()
print("== 2. Described: the same model at every threshold ==")
print(f"{'minmag':>6} {'mean/yr':>8} {'1973-82':>8} {'2016-25':>8} {'ratio':>6} {'per decade':>11} {'95 % CI':>18} {'phi':>6}")
for m in sorted(yearly):
    ys = sorted(yearly[m])
    if len(ys) != 53:
        print(f"{m:>6} incomplete ({len(ys)} years)")
        continue
    tt = trend(ys, [yearly[m][y] for y in ys])
    early = sum(yearly[m][y] for y in range(1973, 1983)) / 10
    late = sum(yearly[m][y] for y in range(2016, 2026)) / 10
    print(
        f"{m:>6} {tt['mean_per_year']:>8.1f} {early:>8.1f} {late:>8.1f} {late / early:>6.2f} "
        f"{tt['change']:>+11.1%} {tt['lo']:>+8.1%} to {tt['hi']:>+6.1%} {tt['phi']:>6.1f}"
    )

hist_path = HERE / "hist.csv"
if hist_path.exists():
    print()
    print("== 3. Described: completeness by window (maximum curvature) ==")
    hist = defaultdict(dict)
    with hist_path.open() as f:
        for row in csv.DictReader(f):
            hist[row["window"]][float(row["bin"])] = int(row["count"])
    mcs = {}
    for w in sorted(hist):
        bins = hist[w]
        mc = max(bins, key=lambda b: bins[b])
        mcs[w] = mc
        span = int(w[-4:]) - int(w[:4]) + 1
        total = sum(bins.values())
        print(f"{w}: Mc {mc:.1f}, {total / span:,.0f} events/yr in 2.5-7.4, peak bin {bins[mc]:,}")
    # Checked 2026-09-11 and it does not work here: in five of six windows
    # the histogram's mode is the 2.5 floor, because the global catalog is
    # dominated by dense regional networks (California, Alaska, ...) that
    # record far below global completeness. Maximum curvature on a mixed
    # catalog measures those networks, not the world. Printed so the failure
    # is visible; not reported as a finding.
    at_floor = sum(1 for mc in mcs.values() if mc == 2.5)
    print(f"NOT A RESULT: {at_floor} of {len(mcs)} windows put Mc at the 2.5 floor; "
          "the method needs a catalog of one network's coverage.")