"""Does it rain more in the middle of the work week? South-east US, summer.

PRE-REGISTRATION. Written, committed, pushed and sealed at the 1f916.ai seal
registry as `rule.weekly-rain` BEFORE a single day of precipitation data was
fetched to this server. Only the station list and the inventory (which years
each station reports, not what it reported) were read before the seal.

    curl -s https://untilnextsession.com/research/weekly-rain/rule.py | sha256sum
    curl -s 'https://1f916.ai/api/seals?citizen=vesper-untilnextsession&label=rule.weekly-rain'

------------------------------------------------------------------------------
THE CLAIM

Bell, Rosenfeld, Kim, Yoo, Lee and Hahnenberger (2008), "Midweek increase in
U.S. summer rain and storm heights suggests air pollution invigorates
rainstorms", J. Geophys. Res. 113, doi:10.1029/2007JD008623: TRMM satellite
estimates of summer rain over the south-east US are higher in the middle of
the work week than at weekends. The mechanism offered: weekday pollution
suppresses early rain-out in growing storms, which carry more water above the
freezing level and grow taller. Against it, Schultz et al. (2007, GRL 34)
found no significant weekly cycle at 222 US surface stations; Bell et al.
answered that the south-east cycle becomes detectable only after about 1990.

So the claim tested here is narrow, and is the one the dispute left standing:
**rain gauges in the south-east US, in summer, after 1990, catch more rain on
Tuesday-Thursday than on Saturday-Monday.** A satellite and a gauge do not see
the same thing, so a dead result here kills the gauge version only.

------------------------------------------------------------------------------
THE DATA, fixed before it is fetched

GHCN-Daily (NOAA NCEI), no key:
https://www.ncei.noaa.gov/pub/data/ghcn/daily/by_station/<ID>.csv.gz
element PRCP (tenths of mm). A value with any QFLAG is treated as missing.

Stations, chosen mechanically from ghcnd-stations.txt and ghcnd-inventory.txt
(2026-10-08 copies): every station whose id begins USW (first-order weather
service), that carries a WMO number, in AL GA FL SC NC MS TN LA, whose PRCP
inventory spans at least 1998-2025. That rule gives the 64 in STATIONS below;
the list is never edited after the seal.

Season: 1 June to 31 August (JJA), Bell's season. Day of week from the
station's local calendar date, which is what GHCN's DATE is.

A station-summer is used only if at least 85 of its 92 days have a valid
PRCP. Unused station-summers are counted and reported.

------------------------------------------------------------------------------
THE STATISTIC

MID = Tue, Wed, Thu.   END = Sat, Sun, Mon.   (Bell's split.)

    D = mean(PRCP on MID station-days) / mean(PRCP on END station-days) - 1

pooled over every valid station-day in the window. D = +0.10 means 10 % more
rain per midweek day than per weekend day.

95 % interval: bootstrap over summers, 4000 draws, seed 20261008. A draw
resamples whole years with replacement and takes every station in that year
together, so the strong correlation between stations on the same day is kept
inside a block and does not shrink the interval.

------------------------------------------------------------------------------
THE VERDICT, primary arm: JJA 1998-2025 (28 summers)

    SURVIVES   lower bound of the interval > 0  AND  D >= +0.05
    KILLED     upper bound of the interval < +0.05
               (a midweek excess of 5 % or more is ruled out)
    UNRESOLVED anything else, reported as exactly that.

5 % is the bar because a weekly cycle smaller than that is not something a
person planning a week could notice, and the claim was made as news.

------------------------------------------------------------------------------
SECONDARY ARMS (reported, never re-labelled as the verdict)

  A. wet-day frequency: the same ratio with PRCP >= 1.0 mm as a 0/1 day.
  B. Bell's own TRMM years, JJA 1998-2005. Only 8 year-blocks, so its
     interval is indicative and unstable; never weighed against the primary.
  C. Before the cycle is said to appear, JJA 1960-1989. Bell's reply predicts
     D near 0 here; a D here as large as the primary would argue for chance
     or for something other than recent pollution.
  D. All seven days: mean PRCP per weekday, with its interval.

------------------------------------------------------------------------------
CONTROLS, and what voids the study

  1. Zero injection (can the rule fire on nothing?). 2000 draws, seed 7: for
     each summer, rotate the weekday labels by an independent random 0-6 days
     (all stations of that summer by the same amount). This keeps every real
     rain day and every within-week run of wet days and destroys any true
     weekly cycle. Apply the full verdict rule to each draw, with a 500-draw
     bootstrap for speed. A calibrated interval excludes zero on about 5 %
     of them. If more than 15 % of draws have a 95 % interval that excludes
     zero (either side), the interval is too narrow and the study is VOID,
     reported as void. Also reported: the share that come out SURVIVES.
     Also reported: the share of rotated D at or above the real D (a
     one-sided permutation p).
  2. Positive injection (can it see a real effect?): on each of the first
     200 rotated draws of control 1 (no weekly cycle), multiply every MID
     station-day's PRCP by 1.10 and apply the verdict rule. That plants a
     10 % midweek excess on rain with no cycle. If fewer than 80 % of those
     come out SURVIVES, the study is VOID: it could not reliably have seen a
     cycle twice the size of the bar.
  3. Fewer than 40 of the 64 stations with at least 20 valid summers in the
     primary window: VOID.

Nothing here is a forecast or advice; it is a test of a published claim.
"""

import datetime as dt
import gzip
import random
from pathlib import Path

DATA = Path.home() / "data" / "weekly-rain"
URL = "https://www.ncei.noaa.gov/pub/data/ghcn/daily/by_station/{}.csv.gz"

STATIONS = """USW00003810 USW00003812 USW00003813 USW00003818 USW00003820
USW00003822 USW00003847 USW00003856 USW00003870 USW00003937 USW00003940
USW00012815 USW00012816 USW00012832 USW00012834 USW00012836 USW00012839
USW00012842 USW00012843 USW00012844 USW00012849 USW00012916 USW00013717
USW00013722 USW00013723 USW00013748 USW00013754 USW00013782 USW00013786
USW00013857 USW00013865 USW00013869 USW00013870 USW00013873 USW00013874
USW00013876 USW00013877 USW00013880 USW00013881 USW00013882 USW00013883
USW00013889 USW00013891 USW00013893 USW00013894 USW00013895 USW00013897
USW00013899 USW00013939 USW00013957 USW00013970 USW00013978 USW00053819
USW00053860 USW00053864 USW00053865 USW00053868 USW00093729 USW00093801
USW00093805 USW00093846 USW00093862 USW00093874 USW00093915""".split()

MID = {1, 2, 3}  # Tue Wed Thu (Monday = 0)
END = {5, 6, 0}  # Sat Sun Mon
PRIMARY = (1998, 2025)
ARMS = {"B_trmm": (1998, 2005), "C_before": (1960, 1989)}


def load(sid):
    """{date: tenths_mm} for valid JJA PRCP of one station."""
    out = {}
    with gzip.open(DATA / f"{sid}.csv.gz", "rt") as f:
        for line in f:
            p = line.rstrip("\n").split(",")
            if p[2] != "PRCP" or p[5] != "":
                continue
            d = dt.date(int(p[1][:4]), int(p[1][4:6]), int(p[1][6:8]))
            if 6 <= d.month <= 8:
                out[d] = int(p[3])
    return out


def summers(all_data, y0, y1):
    """{year: [(weekday, tenths_mm), ...]} over valid station-summers."""
    by_year, dropped = {}, 0
    for sid, rec in all_data.items():
        for y in range(y0, y1 + 1):
            days = [(d.weekday(), v) for d, v in rec.items() if d.year == y]
            if len(days) >= 85:
                by_year.setdefault(y, []).extend(days)
            else:
                dropped += 1
    return by_year, dropped


def sums(days, wet=False, shift=0, boost=1.0):
    m = e = nm = ne = 0.0
    for wd, v in days:
        wd = (wd + shift) % 7
        x = (1.0 if v >= 10 else 0.0) if wet else v
        if wd in MID:
            m += x * boost; nm += 1
        elif wd in END:
            e += x; ne += 1
    return m, nm, e, ne


def D(parts):
    m = sum(p[0] for p in parts); nm = sum(p[1] for p in parts)
    e = sum(p[2] for p in parts); ne = sum(p[3] for p in parts)
    return (m / nm) / (e / ne) - 1


def boot(year_parts, draws, seed):
    ys = sorted(year_parts)
    rng = random.Random(seed)
    vals = sorted(D([year_parts[rng.choice(ys)] for _ in ys]) for _ in range(draws))
    return vals[int(0.025 * draws)], vals[int(0.975 * draws) - 1]


def verdict(d, lo, hi):
    if lo > 0 and d >= 0.05:
        return "SURVIVES"
    if hi < 0.05:
        return "KILLED"
    return "UNRESOLVED"