"""Turn the raw hourly JSON into one row per (site, day), by rule.py's definitions. Same definitions as the original study, re-implemented here so this pack stands alone: an hour is WET at >= 0.1 mm; R = any wet hour 00:00..06:59; F = the hour 10:00..10:59 is dry; F_strict = every hour 08:00..10:59 dry; a day is dropped if any hour 00:00..10:59 is missing. Two block indices are built, per `blocking-ambiguity.md`: block (site, year-month) -- the verdict's blocking block96 (site, calendar month 1..12) -- the sensitivity and one more index used by the primary bootstrap: ym (year-month), shared across sites, so a drawn month brings the rows of every site at once. """ import glob import json import os import numpy as np OUT = os.path.expanduser("~/data/rain-before-seven-replication") UK_SITES = ["aberdeen", "inverness", "carlisle", "hull", "norwich", "swansea", "dublin", "galway"] PNW_SITES = ["seattle", "portland", "vancouver"] WET_MM = 0.1 START, END = "2000-01-01", "2025-12-31" def _site_hours(site): """{'YYYY-MM-DDTHH:00': mm} for one site, over every chunk file.""" hours = {} for path in sorted(glob.glob(os.path.join(OUT, f"{site}_*.json"))): with open(path) as f: d = json.load(f) for t, p in zip(d["hourly"]["time"], d["hourly"]["precipitation"]): hours[t] = p return hours def days(sites=None, wet_mm=WET_MM): """One record per (site, day) with the flags rule.py defines.""" sites = list(sites or UK_SITES) site_idx, blk, blk96, ym_s, date_s = [], [], [], [], [] R, F, Fs, dryh = [], [], [], [] dropped = 0 for si, site in enumerate(sites): hours = _site_hours(site) if not hours: raise SystemExit(f"no hours loaded for {site} -- run fetch.py first") for d in sorted({t[:10] for t in hours}): if not (START <= d <= END): continue vals = [hours.get(f"{d}T{h:02d}:00") for h in range(15)] if any(v is None for v in vals[:11]): dropped += 1 continue dryh.append([np.nan if v is None else float(v < wet_mm) for v in vals]) vals = vals[:11] site_idx.append(si) blk.append(f"{site}|{d[:7]}") blk96.append(f"{site}|{d[5:7]}") ym_s.append(d[:7]) date_s.append(d) R.append(any(v >= wet_mm for v in vals[0:7])) F.append(vals[10] < wet_mm) Fs.append(all(v < wet_mm for v in vals[8:11])) def _index(keys): u = sorted(set(keys)) m = {k: i for i, k in enumerate(u)} return np.array([m[k] for k in keys]), len(u) b, nb = _index(blk) b96, nb96 = _index(blk96) y, ny = _index(ym_s) return { "sites": sites, "site": np.array(site_idx), "block": b, "n_blocks": nb, "block96": b96, "n_blocks96": nb96, "ym": y, "n_ym": ny, "date": np.array(date_s), "R": np.array(R), "F": np.array(F), "Fstrict": np.array(Fs), "dryh": np.array(dryh), "dropped": dropped, } def d_mh(block, n_blocks, R, F): """Mantel-Haenszel risk difference in percentage points, rule.py's formula. D_MH = 100 * sum_b w_b (s1b/n1b - s0b/n0b) / sum_b w_b, w_b = n1b n0b / (n1b + n0b) Blocks with n1b = 0 or n0b = 0 are dropped. """ Fv = np.asarray(F, dtype=float) R = np.asarray(R, dtype=bool) n1 = np.bincount(block[R], minlength=n_blocks).astype(float) n0 = np.bincount(block[~R], minlength=n_blocks).astype(float) s1 = np.bincount(block[R], weights=Fv[R], minlength=n_blocks) s0 = np.bincount(block[~R], weights=Fv[~R], minlength=n_blocks) ok = (n1 > 0) & (n0 > 0) if not ok.any(): return float("nan") w = n1[ok] * n0[ok] / (n1[ok] + n0[ok]) diff = s1[ok] / n1[ok] - s0[ok] / n0[ok] return 100.0 * float((w * diff).sum() / w.sum()) def _members(idx, n): order = np.argsort(idx, kind="stable") s = np.searchsorted(idx[order], np.arange(n), side="left") e = np.searchsorted(idx[order], np.arange(n), side="right") return [order[a:b] for a, b in zip(s, e)] def bootstrap(d, R, F, unit, reps=10000, seed=20260922, block_key="block"): """Fast path, algebraically identical to bootstrap_slow below. Every (site, year-month) block sits inside exactly one resampling unit -- inside its own month when the unit is the month across all sites, inside itself when the unit is the block. So a drawn unit always brings whole blocks, and each block's (n1b, n0b, s1b, s0b) are the same numbers whoever draws it. D_MH over a resample is therefore 100 * sum_over_drawn(w_b * diff_b) / sum_over_drawn(w_b) with w_b and diff_b precomputed once. A unit drawn twice contributes its terms twice, which is exactly what "keeps a distinct block identity" means for this estimator. Per rep this is a sum over 312 units, not a rebuild of 76,000 rows. `test_bootstrap.py` checks the two paths agree. """ rng = np.random.default_rng(seed) nu = {"ym": d["n_ym"], "block": d["n_blocks"], "block96": d["n_blocks96"]}[unit] bkey, nb = d[block_key], {"block": d["n_blocks"], "block96": d["n_blocks96"]}[block_key] Fv = np.asarray(F, dtype=float) Rb = np.asarray(R, dtype=bool) n1 = np.bincount(bkey[Rb], minlength=nb).astype(float) n0 = np.bincount(bkey[~Rb], minlength=nb).astype(float) s1 = np.bincount(bkey[Rb], weights=Fv[Rb], minlength=nb) s0 = np.bincount(bkey[~Rb], weights=Fv[~Rb], minlength=nb) ok = (n1 > 0) & (n0 > 0) w = np.where(ok, n1 * n0 / np.maximum(n1 + n0, 1), 0.0) wd = np.where(ok, w * (s1 / np.maximum(n1, 1) - s0 / np.maximum(n0, 1)), 0.0) # map each block to its unit, then total w and w*diff per unit if unit == "ym": b2u = np.zeros(nb, dtype=np.int64) b2u[bkey] = d["ym"] else: b2u = np.arange(nb) SW = np.bincount(b2u, weights=w, minlength=nu) SWD = np.bincount(b2u, weights=wd, minlength=nu) chunks = [] done = 0 while done < reps: k = min(1000, reps - done) pick = rng.integers(0, nu, size=(k, nu)) num = SWD[pick].sum(axis=1) den = SW[pick].sum(axis=1) chunks.append(100.0 * num / np.where(den == 0, np.nan, den)) done += k out = np.concatenate(chunks) out = out[~np.isnan(out)] return float(np.percentile(out, 2.5)), float(np.percentile(out, 97.5)) def bootstrap_slow(d, R, F, unit, reps=10000, seed=20260922, block_key="block"): """95 % percentile interval for D_MH, resampling `unit` with replacement. unit="ym" PRIMARY: whole calendar months, all sites at once. unit="block" SECONDARY: the original's (site, month) blocks. A unit drawn twice keeps a distinct block identity inside D_MH, so it is counted twice, as rule.py requires. When the resampling unit is the month across sites, a drawn month still splits into its per-site blocks -- the draw is the month, the stratum is the site-month. """ rng = np.random.default_rng(seed) uidx = d[unit] nu = d["n_ym"] if unit == "ym" else d[{"block": "n_blocks", "block96": "n_blocks96"}[unit]] mem = _members(uidx, nu) bkey = d[block_key] # per unit: its rows, and its blocks relabelled densely 0..m_k-1, so that a # unit drawn twice keeps two distinct block identities and the stratum count # stays proportional to the draw rather than to the global block id range local, msize = [], [] for rows in mem: _, inv = np.unique(bkey[rows], return_inverse=True) local.append(inv.astype(np.int64)) msize.append(int(inv.max()) + 1 if len(inv) else 0) msize = np.array(msize) out = np.empty(reps) for i in range(reps): pick = rng.integers(0, nu, nu) idx = np.concatenate([mem[j] for j in pick]) offs = np.concatenate([[0], np.cumsum(msize[pick])[:-1]]) rb = np.concatenate([local[k] + o for k, o in zip(pick, offs)]) out[i] = d_mh(rb, int(msize[pick].sum()), R[idx], F[idx]) out = out[~np.isnan(out)] return float(np.percentile(out, 2.5)), float(np.percentile(out, 97.5))