"""The sealed rule for the Ice Saints.

Sealed before any row of March to July reached this server. Before the seal I read only
August to October rows of seven of the 35 stations, in design.py, to see how noisy a
5-day window is against its flanks (year SD about 2.0 C, SE about 0.2 C on 97 years).

Claim tested. Central European weather lore: a cold spell returns on the Ice Saints,
11 to 15 May (Mamertus, Pancras, Servatius, Boniface, Sophie), and planting is safe after.

Frame, fixed from the station list alone: every DWD daily KL station whose record runs
from 1 January 1925 or earlier to 31 December 2025 or later, at 800 m or lower, 35
stations, listed in FRAME. Years 1925 to 2025, element TMK (daily mean air temperature,
C). DWD's historical files, not homogenised; a contrast inside one month does not care
about a station move years apart.

1. For a window starting on date s in year y at one station:
   C(s, y) = 1/2 (mean T of s-10..s-1 + mean T of s+5..s+14) - mean T of s..s+4.
   Positive means the window is colder than the days either side. The station counts
   if the window has all 5 days and the two flanks together miss at most 2 of 20.
2. A year's value is the mean over stations that count, and the year counts if at least
   18 stations do. CBAR(s) is the mean over years that count.
3. Target: s = 11 May. C = CBAR(11 May).
4. Placebos: s every 5 days from 1 March to 31 October, leaving out any s whose span
   s-10..s+14 touches 1 to 25 May or 4 to 20 June.
5. SE = the larger of the bootstrap SE of the target's year values (10,000 draws, seed
   SEED) and the sample SD of the placebos' CBAR. Interval C +/- 1.96 SE.
6. Verdict, bar 1.0 C (a cold spell is at least a degree colder than the days around it):
   killed if the upper bound is below 1.0; survived if the lower bound is above 0 and
   C is at least 1.0; inconclusive otherwise.

Void if any of:
V1 2.8 x SE >= 1.0.
V2 February 1929, the great cold wave: the 1929 year value at s = 10 February must
   exceed 1.0 C (the sign and the date alignment).
V3 Over all stations and years mean July T minus mean January T must exceed 12 C.
V4 1.5 C planted (subtracted from all stations and years) in the window of each of ten
   placebos, every len//10-th from the first: at least 8 of 10 survived, each judged
   with the SE of step 5.
V5 At least 90 years count and at least 25 stations have data; every T in -40..45.

Printed only, cannot move the verdict: 1951 to 2025 and 1996 to 2025 alone; stations
north and south of 51 N; the target's rank among all 5-day starts of the placebos' list;
s = 9 June (the sheep cold); the share of station-years with a frost (TNK < 0) on or after
16 May, by 30-year period, at stations at or below 400 m.
"""

import datetime as dt

import numpy as np

FRAME = """00164 00314 00320 00445 00691 00853 00880 01443 01526 01580 01684 01691 01757 02444
02486 02578 02629 02928 03015 03126 03196 03631 03668 03946 03987 04024 04104 04445 04625
04748 05099 05349 05440 05797 05906""".split()
Y0, Y1 = 1925, 2025
BAR = 1.0
DRAWS = 10_000
SEED = 20261105
PLANT = 1.5
Z = 1.959964
MIN_ST = 18
ORIGIN = dt.date(Y0 - 1, 12, 1).toordinal()


def station_value(a, s_ord):
    """C(s, y) on one station's array a, NaN for missing, indexed by ordinal - ORIGIN."""
    i = s_ord - ORIGIN
    w = a[i : i + 5]
    f1, f2 = a[i - 10 : i], a[i + 5 : i + 15]
    if np.isnan(w).any() or np.isnan(f1).sum() + np.isnan(f2).sum() > 2:
        return np.nan
    if np.isnan(f1).all() or np.isnan(f2).all():
        return np.nan
    return 0.5 * (np.nanmean(f1) + np.nanmean(f2)) - w.mean()


def year_values(arrays, month, day, y0=Y0, y1=Y1, plant=0.0):
    out = {}
    for y in range(y0, y1 + 1):
        s = dt.date(y, month, day).toordinal()
        v = [station_value(a, s) for a in arrays]
        v = np.array(v, float)
        ok = ~np.isnan(v)
        if ok.sum() >= MIN_ST:
            out[y] = v[ok].mean() + plant
    return out


def cbar(arrays, month, day, **kw):
    return float(np.mean(list(year_values(arrays, month, day, **kw).values())))


def placebos():
    res = []
    d = dt.date(2025, 3, 1)
    end = dt.date(2025, 10, 31)
    bad = [(dt.date(2025, 5, 1), dt.date(2025, 5, 25)), (dt.date(2025, 6, 4), dt.date(2025, 6, 20))]
    while d <= end:
        a, b = d - dt.timedelta(10), d + dt.timedelta(14)
        if not any(a <= q and p <= b for p, q in bad):
            res.append((d.month, d.day))
        d += dt.timedelta(5)
    return res


def bootstrap_se(values):
    v = np.array(values, float)
    rng = np.random.default_rng(SEED)
    idx = rng.integers(0, len(v), (DRAWS, len(v)))
    return float(v[idx].mean(axis=1).std(ddof=1))


def verdict(c, se):
    lo, hi = c - Z * se, c + Z * se
    if hi < BAR:
        return "killed", lo, hi
    if lo > 0 and c >= BAR:
        return "survived", lo, hi
    return "inconclusive", lo, hi