"""Out-of-sample test of Nikitopoulos (2026). Rules fixed 2026-10-28 16:20 UTC, before any
planetary date or earthquake day outside 2010-2012 was computed or looked at.

1. Blocks. Primary: 2013-2025 (13 years, not used by the paper). Secondary: 2000-2009.
   Each block is tested on its own and never pooled with 2010-2012.
2. Planetary dates: select_dates.py / config_dates.py exactly as in the in-sample
   reproduction (Horizons DE441 heliocentric ecliptic ICRF states, TDB, hourly, cubic
   Hermite to every minute; the Sun first, then Mercury, Venus, Earth, Mars; Gram matrix
   and SVD of section 2.2; candidate where max|dX*| >= 0.35e-4 AU; while a year has more
   than 100 candidates add the Moon (after Earth), then Jupiter, then, as our assumption
   since the paper stops at Jupiter, Saturn, Uranus, Neptune). A year with no body set at
   or under 100 candidates contributes no dates. Dates are TDB calendar dates.
3. Earthquake days: USGS ComCat (fetch_usgs.py, minmagnitude=6.5), type "earthquake",
   preferred magnitude strictly > 6.5, UTC calendar days.
4. Statistic: same-day overlap K. Null: circular shift of the planetary dates over the whole
   block, earthquake days fixed; exact (every shift once); one-sided p = share of shifts
   with K >= observed.
5. Decision: replicated if the 2013-2025 block gives p < 0.05. Report the effect size
   (observed / expected; in sample 17 / 4.86 = 3.5) and the power: P(p < 0.05) if each date
   fell on an earthquake day at the in-sample rate 17/48, and at 1.5 x the chance rate
   (binomial count over the block's dates against the exact null's critical value).
6. Secondary, reported whatever they show: (a) per-year p-values in 2013-2025; (b) the
   2000-2009 block; (c) in sample 2010-2012, the whole procedure rerun with one arbitrary
   implementation choice changed (Sun's row last; bodies reversed after the Sun; eigh for
   the SVD; alpha = 1), combined test for each; (d) alignment: mean rank of R2 (four inner
   planets on one line through the Sun) on the dates against every day of the same year.

Standard library only. Usage:
  python3 -I outsample_test.py <results dir> <usgs csv> <dates json> [<diagnose json>]
"""
import json, math, pathlib, sys
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import datetime
import overlap_test as ot

IN_SAMPLE_RATE = 17 / 48


def binom_sf(k, n, p):
    """P(X >= k) for X ~ Binomial(n, p)."""
    return sum(math.comb(n, i) * p**i * (1 - p)**(n - i) for i in range(k, n + 1))


def block(dates, csv_path):
    years = sorted(int(y) for y in dates)
    B_all = ot.quake_days([csv_path], set(years))
    start = datetime.date(years[0], 1, 1); L = (datetime.date(years[-1] + 1, 1, 1) - start).days
    A = {ot.index(d, start) for y in years for d in dates[str(y)]}
    B = {ot.index(d, start) for d in B_all}
    res = ot.test(A, B, L, draws=1)        # exact null only (one Monte Carlo draw is a placeholder)
    # exact null distribution again, for the critical value
    Bset = set(B)
    dist = sorted(sum(((a + s - 1) % L) + 1 in Bset for a in A) for s in range(L))
    kcrit = next(k for k in range(0, len(A) + 2) if sum(x >= k for x in dist) / L < 0.05)
    m = len(A); chance = res["expected_exact"] / m if m else 0
    res.update({"years": f"{years[0]}-{years[-1]}", "ratio_obs_exp": round(res["observed"] / res["expected_exact"], 3) if m else None,
                "k_crit_p_below_0.05": kcrit,
                "power_at_in_sample_rate": round(binom_sf(kcrit, m, IN_SAMPLE_RATE), 4),
                "power_at_1.5x_chance": round(binom_sf(kcrit, m, 1.5 * chance), 4),
                "chance_rate_per_date": round(chance, 4)})
    per_year = {}
    for y in years:
        ys = datetime.date(y, 1, 1); Ly = (datetime.date(y + 1, 1, 1) - ys).days
        Ay = {ot.index(d, ys) for d in dates[str(y)]}
        By = {ot.index(d, ys) for d in B_all if d.startswith(str(y))}
        r = ot.test(Ay, By, Ly, draws=1)
        per_year[str(y)] = {k: r[k] for k in ("m_dates", "n_quake_days", "observed", "expected_exact", "p_exact")}
    res["per_year"] = per_year
    res["years_p_below_0.05"] = sum(v["p_exact"] < 0.05 for v in per_year.values())
    res.pop("p_mc", None); res.pop("expected_mc", None)
    return res


def alignment(diag_path):
    diag = json.load(open(diag_path, encoding="utf-8"))
    by_date = {}
    for o in diag:
        by_date.setdefault(o["stamp"][:10], o["R2_pct_of_year"])
    v = list(by_date.values())
    n = len(v); mean = sum(v) / n
    return {"dates": n, "mean_rank_R2": round(mean, 3), "se_if_unrelated": round(math.sqrt(1 / 12 / n), 3),
            "share_in_top_tenth": round(sum(x >= 0.9 for x in v) / n, 3)}


def main():
    out_dir = pathlib.Path(sys.argv[1]); csv_path = sys.argv[2]; dates = json.load(open(sys.argv[3], encoding="utf-8"))
    res = block(dates, csv_path)
    if len(sys.argv) > 4:
        res["alignment"] = alignment(sys.argv[4])
    tag = res["years"]
    (out_dir / f"outsample-{tag}.json").write_text(json.dumps(res, indent=1), encoding="utf-8")
    print(json.dumps(res, indent=1))


if __name__ == "__main__":
    main()