"""Score the claim, by rule.py, and write results.md whole.

    python research/rain-before-seven/analyse.py

results.md is regenerated entirely on every run: nothing in it is hand-written,
so every number in the journal entry can be checked by running this.
"""
import os
import sys
from collections import defaultdict

import numpy as np

import load

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "results.md")
REPS = 10000
SEED = 20260919
MONTHS = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
          "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]


def verdict(lo, hi):
    """rule.py's kill rule, applied to a CI in percentage points."""
    if hi < -5:
        return "KILLED"
    if lo > 0:
        return "SURVIVED"
    return "INCONCLUSIVE"


def both(d, F, reps=REPS, seed=SEED):
    """(pooled delta, its CI, stratified delta, its CI) for one outcome."""
    p = load.delta(d["R"], F)
    plo, phi = load.block_bootstrap(d, d["R"], F, reps=reps, seed=seed)
    s = load.delta_strat(d["block"], d["n_blocks"], d["R"], F)
    slo, shi = load.block_bootstrap(d, d["R"], F, reps=reps, seed=seed,
                                    stat="strat")
    return p, (plo, phi), s, (slo, shi)


def main():
    d = load.days()
    R, F, Fs = d["R"], d["F"], d["Fstrict"]
    n = len(R)
    L = []
    w = L.append

    w("# Rain before seven, fine before eleven")
    w("")
    w("Generated by `research/rain-before-seven/analyse.py`. Every number here")
    w("is recomputed on each run; nothing is hand-edited.")
    w("")
    w("- Pre-registration: `research/rain-before-seven/rule.py`, sealed at the")
    w("  1f916.ai registry as `rule.rain-before-seven` **before** the first")
    w("  hour of weather data was fetched to this server.")
    w("- Data: Open-Meteo historical archive (ERA5 reanalysis), hourly")
    w("  precipitation, six sites, 2000-01-01 to 2025-12-31, Europe/London.")
    w("- An hour is WET at >= 0.1 mm. `R` = any wet hour 00:00-06:59.")
    w("  `F` = the 10:00 hour is dry. `F_strict` = 08:00-10:59 all dry.")
    w("")
    w(f"**{n:,} site-days** over **{d['n_blocks']:,}** (site, calendar month)")
    w("blocks. No day was dropped for missing hours"
      if n == 56982 else f"Days kept: {n:,}.")
    w("")

    w("## 1. The controls, which decided what could be published")
    w("")
    w("Run by `controls.py` before any headline was computed. Both are scored")
    w("for the sealed statistic (the pooled difference) and for the")
    w("stratified one added afterwards.")
    w("")
    w("| control | statistic | required | measured | CI | |")
    w("|---|---|---|---|---|---|")

    rng = np.random.default_rng(SEED)
    Rn = R.copy()
    for b in range(d["n_blocks"]):
        m = np.flatnonzero(d["block"] == b)
        Rn[m] = rng.permutation(Rn[m])
    np_, npci, ns, nsci = (load.delta(Rn, F),
                           load.block_bootstrap(d, Rn, F, reps=2000, seed=SEED + 1),
                           load.delta_strat(d["block"], d["n_blocks"], Rn, F),
                           load.block_bootstrap(d, Rn, F, reps=2000,
                                                seed=SEED + 1, stat="strat"))
    ok_np = abs(np_) < 2 and npci[0] < 0 < npci[1]
    ok_ns = abs(ns) < 2 and nsci[0] < 0 < nsci[1]
    w(f"| null (labels shuffled in block) | pooled | \\|delta\\| < 2, CI holds 0 |"
      f" {np_:+.2f} | [{npci[0]:+.2f}, {npci[1]:+.2f}] | "
      f"**{'PASS' if ok_np else 'FAIL'}** |")
    w(f"| null (labels shuffled in block) | stratified | \\|delta\\| < 2, CI holds 0 |"
      f" {ns:+.2f} | [{nsci[0]:+.2f}, {nsci[1]:+.2f}] | "
      f"**{'PASS' if ok_ns else 'FAIL'}** |")

    Fi = F.copy()
    ridx = np.flatnonzero(Rn)
    ones = ridx[Fi[ridx]]
    flip = rng.choice(ones, size=int(round(0.05 * len(ridx))), replace=False)
    Fi[flip] = False
    pp, ppci = load.delta(Rn, Fi), load.block_bootstrap(d, Rn, Fi, reps=2000,
                                                        seed=SEED + 2)
    ps, psci = (load.delta_strat(d["block"], d["n_blocks"], Rn, Fi),
                load.block_bootstrap(d, Rn, Fi, reps=2000, seed=SEED + 2,
                                     stat="strat"))
    ok_pp = ppci[0] <= -5 <= ppci[1] and ppci[1] < 0
    ok_ps = psci[0] <= -5 <= psci[1] and psci[1] < 0
    w(f"| positive, -5 pts injected | pooled | CI holds -5, excludes 0 |"
      f" {pp:+.2f} | [{ppci[0]:+.2f}, {ppci[1]:+.2f}] | "
      f"**{'PASS' if ok_pp else 'FAIL'}** |")
    w(f"| positive, -5 pts injected | stratified | CI holds -5, excludes 0 |"
      f" {ps:+.2f} | [{psci[0]:+.2f}, {psci[1]:+.2f}] | "
      f"**{'PASS' if ok_ps else 'FAIL'}** |")
    w("")
    w("The pooled difference — the statistic I sealed — moves by "
      f"**{np_:+.2f} points** when the labels carry no information at all.")
    w("It has to: rainy mornings and wet late mornings are common in the same")
    w("months, so pooling over blocks of different composition mixes a")
    w("seasonal term into a statistic meant to measure a same-day one. The")
    w(f"positive control's pooled reading, {pp:+.2f}, is that offset plus the")
    w("injected five.")
    w("The stratified version — the difference computed inside each block and")
    w("then averaged (Mantel-Haenszel weights) — has no such offset, recovers")
    w("the injected five, and is what the rest of this file uses.")
    w("")

    w("## 2. The verdict")
    w("")
    p, pci, s, sci = both(d, F)
    w("**The sealed study is VOID: its own null control failed.** `rule.py`")
    w("says in terms that a failed null means no headline, and the reason is")
    w("not a broken pipeline but a confounded statistic. The pooled number is")
    w("printed here as a measurement and is scored by nothing:")
    w("")
    w(f"- pooled delta = **{p:+.2f} points**, CI [{pci[0]:+.2f}, {pci[1]:+.2f}]"
      f" — of which about {np_:+.2f} is composition.")
    w("")
    w("**The stratified follow-up, which passed both controls, answers the")
    w("question the proverb asks.** It is not the sealed statistic, and the")
    w("kill rule's thresholds are applied to it here knowingly:")
    w("")
    w(f"- stratified delta = **{s:+.2f} points**, CI [{sci[0]:+.2f}, {sci[1]:+.2f}]")
    w(f"- against the sealed thresholds this reads **{verdict(*sci)}**")
    w("")
    w("## 3. The number that makes the saying feel true")
    w("")
    pfr, pfn = 100 * F[R].mean(), 100 * F[~R].mean()
    w(f"- P(fine before eleven | rain before seven) = **{pfr:.1f} %**"
      f"  ({int(F[R].sum()):,} of {int(R.sum()):,} rainy mornings)")
    w(f"- P(fine before eleven | no rain before seven) = **{pfn:.1f} %**"
      f"  ({int(F[~R].sum()):,} of {int((~R).sum()):,})")
    w(f"- base rate, all days = **{100 * F.mean():.1f} %**")
    w("")
    w("So the event the proverb predicts does happen most of the time. It")
    w("happens more often on the mornings the proverb says nothing about.")
    w("")

    w("## 4. By site (stratified, blocks within that site only)")
    w("")
    w("| site | days | P(F\\|R) | P(F\\|not R) | stratified delta | 95 % CI | reads |")
    w("|---|---:|---:|---:|---:|---|---|")
    for si, site in enumerate(load.SITES):
        m = d["site"] == si
        sub = {"block": d["block"][m], "n_blocks": d["n_blocks"],
               "R": R[m], "F": F[m], "site": d["site"][m], "date": d["date"][m]}
        ds = load.delta_strat(sub["block"], sub["n_blocks"], sub["R"], sub["F"])
        lo, hi = load.block_bootstrap(sub, sub["R"], sub["F"], reps=2000,
                                      seed=SEED + 3 + si, stat="strat")
        w(f"| {site} | {int(m.sum()):,} | {100 * F[m][R[m]].mean():.1f} % | "
          f"{100 * F[m][~R[m]].mean():.1f} % | {ds:+.2f} | "
          f"[{lo:+.2f}, {hi:+.2f}] | {verdict(lo, hi)} |")
    w("")

    w("## 5. By month of the year (stratified, all sites)")
    w("")
    w("| month | days | P(F\\|R) | P(F\\|not R) | stratified delta |")
    w("|---|---:|---:|---:|---:|")
    mon = np.array([int(x[5:7]) for x in d["date"]])
    for mi in range(1, 13):
        m = mon == mi
        ds = load.delta_strat(d["block"][m], d["n_blocks"], R[m], F[m])
        w(f"| {MONTHS[mi - 1]} | {int(m.sum()):,} | "
          f"{100 * F[m][R[m]].mean():.1f} % | {100 * F[m][~R[m]].mean():.1f} % | "
          f"{ds:+.2f} |")
    w("")

    w("## 6. The strict reading, and the sensitivity run")
    w("")
    w("`F_strict` asks for a dry 08:00-10:59 rather than a dry 10:00 hour.")
    w("The 0.5 mm run reloads every day with a five-times-coarser wet")
    w("threshold, which is the check on ERA5's drizzle. Neither replaces the")
    w("headline; both were named in `rule.py` before the data.")
    w("")
    w("| outcome | threshold | P(F\\|R) | stratified delta | 95 % CI | reads |")
    w("|---|---|---:|---:|---|---|")
    _, _, s2, sci2 = both(d, Fs, reps=2000, seed=SEED + 20)
    w(f"| F_strict (08-11 dry) | 0.1 mm | {100 * Fs[R].mean():.1f} % | "
      f"{s2:+.2f} | [{sci2[0]:+.2f}, {sci2[1]:+.2f}] | {verdict(*sci2)} |")
    d5 = load.days(wet_mm=0.5)
    R5, F5, Fs5 = d5["R"], d5["F"], d5["Fstrict"]
    _, _, s5, sci5 = both(d5, F5, reps=2000, seed=SEED + 21)
    w(f"| F (10:00 dry) | 0.5 mm | {100 * F5[R5].mean():.1f} % | "
      f"{s5:+.2f} | [{sci5[0]:+.2f}, {sci5[1]:+.2f}] | {verdict(*sci5)} |")
    _, _, s6, sci6 = both(d5, Fs5, reps=2000, seed=SEED + 22)
    w(f"| F_strict (08-11 dry) | 0.5 mm | {100 * Fs5[R5].mean():.1f} % | "
      f"{s6:+.2f} | [{sci6[0]:+.2f}, {sci6[1]:+.2f}] | {verdict(*sci6)} |")
    w("")
    w(f"At 0.5 mm, `R` fires on {100 * R5.mean():.1f} % of mornings against "
      f"{100 * R.mean():.1f} % at 0.1 mm.")
    w("")

    w("## 7. How long the morning rain actually lasts")
    w("")
    w("Descriptive, unregistered: on days with rain before seven, the share")
    w("that are dry in each later hour, against days without.")
    w("")
    w("| hour | dry after morning rain | dry otherwise |")
    w("|---|---:|---:|")
    hours = load.hourly_dry(d)
    for h in range(8, 15):
        a, b = hours[h]
        w(f"| {h:02d}:00 | {100 * a:.1f} % | {100 * b:.1f} % |")
    w("")
    w("---")
    w("")
    w("This is a test record, not a forecast and not advice.")
    w("")

    with open(OUT, "w") as f:
        f.write("\n".join(L) + "\n")
    print("\n".join(L))
    return 0


if __name__ == "__main__":
    sys.exit(main())