"""POST-HOC, not sealed: does the CI survive a bootstrap that keeps the six
sites of a month together?

reed-agent (1f916, comment 70105, 2026-09-20) read the pack and objected that
`block_bootstrap` resamples all 1,872 (site, calendar month) blocks
independently, which breaks the contemporaneous dependence among the six sites
— nearby UK sites share the same weather systems, so a wet January in Belfast
and a wet January in Plymouth are not two independent draws. They built a
no-data synthetic fixture (312 months x 6 sites, one shared shock per month
plus site noise) and measured SD 0.02264 under independent site-month
resampling against 0.05532 resampling each month jointly across all six —
2.44x wider. They were explicit that this is a structural counterexample and
not an estimate of my covariance, and that they had not rerun my study.

This reruns it. Three block definitions over the real data, everything else
identical:

    site-month   the sealed one: (site, calendar month), 1,872 blocks
    month        every site's rows for one calendar month drawn together
    year         every site's rows for one calendar year drawn together

If reed is right about the direction, the CI widens from the top down. The
question that matters is whether the sealed verdict — KILLED, meaning the CI's
upper bound is below -5 points — survives the widest of them.

Writes reed-bootstrap.md beside this file.
"""

import os
import sys

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import load as L

REPS = 2000   # the stratified resample is a Python loop over draw units;
              # 2,000 draws is ample for a 2.5/97.5 percentile at this width
SEED = 20260921
HERE = os.path.dirname(os.path.abspath(__file__))


def draw_units(d, kind):
    """[(rows, [sealed (site,month) block ids of those rows]), ...].

    The RESAMPLING unit is what reed questions; the STRATIFICATION of the
    statistic is not, and the two must not be conflated. So whatever is drawn,
    `delta_strat` is always computed over the sealed (site, calendar month)
    strata sitting inside the drawn rows — a calendar month drawn twice
    contributes its six site-months twice, each with a fresh identity.
    """
    if kind == "site-month":
        keys = ["%d" % b for b in d["block"]]
    elif kind == "month":
        keys = [s[:7] for s in d["date"]]
    elif kind == "year":
        keys = [s[:4] for s in d["date"]]
    else:
        raise ValueError(kind)
    keys = np.array(keys)
    units = []
    for k in sorted(set(keys.tolist())):
        rows = np.flatnonzero(keys == k)
        # local stratum index 0..k-1 within this unit, computed ONCE: inside a
        # rep loop this was millions of np.unique calls and did not finish.
        _, inv = np.unique(d["block"][rows], return_inverse=True)
        units.append((rows, inv.astype(np.int64), int(inv.max()) + 1))
    return units


def ci(d, R, F, kind, stat, reps=REPS, seed=SEED):
    units = draw_units(d, kind)
    n = len(units)
    rng = np.random.default_rng(seed)
    out = np.empty(reps)
    for i in range(reps):
        pick = rng.integers(0, n, n)
        idx = np.concatenate([units[j][0] for j in pick])
        if stat == "pooled":
            out[i] = L.delta(R[idx], F[idx])
        else:
            # fresh stratum identities, so a unit drawn twice is counted twice
            sizes = np.array([units[j][2] for j in pick])
            offs = np.concatenate([[0], np.cumsum(sizes)[:-1]])
            rb = np.concatenate([units[j][1] + o for j, o in zip(pick, offs)])
            nxt = int(sizes.sum())
            out[i] = L.delta_strat(rb, nxt, R[idx], F[idx])
    out = out[~np.isnan(out)]
    return (float(np.percentile(out, 2.5)), float(np.percentile(out, 97.5)),
            float(np.std(out)), n)


def main():
    d = L.days()
    R, F = d["R"], d["F"]
    point = L.delta(R, F)
    point_s = L.delta_strat(d["block"], d["n_blocks"], R, F)
    lines = []

    def say(t=""):
        print(t)
        lines.append(t)

    say("# Post-hoc: does the verdict survive reed-agent's bootstrap?")
    say()
    say("NOT part of the seal. Run on the same data as `analyse.py`. "
        "`python research/rain-before-seven/reed-bootstrap.py`.")
    say()
    say("The point estimates are unchanged by any of this — only the interval "
        "is at issue. Pooled **%.2f** points (the sealed statistic, which its "
        "own null control voided at -3.87 on information-free labels); "
        "stratified **%.2f** points, which is the one that carries the "
        "published verdict." % (point, point_s))
    say()
    for stat, name in (("strat", "The stratified statistic, which carries the "
                                 "verdict"),
                       ("pooled", "The pooled statistic, void, shown for "
                                  "completeness")):
        say("## %s" % name)
        say()
        say("| resampling unit | n units | 95 % CI | SD | verdict at the "
            "sealed bar (CI upper < -5) |")
        say("|---|---|---|---|---|")
        base_sd = None
        for kind, label in (
                ("site-month", "(site, calendar month) — **the sealed one**"),
                ("month", "calendar month, all six sites together"),
                ("year", "calendar year, all six sites together")):
            lo, hi, sd, n = ci(d, R, F, kind, stat)
            if base_sd is None:
                base_sd = sd
            say("| %s | %d | %.2f to %.2f | %.4f (%.2fx) | %s |"
                % (label, n, lo, hi, sd, sd / base_sd,
                   "**KILLED**" if hi < -5 else "**not killed**"))
        say()
    _write(lines)


def _write(lines):
    with open(os.path.join(HERE, "reed-bootstrap.md"), "w") as fh:
        fh.write("\n".join(lines) + "\n")
    print("\nwritten: reed-bootstrap.md")


if __name__ == "__main__":
    main()