"""
Out-of-sample test of Davidson, U-Yen & Holloman (2015): do M8+ earthquakes after January 2014
fall inside the paper's 'Significant Windows' more often than chance?

Needs, in the same folder as this script and spf_windows.py:
  wso_polar.html  - saved from http://wso.stanford.edu/Polar.html
  usgs_m75.csv    - saved from https://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=1976-01-01&endtime=2026-10-06&minmagnitude=7.5&orderby=time-asc
Rules were frozen ({'c13': 'ge9'}) after checking the rebuild against the paper's 1976-2014 statistics,
before any post-2014 earthquake was scored.
"""
import datetime as dt, numpy as np
from math import comb
from spf_windows import *


def sf(k, n, p): return sum(comb(n, i) * p**i * (1-p)**(n-i) for i in range(k, n+1))
def pmf(k, n, p): return comb(n, k) * p**k * (1-p)**(n-k)
FROZEN = {'c13': 'ge9'}
quakes = load_quakes('usgs_m75.csv')

def run(period, data_end, opts=FROZEN, offset=0, report=False):
    day0, day1 = period
    dates, N, S, nint = load_wso('wso_polar.html', end=data_end)
    dates = np.array([d + dt.timedelta(days=offset) for d in dates])
    inst, ph, poles, _ = window_instances(dates, N, S, opts)
    events = delay_events(quakes)
    cov = build_coverage(dates, inst, events, day0, day1)
    tgt = [q for q in quakes if q[1] >= 8.0 and day0 <= q[0] <= day1]
    h = hits(cov, day0, tgt)
    if report:
        print(f"  WSO data used: {dates[0]} .. {dates[-1]} (interpolated values: {nint})")
        for d, m, p, ok in h:
            labs = [f"{lab} centred {dates[r]} ±{hw}d" for r, hw, lab in inst
                    if abs((d - dates[r]).days) <= hw + 20]
            print(f"   {d}  M{m:<4} {'INSIDE ' if ok else 'outside'}  {p[:38]:38s} | nearest: " + "; ".join(labs[:2]))
    return cov, h, inst, dates, events

def circular_shift_p(cov_builder, n_days, tgt_offsets, k_obs):
    """Shift all quakes together by every possible offset (wrap-around); count hits."""
    ks = []
    for o in range(n_days):
        cov = cov_builder(o)
        ks.append(sum(cov[(t + o) % n_days] for t in tgt_offsets))
    ks = np.array(ks)
    return (ks >= k_obs).mean(), ks

def shifted_cov_builder(inst, dates, day0, day1, o):
    n = (day1 - day0).days + 1
    qs = []
    for d, m, p in quakes:
        if day0 <= d <= day1:
            qs.append((day0 + dt.timedelta(days=((d - day0).days + o) % n), m, p))
    ev = delay_events(qs)
    return build_coverage(dates, inst, ev, day0, day1)

print("=== OUT-OF-SAMPLE (frozen rules, all WSO data to its last value) ===")
P_OUT = (dt.date(2014, 1, 18), dt.date(2026, 1, 9))
cov, h, inst, dates, events = run(P_OUT, None, report=True)
p_out = cov.mean(); n = len(h); k = sum(x[3] for x in h)
print(f"\n  Days inside windows, Jan 2014 - Jan 2026: {p_out*100:.1f}%")
print(f"  M8+ quakes inside windows: {k} of {n}   (chance expectation {n*p_out:.1f})")
print(f"  Binomial P(>= {k} | chance) = {sf(k, n, p_out):.3f}")
print(f"  P(exactly {k} | claimed 78.8% rate) = {pmf(k, n, 0.788):.4f};  P(<= {k} | 78.8%) = {1-sf(k+1, n, 0.788):.4f}")
lr = (0.788/p_out)**k * ((1-0.788)/(1-p_out))**(n-k)
print(f"  Likelihood ratio, claimed rate vs chance: {lr:.3g}  (i.e. {1/lr:.0f}:1 in favour of chance)" if lr < 1 else f"  LR {lr:.2f}")

day0, day1 = P_OUT; nd = (day1 - day0).days + 1
tgt_off = [(d - day0).days for d, m, p, ok in h]
pc, ks = circular_shift_p(lambda o: shifted_cov_builder(inst, dates, day0, day1, o), nd, tgt_off, k)
print(f"  Clustering-preserving test (all {nd} circular shifts): P(>= {k}) = {pc:.3f}; mean hits under shifts {ks.mean():.2f}")

print("\n=== IN-SAMPLE, same frozen rules, data to Jan 2014 (for comparison) ===")
P_IN = (dt.date(1976, 7, 5), dt.date(2014, 1, 17))
cov_i, h_i, inst_i, dates_i, _ = run(P_IN, dt.date(2014, 1, 17))
k_i = sum(x[3] for x in h_i); n_i = len(h_i)
print(f"  coverage {cov_i.mean()*100:.1f}%  hits {k_i}/{n_i}  binomial p = {sf(k_i, n_i, cov_i.mean()):.4f}")
d0, d1 = P_IN; ndi = (d1 - d0).days + 1
toff = [(d - d0).days for d, m, p, ok in h_i]
pci, ksi = circular_shift_p(lambda o: shifted_cov_builder(inst_i, dates_i, d0, d1, o), ndi, toff, k_i)
print(f"  clustering-preserving test: P(>= {k_i}) = {pci:.4f}; mean hits under shifts {ksi.mean():.2f}")

print("\n=== Robustness: other readings of the ambiguous rules (out-of-sample) ===")
for c13 in ('ge9', 'le9', 'none'):
    for off in (0, 4):
        cv, hh, *_ = run(P_OUT, None, opts={'c13': c13}, offset=off)
        kk = sum(x[3] for x in hh)
        print(f"  C13={c13:4s} offset+{off}: coverage {cv.mean()*100:5.1f}%  hits {kk}/{len(hh)}  p={sf(kk, len(hh), cv.mean()):.3f}")
np.save('outsample_cov.npy', cov)