"""The mate, written AFTER the sealed verdict came back `survived`.

Question it answers: is the rise in ComCat's M >= 7 count a rise in the Earth,
or a change of ruler? ComCat sized 72 of its 78 M >= 7 events of the 1970s on
Ms (surface-wave magnitude), which under-reads deep earthquakes, and every one
from the 1990s on moment magnitude (Mw).

Second surface: the Global CMT catalog, which computes Mw from the seismic
moment for every event it solves, 1976 onward, with one method throughout.
The same `rule.py`, unchanged, is applied to both catalogs on the years both
cover: 1976-2020 (the published table) and 1976-2025 (added the same session,
once the monthly files for 2021-2025 were fetched). This is a post-hoc check
and is reported as one; it does not re-file the sealed verdict.

    python mate.py > mate.txt   # also writes gcmt_yearly.csv

Sources, kept in ~/data, not committed (Dziewonski et al. 1981; Ekström et al. 2012):
  https://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/jan76_dec20.ndk
  https://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/NEW_MONTHLY/<yyyy>/<mon><yy>.ndk
NDK: five lines an event. Line 1 hypocentre (year at [5:9]), line 3 centroid
(depth at [48:53]), line 4 exponent at [0:2], line 5 scalar moment at [49:56]
in 10**exponent dyne-cm. Mw = 2/3 (log10 M0 - 16.1).
"""

import csv
import math
import pathlib
from collections import Counter

import pandas as pd

from rule import trend, verdict

HERE = pathlib.Path(__file__).parent
DATA = pathlib.Path.home() / "data"
FILES = [DATA / "gcmt_jan76_dec20.ndk", *sorted((DATA / "gcmt_monthly").glob("*.ndk"))]
WINDOWS = {"1976-2020": list(range(1976, 2021)), "1976-2025": list(range(1976, 2026))}

m7, deep7, total = Counter(), Counter(), 0
for path in FILES:
    lines = path.read_text().splitlines()
    assert len(lines) % 5 == 0, f"{path} is not in five-line records"
    for i in range(0, len(lines), 5):
        rec = lines[i : i + 5]
        year = int(rec[0][5:9])
        exponent = int(rec[3][0:2])
        moment = float(rec[4][49:56])
        mw = round((2.0 / 3.0) * (math.log10(moment) + exponent - 16.1), 1)
        depth = float(rec[2][48:53])
        total += 1
        if mw >= 7.0:
            m7[year] += 1
            if depth > 300:
                deep7[year] += 1

print(f"GCMT files: {len(FILES)}; records parsed: {total:,}; Mw >= 7.0: {sum(m7.values())}")
years_seen = sorted(m7)
assert years_seen[0] == 1976 and years_seen[-1] == 2025, years_seen[-1]

with (HERE / "gcmt_yearly.csv").open("w", newline="") as f:
    w = csv.writer(f)
    w.writerow(["year", "gcmt_mw7", "gcmt_mw7_deep300"])
    for y in WINDOWS["1976-2025"]:
        w.writerow([y, m7[y], deep7[y]])

yearly = pd.read_csv(HERE / "yearly.csv")
comcat = yearly[yearly.minmag == 7.0].set_index("year")["count"]

events = pd.read_csv(DATA / "comcat_m65.csv")
events["year"] = pd.to_datetime(events["time"]).dt.year
cc_deep = events[(events.mag >= 7.0) & (events.depth > 300)].groupby("year").size()

for label, years in WINDOWS.items():
    print()
    print(f"== The sealed rule, unchanged, on {label} in both catalogs ==")
    for name, series in (("ComCat M>=7", [comcat[y] for y in years]), ("GCMT Mw>=7", [m7[y] for y in years])):
        t = trend(years, series)
        print(
            f"{name:<12} mean {t['mean_per_year']:5.2f}/yr  per decade {t['change']:+6.1%} "
            f"(95 % CI {t['lo']:+.1%} to {t['hi']:+.1%})  phi {t['phi']:.2f}  -> {verdict(t)}"
        )

print()
print("== Five-year sums: ComCat, GCMT, and deep (>300 km) events in each ==")
print(f"{'years':<10} {'ComCat':>7} {'GCMT':>6} {'CC deep':>8} {'GCMT deep':>10}")
for a in range(1976, 2026, 5):
    span = range(a, a + 5)
    print(
        f"{a}-{span[-1]:<5} {sum(comcat[y] for y in span):>7} {sum(m7[y] for y in span):>6} "
        f"{sum(int(cc_deep.get(y, 0)) for y in span):>8} {sum(deep7[y] for y in span):>10}"
    )