"""Post-hoc, described only: is the busiest ten years unusual for a world with
a constant rate?

In both catalogs the busiest stretch is 2006-2015. Picking the busiest window
after looking is a multiple-comparisons trap, so the null keeps the picking:
simulate a constant-rate Poisson world at each catalog's own mean, take the
busiest ten-year window in each simulated record, and ask how often it is at
least as busy as the observed busiest window.

This asks about annual counts only. Aftershock sequences put several M7s in
one year, and the measured dispersion of the annual counts is 1.00, so a
Poisson null is the fair one at this resolution; it says nothing about
clustering within a year.

    python cluster.py > cluster.txt
"""

import csv
import pathlib

import numpy as np

HERE = pathlib.Path(__file__).parent
SIMS = 100_000
WINDOW = 10
rng = np.random.default_rng(20260911)


def busiest(counts):
    c = np.convolve(counts, np.ones(WINDOW, dtype=int), mode="valid")
    return int(c.max()), int(c.argmax())


def load():
    comcat, gcmt = {}, {}
    with (HERE / "yearly.csv").open() as f:
        for row in csv.DictReader(f):
            if float(row["minmag"]) == 7.0:
                comcat[int(row["year"])] = int(row["count"])
    with (HERE / "gcmt_yearly.csv").open() as f:
        for row in csv.DictReader(f):
            gcmt[int(row["year"])] = int(row["gcmt_mw7"])
    return comcat, gcmt


comcat, gcmt = load()
for name, series in (("ComCat M>=7", comcat), ("GCMT Mw>=7", gcmt)):
    years = sorted(series)
    counts = np.array([series[y] for y in years])
    top, at = busiest(counts)
    mean = counts.mean()
    sims = rng.poisson(mean, size=(SIMS, len(counts)))
    windows = np.lib.stride_tricks.sliding_window_view(sims, WINDOW, axis=1).sum(axis=2)
    p = float((windows.max(axis=1) >= top).mean())
    print(
        f"{name:<12} {years[0]}-{years[-1]}: busiest {WINDOW} years {years[at]}-{years[at] + WINDOW - 1} "
        f"= {top} ({top / WINDOW:.1f}/yr against a mean of {mean:.2f}); "
        f"P(busiest window >= {top} | constant Poisson rate) = {p:.3f}"
    )