"""Same-day overlap between "planetary configuration dates" and earthquake days, with the
preprint's circular-shift null (section 2.4, pp. 6-8). Standard library only.

Earthquake days: UTC calendar days with at least one USGS ComCat event of preferred
magnitude > 6.5 (strictly greater, which reproduces the paper's 39, 38 and 34 days).
The planetary dates are shifted circularly within the analysis interval of L days;
the earthquake days stay fixed. We give the exact null (every one of the L shifts once)
and the paper's Monte Carlo version (200,000 uniform shifts, seeded).

Usage: python3 -I overlap_test.py <dates.json> <usgs.csv> [<usgs.csv> ...]
dates.json: {"2010": ["2010-01-06", ...], ...}; years must be consecutive for the combined test.
"""
import csv, datetime, json, random, sys


def quake_days(paths, years, mmin=6.5):
    days = set()
    for p in paths:
        with open(p, encoding="utf-8") as f:
            for r in csv.DictReader(f):
                if r["type"] != "earthquake" or not r["mag"]:
                    continue
                if float(r["mag"]) > mmin and int(r["time"][:4]) in years:
                    days.add(r["time"][:10])
    return days


def index(day, start):
    return (datetime.date.fromisoformat(day) - start).days + 1      # 1..L


def test(A, B, L, draws=200_000, seed=20261028):
    """A, B: sets of day indices in 1..L. Returns observed, exact and Monte Carlo statistics."""
    Bset = set(B)
    def overlap(delta):
        return sum(((a + delta - 1) % L) + 1 in Bset for a in A)
    k_obs = overlap(0)
    dist = [overlap(d) for d in range(L)]                 # every shift once: the exact null
    exact_p = sum(k >= k_obs for k in dist) / L
    exact_mean = sum(dist) / L
    srt = sorted(dist); q95 = srt[int(0.95 * L) - 1] if L >= 20 else None
    rng = random.Random(seed)
    mc = [dist[rng.randrange(L)] for _ in range(draws)]   # same as drawing shifts uniformly
    mc_p = sum(k >= k_obs for k in mc) / draws
    mc_mean = sum(mc) / draws
    shifts_ge = sum(k >= k_obs for k in dist)
    return {"L": L, "m_dates": len(A), "n_quake_days": len(Bset), "observed": k_obs,
            "expected_exact": round(exact_mean, 4), "expected_mc": round(mc_mean, 4),
            "p_exact": round(exact_p, 5), "shifts_at_or_above_observed": shifts_ge,
            "p_mc": round(mc_p, 5), "q95_exact": q95, "max_over_shifts": max(dist)}


def main():
    dates = json.load(open(sys.argv[1], encoding="utf-8"))
    years = sorted(int(y) for y in dates)
    B_all = quake_days(sys.argv[2:], set(years))
    out = {}
    for y in years:
        start = datetime.date(y, 1, 1); L = (datetime.date(y + 1, 1, 1) - start).days
        A = {index(d, start) for d in dates[str(y)]}
        B = {index(d, start) for d in B_all if d.startswith(str(y))}
        r = test(A, B, L)
        r["overlap_days"] = sorted(d for d in dates[str(y)] if d in B_all)
        out[str(y)] = r
    if len(years) > 1:
        start = datetime.date(years[0], 1, 1); L = (datetime.date(years[-1] + 1, 1, 1) - start).days
        A = {index(d, start) for y in years for d in dates[str(y)]}
        B = {index(d, start) for d in B_all}
        out[f"{years[0]}-{years[-1]}"] = test(A, B, L)
    print(json.dumps(out, indent=1))


if __name__ == "__main__":
    main()