"""Daily EHD, exceedance, weekly aggregates and heat-wave episode accounting. All definitions here are LOCKED (see docs/decision_log.md). Key distinctions that fail quietly if confused, and are therefore enforced and unit-tested: * HW2_WEEK (>=2 EHDs anywhere in the week, no consecutive-day condition) is a DIFFERENT concept from a heat-wave EPISODE (>=2 CONSECUTIVE EHDs). They are kept in separate functions and never conflated. * Episodes are detected on the CONTINUOUS daily sequence BEFORE any weekly aggregation. An episode crossing a week boundary is ONE event. * Cross-week accounting never duplicates full duration or full severity across overlapping weeks; those belong to the onset week only, and never to a right-censored episode. """ import numpy as np # Lead-day (LT) ranges, inclusive. W7 is a 4-day PARTIAL week, verified # separately and never pooled with W1-W6. DEFAULT_WEEKS = [ ("W1", 0, 6), ("W2", 7, 13), ("W3", 14, 20), ("W4", 21, 27), ("W5", 28, 34), ("W6", 35, 41), ("W7", 42, 45), # partial (4 days) ] def daily_ehd(tx, thr): """EHD(s,t) = 1 if Tx(s,t) >= EHD_THRESHOLD(s,d) else 0.""" return (np.asarray(tx, dtype="float64") >= np.asarray(thr, dtype="float64")).astype(np.int8) def exceedance(tx, thr): """exceedance(s,t) = max(Tx - EHD_THRESHOLD, 0).""" return np.maximum(np.asarray(tx, dtype="float64") - np.asarray(thr, dtype="float64"), 0.0) def _week_of_lead(lead, weeks): for name, lo, hi in weeks: if lo <= lead <= hi: return name return None def weekly_aggregate(lead_days, ehd, exc, weeks=DEFAULT_WEEKS): """Per-week EHD_COUNT, EHE_WEEK, HW2_WEEK, SEV_WEEK for one 1-D sequence. Parameters ---------- lead_days : 1-D int array of lead days (LT) for each daily step. ehd, exc : 1-D arrays aligned with lead_days. Returns a dict keyed by week name with sub-keys EHD_COUNT, EHE_WEEK, HW2_WEEK, SEV_WEEK. Buffer days (lead outside any week range) are ignored. """ lead_days = np.asarray(lead_days, dtype=int) ehd = np.asarray(ehd) exc = np.asarray(exc, dtype="float64") out = {} for name, lo, hi in weeks: m = (lead_days >= lo) & (lead_days <= hi) count = int(ehd[m].sum()) out[name] = { "EHD_COUNT": count, "EHE_WEEK": int(count >= 1), "HW2_WEEK": int(count >= 2), "SEV_WEEK": float(exc[m][ehd[m] == 1].sum()) if m.any() else 0.0, } return out def detect_episodes(ehd, exc, min_len=2): """Detect heat-wave episodes as maximal runs of >= min_len consecutive EHDs. Operates on a CONTINUOUS daily sequence (including any buffer). Returns a list of dicts with start/end indices (into the given arrays), duration, mean/max intensity and summed severity (degC and degC-day). """ ehd = np.asarray(ehd).astype(int) exc = np.asarray(exc, dtype="float64") n = len(ehd) episodes = [] i = 0 while i < n: if ehd[i] == 1: j = i while j + 1 < n and ehd[j + 1] == 1: j += 1 if (j - i + 1) >= min_len: e = exc[i:j + 1] episodes.append({ "start": i, "end": j, "duration": int(j - i + 1), "i_mean": float(np.mean(e)), "i_max": float(np.max(e)), "sev": float(np.sum(e)), }) i = j + 1 else: i += 1 return episodes def analyze_forecast_sequence(lead_days, ehd, exc, weeks=DEFAULT_WEEKS, horizon=(0, 45)): """Full episode + cross-week accounting for one continuous daily sequence. ``lead_days`` may extend below 0 (pre-LT0 buffer) and above the horizon end (post-LT45 buffer); those buffer days are used for episode detection and continuity ONLY. The official horizon stays ``horizon``. Censoring: an episode touching the earliest available day is flagged ``left_censored``; one touching the latest available day is ``right_censored``. A left-censored episode contributes to HW_ACTIVE and HW_DAYS but is NOT an onset. A right-censored episode contributes to onset and active-day statistics but never to full-duration statistics. Returns (episodes, weekly) where ``weekly`` has per-week HW_ONSET, HW_ONSET_COUNT, HW_ACTIVE, HW_DAYS, HW_SEV_ACTIVE. """ lead_days = np.asarray(lead_days, dtype=int) exc = np.asarray(exc, dtype="float64") lo_h, hi_h = horizon n = len(lead_days) raw = detect_episodes(ehd, exc) weekly = {name: {"HW_ONSET": 0, "HW_ONSET_COUNT": 0, "HW_ACTIVE": 0, "HW_DAYS": 0, "HW_SEV_ACTIVE": 0.0} for name, _, _ in weeks} episodes = [] for e in raw: ld_start = int(lead_days[e["start"]]) ld_end = int(lead_days[e["end"]]) left_censored = e["start"] == 0 right_censored = e["end"] == n - 1 # onset: counted once, in the week it starts; only if the start is # observed (not left-censored) and lies inside the official horizon. onset_week = None if not left_censored and lo_h <= ld_start <= hi_h: onset_week = _week_of_lead(ld_start, weeks) # per-day active accounting, only for days inside the official horizon active_weeks = set() for k in range(e["start"], e["end"] + 1): ld = int(lead_days[k]) if not (lo_h <= ld <= hi_h): continue wk = _week_of_lead(ld, weeks) if wk is None: continue active_weeks.add(wk) weekly[wk]["HW_DAYS"] += 1 # every day inside a detected episode is an EHD by construction weekly[wk]["HW_SEV_ACTIVE"] += float(exc[k]) for wk in active_weeks: weekly[wk]["HW_ACTIVE"] = 1 if onset_week is not None: weekly[onset_week]["HW_ONSET"] = 1 weekly[onset_week]["HW_ONSET_COUNT"] += 1 episodes.append({ **e, "lead_start": ld_start, "lead_end": ld_end, "onset_week": onset_week, "left_censored": bool(left_censored), "right_censored": bool(right_censored), # full duration/severity are valid onset-based statistics only for # complete (not right-censored) episodes with an onset in-horizon. "full_stats_valid": bool(onset_week is not None and not right_censored), }) return episodes, weekly