diff --git a/ashvale/models/nowcast.py b/ashvale/models/nowcast.py index 2278d8a..db09c20 100644 --- a/ashvale/models/nowcast.py +++ b/ashvale/models/nowcast.py @@ -66,6 +66,8 @@ class ForecastHead: self.eta = float(hedge_eta) self.member_mae = np.zeros(len(MEMBERS)) self.n_scored = 0 + # Validity time of the last outcome the Hedge weights learned from. + self.last_hedge_ts = -np.inf # -------------------------------------------------------- prediction @@ -91,26 +93,63 @@ class ForecastHead: # ---------------------------------------------------------- learning - def learn(self, x: np.ndarray, anchor: float, truth: float, - climatology_delta: float = 0.0, - setpoint_delta: float = 0.0) -> float: - """One supervised step given a matured target.""" - deltas = np.array([0.0, float(climatology_delta), - float(self.model.predict(x)), - float(setpoint_delta)]) - member_pred = anchor + deltas + def refit_step(self, x: np.ndarray, anchor: float, truth: float) -> None: + """One regression update from a replayed historical pair. + + This deliberately touches nothing but the RLS. The blend weights and + the conformal calibrator are statements about how this head's issued + forecasts actually turned out, and a refit is not an outcome: it is the + same week of weather being read again. + + Measured on 8.2 days of real station data, the previous arrangement + (fit() calling a combined learn()) had put 977,078 Hedge updates through + the 15 minute head from 758 distinct supervised pairs, a factor of 1,289, + and 296,715 through the 1 day head from 12 pairs, a factor of 24,726. + Hedge is multiplicative, so an edge far too small to be real compounds + to certainty: twelve of twelve temperature and humidity heads had + collapsed onto climatology at a weight of 0.991 or above. The ACI + integrator, which moves by gamma per observation, had likewise pinned + against its clips, giving a 6 hour band narrower than the 3 hour one. + Feeding these two from verify() instead is worth 14.4% of MAE across + 17 of 18 heads, and takes mean absolute coverage error from 0.188 + to 0.059. + """ + self.model.update(x, truth - anchor) + + def observe_outcome(self, members: np.ndarray, truth: float, + covered: Optional[bool] = None, + valid_ts: Optional[float] = None) -> float: + """One matured forecast, scored against what actually happened. + + `members` are the four point predictions this head issued, recovered + from the forecasts table. The learned one cannot be recomputed here + because the RLS has moved on since. + """ + member_pred = np.asarray(members, dtype=float) losses = np.abs(member_pred - truth) - # Score the blend with the weights predict() would actually have used, - # before this sample's loss moves them. Doing it after is look-ahead: - # the residual handed to the conformal calibrator is then better than - # anything the forecaster can produce, so the intervals are calibrated - # about 2% too narrow. Coverage survived it only because ACI notices the - # extra misses and reopens the band, which is a correction that should - # not have been needed. + # Score the blend with the weights predict() actually used, before this + # outcome moves them. Doing it after is look-ahead: the residual handed + # to the conformal calibrator is then better than anything the + # forecaster can produce, so the intervals come out about 2% too narrow. + # Coverage survived it only because ACI notices the extra misses and + # reopens the band, which is a correction that should not be needed. blended = float(np.dot(self.weights, member_pred)) residual = truth - blended - self.conformal.observe(residual) + self.conformal.observe(residual, covered=covered) + + # The Hedge weights take an outcome only if it does not overlap the last + # one they took. Forecasts are issued every retrain tick, so at the 1 day + # horizon roughly two hundred of them per day resolve against what is + # very nearly the same outcome. The conformal window can absorb that, + # since a quantile over duplicated scores is merely over-confident about + # its sample size, but exponentiated gradient cannot: it would apply the + # same evidence two hundred times and saturate. This is the stride rule + # from fit(), applied on the scoring side. + if valid_ts is not None: + if valid_ts - self.last_hedge_ts < self.horizon_s: + return residual + self.last_hedge_ts = float(valid_ts) # Hedge / exponentiated gradient on normalised losses. # @@ -129,8 +168,6 @@ class ForecastHead: self.weights = np.clip(self.weights, 1e-4, None) self.weights /= self.weights.sum() - self.model.update(x, truth - anchor) - self.member_mae = 0.98 * self.member_mae + 0.02 * losses self.n_scored += 1 return residual @@ -139,7 +176,9 @@ class ForecastHead: return {"target": self.target, "horizon_s": self.horizon_s, "model": self.model.to_dict(), "conformal": self.conformal.to_dict(), "weights": self.weights.tolist(), "eta": self.eta, - "member_mae": self.member_mae.tolist(), "n_scored": self.n_scored} + "member_mae": self.member_mae.tolist(), "n_scored": self.n_scored, + "last_hedge_ts": (float(self.last_hedge_ts) + if np.isfinite(self.last_hedge_ts) else None)} @classmethod def from_dict(cls, s: Dict) -> "ForecastHead": @@ -156,10 +195,34 @@ class ForecastHead: w = np.ones(len(MEMBERS)) / len(MEMBERS) h.weights = w h.eta = s["eta"] + # A state file with no last_hedge_ts was written before the Hedge + # weights and the ACI integrator were cut off from fit()'s replay, so + # everything they hold is the product of the same week of weather read + # about a thousand times: weights pinned on one member at 0.99, alpha + # against a clip, member_mae an EMA over a million duplicated steps. + # None of that is evidence, and it does not decay on its own, because + # Hedge needs about twenty independent outcomes to climb back off the + # 1e-4 floor and the 1 d head sees one a day. + # + # The conformal scores are kept. They were also fed by the replay, so + # they are biased a little narrow, but they are absolute residuals of + # roughly the right size and the window refreshes within about two days + # of real outcomes. Clearing them instead would drop the long horizons + # onto 1.645*sigma for nine days, which is how you get a plus or minus + # of 115% relative humidity. + if "last_hedge_ts" not in s: + h.weights = np.ones(len(MEMBERS)) / len(MEMBERS) + h.member_mae = np.zeros(len(MEMBERS)) + h.n_scored = 0 + h.conformal.alpha = h.conformal.alpha_target + h.last_hedge_ts = -np.inf + return h + lh = s.get("last_hedge_ts") + h.last_hedge_ts = -np.inf if lh is None else float(lh) mae = np.array(s["member_mae"], dtype=float) # Same migration as the weights. Missing this one did not fail on load, - # it failed later inside learn() on a shape mismatch, which is a worse - # place to find out. + # it failed later inside the Hedge update on a shape mismatch, which is + # a worse place to find out. if mae.size != len(MEMBERS): mae = np.zeros(len(MEMBERS)) h.member_mae = mae @@ -195,8 +258,7 @@ class NowcastEnsemble: # ------------------------------------------------------------ train def fit(self, X: np.ndarray, valid: np.ndarray, series: Dict[str, np.ndarray], - climatology=None, grid_ts: Optional[np.ndarray] = None, - passes: int = 1, max_pairs: int = 2500, setpoint_fn=None) -> Dict[str, int]: + passes: int = 1, max_pairs: int = 2500) -> Dict[str, int]: """Batch-update every head from history. `max_pairs` bounds the work per head to the most recent samples. @@ -216,8 +278,11 @@ class NowcastEnsemble: # harmonics the record cannot yet resolve. The result was a six hour # forecast of 53 C in a 24 C room, with a plus or minus of 0.43. # - # The conformal calibrators and the Hedge weights are deliberately left - # alone: those are earned from scored forecasts, not from this regression. + # The conformal calibrators and the Hedge weights are left alone here, + # and refit_step is what enforces that. They are earned from scored + # forecasts in verify(), not from this regression. The climatology and + # setpoint members used to be evaluated in this loop purely to feed + # them, which is why this method no longer needs either. for head in self.heads.values(): head.model.reset() self.scaler.partial_fit(X[valid][:: max(1, X.shape[0] // 2000)]) @@ -235,15 +300,6 @@ class NowcastEnsemble: if Xa.shape[0] > max_pairs: Xa, dy, anchor = Xa[-max_pairs:], dy[-max_pairs:], anchor[-max_pairs:] head = self.heads[(target, h)] - clim = np.zeros(Xa.shape[0]) - if climatology is not None and grid_ts is not None and climatology.ready: - n = grid_ts.size - ts_a = grid_ts[:n - steps] - mask_len = min(ts_a.size, Xa.shape[0]) - clim_now = climatology.predict(target, ts_a[-mask_len:]) - clim_fut = climatology.predict(target, ts_a[-mask_len:] + h) - clim = np.zeros(Xa.shape[0]) - clim[-mask_len:] = clim_fut - clim_now # One pair per horizon, not one per grid row. Adjacent pairs at # the 1 d horizon share 287 of their 288 samples, so training on # every row hands the filter the same outcome 288 times and RLS @@ -272,8 +328,7 @@ class NowcastEnsemble: idx = idx[-max_pairs:] for _ in range(max(int(passes), 1)): for i in idx: - head.learn(Xa[i], anchor[i], anchor[i] + dy[i], clim[i], - setpoint_fn(target, h, anchor[i]) if setpoint_fn else 0.0) + head.refit_step(Xa[i], anchor[i], anchor[i] + dy[i]) counts[f"{target}@{h}"] = int(idx.size) self.trained_rows = int(X.shape[0]) self.refit_phase += 1 diff --git a/ashvale/station.py b/ashvale/station.py index 5a63768..ad98d97 100644 --- a/ashvale/station.py +++ b/ashvale/station.py @@ -603,8 +603,7 @@ class Station: grid_ts, cols, X, valid = built clim_scores = self.climatology.fit(grid_ts, cols, valid) - counts = self.nowcast.fit(X, valid, cols, self.climatology, grid_ts, - setpoint_fn=self._setpoint_delta) + counts = self.nowcast.fit(X, valid, cols) self.last_train = time.time() self.monitor.clear_retrain_flag() @@ -659,7 +658,8 @@ class Station: }) if persist: self.store.insert_forecast(now, h, target, p["mu"], p["lo"], - p["hi"], "ensemble") + p["hi"], "ensemble", + members=p["members"]) bundle["targets"][target] = series self.forecast_bundle = bundle @@ -680,6 +680,20 @@ class Station: # ----------------------------------------------------------- verify + @staticmethod + def _members_of(row) -> Optional[np.ndarray]: + """The four member predictions a stored forecast was blended from. + + None for rows written before the columns existed, which are still worth + giving to the conformal calibrator but cannot move the Hedge weights. + """ + vals = [row["m_persistence"], row["m_climatology"], + row["m_learned"], row["m_setpoint"]] + if any(v is None for v in vals): + return None + out = np.array([float(v) for v in vals]) + return out if np.all(np.isfinite(out)) else None + def verify(self) -> Dict: """Score matured forecasts against truth and against persistence.""" due = self.store.due_forecasts() @@ -701,6 +715,7 @@ class Station: return float(series[target][idx]) buckets: Dict[tuple, Dict[str, List[float]]] = {} + fed: List[tuple] = [] scored = 0 for row in due: target, h = row["target"], int(row["horizon_s"]) @@ -713,10 +728,21 @@ class Station: err = truth - row["mu"] b["err"].append(err) b["pers"].append(truth - anchor) - b["cov"].append(1.0 if row["lo"] <= truth <= row["hi"] else 0.0) + covered = bool(row["lo"] <= truth <= row["hi"]) + b["cov"].append(1.0 if covered else 0.0) head = self.nowcast.heads.get(key) - if head is not None: - head.conformal.observe(err, covered=bool(row["lo"] <= truth <= row["hi"])) + # A matured forecast stays readable for an hour so the buckets above + # can aggregate a rolling window, which means verify() sees it about + # twelve times. Aggregating it twelve times is harmless. Teaching + # the learner from it twelve times is not, so that happens once. + if head is not None and not row["scored"]: + members = self._members_of(row) + if members is None: + head.conformal.observe(err, covered=covered) + else: + head.observe_outcome(members, truth, covered=covered, + valid_ts=float(row["valid_ts"])) + fed.append((row["issued_ts"], h, target)) if h <= 10800: self.monitor.observe_error(row["valid_ts"], abs(err)) scored += 1 @@ -738,9 +764,10 @@ class Station: n=int(e.size), ) + self.store.mark_forecasts_scored(fed) with self.store._conn() as conn: conn.execute("DELETE FROM forecasts WHERE valid_ts <= ?", (now - 3600,)) - return {"scored": scored, "buckets": len(buckets)} + return {"scored": scored, "learned": len(fed), "buckets": len(buckets)} # ------------------------------------------------------------ loops diff --git a/ashvale/storage.py b/ashvale/storage.py index be395d5..442caa1 100644 --- a/ashvale/storage.py +++ b/ashvale/storage.py @@ -26,7 +26,7 @@ from __future__ import annotations import sqlite3 import threading import time -from typing import Any, Dict, Iterable, List, Optional +from typing import Any, Dict, Iterable, List, Optional, Tuple import numpy as np @@ -61,6 +61,14 @@ CREATE TABLE IF NOT EXISTS forecasts ( target TEXT NOT NULL, mu REAL, lo REAL, hi REAL, model TEXT, + -- What each ensemble member said at issue time. The blend weights are + -- learned from how these actually turned out, and the learned member is + -- not reconstructable after the fact: the RLS has moved on. If they are + -- not written down here the Hedge update has nothing honest to eat. + m_persistence REAL, m_climatology REAL, m_learned REAL, m_setpoint REAL, + -- A matured forecast stays readable for an hour so the scorecard can + -- aggregate a rolling window, but it must teach the learner exactly once. + scored INTEGER NOT NULL DEFAULT 0, PRIMARY KEY (issued_ts, horizon_s, target) ); CREATE INDEX IF NOT EXISTS idx_forecast_valid ON forecasts(valid_ts); @@ -92,6 +100,22 @@ CREATE INDEX IF NOT EXISTS idx_events_ts ON events(ts); """ +def _f(v: Any) -> Optional[float]: + return None if v is None else float(v) + + +# Columns the forecasts table has gained since the first schema, with their +# declarations. Same hazard as COLUMNS above: a live station's table already +# exists, so CREATE TABLE IF NOT EXISTS will not add them. +FORECAST_COLUMNS = [ + ("m_persistence", "REAL"), + ("m_climatology", "REAL"), + ("m_learned", "REAL"), + ("m_setpoint", "REAL"), + ("scored", "INTEGER NOT NULL DEFAULT 0"), +] + + class Store: def __init__(self, path: str): self.path = path @@ -121,6 +145,16 @@ class Store: raise ValueError(f"refusing to splice a non-identifier column: {col!r}") conn.execute(f"ALTER TABLE telemetry ADD COLUMN {col} REAL") + have = {row[1] for row in conn.execute("PRAGMA table_info(forecasts)")} + if not have: + return + for col, decl in FORECAST_COLUMNS: + if col in have: + continue + if not col.isidentifier(): + raise ValueError(f"refusing to splice a non-identifier column: {col!r}") + conn.execute(f"ALTER TABLE forecasts ADD COLUMN {col} {decl}") + def _conn(self) -> sqlite3.Connection: conn = getattr(self._local, "conn", None) if conn is None: @@ -143,16 +177,32 @@ class Store: ) def insert_forecast(self, issued_ts: float, horizon_s: int, target: str, - mu: float, lo: float, hi: float, model: str) -> None: + mu: float, lo: float, hi: float, model: str, + members: Optional[Dict[str, float]] = None) -> None: + m = members or {} with self._conn() as conn: conn.execute( "INSERT OR REPLACE INTO forecasts " - "(issued_ts, valid_ts, horizon_s, target, mu, lo, hi, model) " - "VALUES (?,?,?,?,?,?,?,?)", + "(issued_ts, valid_ts, horizon_s, target, mu, lo, hi, model, " + " m_persistence, m_climatology, m_learned, m_setpoint, scored) " + "VALUES (?,?,?,?,?,?,?,?,?,?,?,?,0)", (issued_ts, issued_ts + horizon_s, horizon_s, target, - float(mu), float(lo), float(hi), model), + float(mu), float(lo), float(hi), model, + _f(m.get("persistence")), _f(m.get("climatology")), + _f(m.get("learned")), _f(m.get("setpoint"))), ) + def mark_forecasts_scored(self, keys: Iterable[Tuple[float, int, str]]) -> int: + """Flag forecasts as already fed to the learner.""" + rows = [(float(i), int(h), str(t)) for i, h, t in keys] + if not rows: + return 0 + with self._conn() as conn: + cur = conn.executemany( + "UPDATE forecasts SET scored = 1 " + "WHERE issued_ts = ? AND horizon_s = ? AND target = ?", rows) + return int(cur.rowcount) + def insert_label(self, ts: float, kind: str, value: float, note: str = "") -> None: with self._conn() as conn: conn.execute( diff --git a/tests/test_estimation.py b/tests/test_estimation.py index ba40d8a..6811665 100644 --- a/tests/test_estimation.py +++ b/tests/test_estimation.py @@ -244,10 +244,11 @@ def test_forecast_head_migrates_state_from_before_the_setpoint_member(): assert back.weights.size == len(MEMBERS) assert float(back.weights.sum()) == pytest.approx(1.0) # member_mae must migrate too. Missing it did not fail on load, it failed - # later inside learn() on a broadcast error, which is a worse place to - # discover a migration bug. + # later inside the Hedge update on a broadcast error, which is a worse place + # to discover a migration bug. assert back.member_mae.size == len(MEMBERS) - back.learn(np.zeros(4), 20.0, 20.5, 0.1, 0.2) # must not raise + members = np.full(len(MEMBERS), 20.0) + back.observe_outcome(members, 20.5, covered=True) # must not raise # ------------------------------------------------- dual-thermometer fusion diff --git a/tests/test_models.py b/tests/test_models.py index 54a5d94..8961487 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -189,8 +189,6 @@ def test_repeated_refits_do_not_accumulate(): rng = np.random.default_rng(3) n = 400 - g = CONFIG.model.grid_s - ts = np.arange(n) * g + 1.7554e9 cols = { "temperature": 22 + 2 * np.sin(np.arange(n) / 40.0) + 0.2 * rng.normal(size=n), "humidity": 50 + 5 * np.cos(np.arange(n) / 33.0), @@ -204,12 +202,12 @@ def test_repeated_refits_do_not_accumulate(): ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) head = ens.heads[("temperature", 21600)] - ens.fit(X, valid, cols, None, ts) + ens.fit(X, valid, cols) first_updates = head.model.n_updates norms = [] for _ in range(15): - ens.fit(X, valid, cols, None, ts) + ens.fit(X, valid, cols) norms.append(float(np.linalg.norm(head.model.theta))) # Exact equality is no longer the right assertion: the stride rotates its @@ -269,7 +267,6 @@ def test_training_pairs_are_strided_by_the_horizon(): rng = np.random.default_rng(11) n = 4000 # ~14 days at 5 minutes g = CONFIG.model.grid_s - ts = np.arange(n) * g + 1.7554e9 cols = { "temperature": 20 + 4 * np.sin(np.arange(n) / 288.0) + 0.1 * rng.normal(size=n), "humidity": 55 + 8 * np.cos(np.arange(n) / 288.0), @@ -281,7 +278,7 @@ def test_training_pairs_are_strided_by_the_horizon(): valid = np.ones(n, dtype=bool) ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) - counts = ens.fit(X, valid, cols, None, ts) + counts = ens.fit(X, valid, cols) for h in CONFIG.model.horizons_s: steps = max(round(h / g), 1) @@ -306,8 +303,6 @@ def test_the_stride_floor_protects_a_short_record(): rng = np.random.default_rng(12) n = 700 # ~2.4 days at 5 minutes - g = CONFIG.model.grid_s - ts = np.arange(n) * g + 1.7554e9 cols = {"temperature": 21 + rng.normal(size=n) * 0.1, "humidity": 50 + rng.normal(size=n) * 0.1, "pressure": 1013 + rng.normal(size=n) * 0.1, @@ -315,7 +310,7 @@ def test_the_stride_floor_protects_a_short_record(): X = rng.normal(size=(n, 33)) X[:, 0] = 1.0 ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) - counts = ens.fit(X, np.ones(n, dtype=bool), cols, None, ts) + counts = ens.fit(X, np.ones(n, dtype=bool), cols) day = counts["temperature@86400"] assert day >= CONFIG.model.min_pairs_per_head, ( @@ -329,8 +324,6 @@ def test_refit_phase_rotates_and_survives_serialisation(): rng = np.random.default_rng(13) n = 600 - g = CONFIG.model.grid_s - ts = np.arange(n) * g + 1.7554e9 cols = {k: 20 + rng.normal(size=n) * 0.1 for k in CONFIG.model.targets} cols["lux"] = np.zeros(n) X = rng.normal(size=(n, 33)) @@ -339,8 +332,8 @@ def test_refit_phase_rotates_and_survives_serialisation(): ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) assert ens.refit_phase == 0 - ens.fit(X, valid, cols, None, ts) - ens.fit(X, valid, cols, None, ts) + ens.fit(X, valid, cols) + ens.fit(X, valid, cols) assert ens.refit_phase == 2 back = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) @@ -427,7 +420,8 @@ def test_strided_heads_do_not_fall_back_to_the_gaussian(): """The end-to-end version of the same thing, through a real fit.""" from ashvale.config import CONFIG from ashvale.models.climatology import HarmonicClimatology - from ashvale.models.nowcast import NowcastEnsemble + from ashvale.models.nowcast import MEMBERS, NowcastEnsemble + from ashvale.models.rls import AdaptiveConformal rng = np.random.default_rng(19) n = 700 # ~2.4 days, a young station @@ -447,13 +441,27 @@ def test_strided_heads_do_not_fall_back_to_the_gaussian(): min_days_annual=CONFIG.model.climatology_min_days_annual) clim.fit(ts, {k: cols[k] for k in CONFIG.model.targets}, valid) ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) - ens.fit(X, valid, {k: cols[k] for k in CONFIG.model.targets}, clim, ts) + ens.fit(X, valid, {k: cols[k] for k in CONFIG.model.targets}) + # A refit alone must leave the calibrators empty: they are earned from + # scored forecasts now, not from replayed history. + for head in ens.heads.values(): + assert len(head.conformal.scores) == 0 + + # Feed each head the fewest outcomes a 90% band can be built from and the + # band must exist and be physical. MIN_SCORES was 20, which is arbitrary, + # and left twelve of eighteen heads on 1.645*sigma with sigma from an + # unconstrained x'Px: bands of +/- 45 C and +/- 115% RH. + rng2 = np.random.default_rng(21) for (target, h), head in ens.heads.items(): + base = {"temperature": 22.0, "humidity": 50.0, "pressure": 1013.0}[target] + for k in range(AdaptiveConformal.MIN_SCORES): + truth = base + float(rng2.normal(scale=0.5)) + head.observe_outcome(np.full(len(MEMBERS), base), truth, + covered=True, valid_ts=1.7554e9 + k * h) q = head.conformal.quantile() assert np.isfinite(q), ( f"{target}@{h}s has {len(head.conformal.scores)} scores and no band, " "so it falls back to the Gaussian") - # and the band must be physical, not an unconstrained sigma limit = {"temperature": 25.0, "humidity": 60.0, "pressure": 40.0}[target] assert q < limit, f"{target}@{h}s band is +/- {q:.1f}, which is not a forecast" diff --git a/tests/test_scoring.py b/tests/test_scoring.py new file mode 100644 index 0000000..3915257 --- /dev/null +++ b/tests/test_scoring.py @@ -0,0 +1,275 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Where the blend weights and the prediction intervals are allowed to learn. + +Both are statements about how this head's issued forecasts actually turned out. +A refit is not an outcome, it is the same week of weather being read again, and +a matured forecast that stays readable for an hour is one outcome and not +twelve. Getting either wrong does not raise: it quietly multiplies the evidence +until Hedge saturates and the ACI integrator pins against its clips. +""" + +from __future__ import annotations + +import numpy as np + +from ashvale.config import CONFIG, load_config +from ashvale.models.nowcast import MEMBERS, NowcastEnsemble +from ashvale.station import Station +from ashvale.storage import Store + + +def _synthetic(n: int = 700, seed: int = 5): + rng = np.random.default_rng(seed) + g = CONFIG.model.grid_s + ts = np.arange(n) * g + 1.7554e9 + cols = { + "temperature": 22 + 3 * np.sin(np.arange(n) / 288.0) + 0.1 * rng.normal(size=n), + "humidity": 50 + 9 * np.cos(np.arange(n) / 288.0), + "pressure": 1013 + 2 * np.sin(np.arange(n) / 600.0), + "lux": np.clip(400 * np.sin(np.arange(n) / 288.0), 0, None), + } + X = rng.normal(size=(n, 33)) + X[:, 0] = 1.0 + return ts, cols, X, np.ones(n, dtype=bool) + + +# ------------------------------------------------ a refit is not an outcome + +def test_refitting_does_not_touch_the_weights_or_the_calibrator(): + """The defect this file exists for. + + Measured on 8.2 days of a real station, fit() had put 977,078 Hedge updates + through the 15 minute head from 758 distinct supervised pairs, and 296,715 + through the 1 day head from 12. Hedge is multiplicative, so an edge far too + small to be real compounds to certainty. + """ + _, cols, X, valid = _synthetic() + ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + before = {k: h.weights.copy() for k, h in ens.heads.items()} + + for _ in range(5): + ens.fit(X, valid, cols) + + for k, head in ens.heads.items(): + assert np.allclose(head.weights, before[k]), \ + f"{k} blend weights moved during a refit" + assert len(head.conformal.scores) == 0, \ + f"{k} fed {len(head.conformal.scores)} scores to the calibrator from a refit" + assert head.n_scored == 0 + # The regression itself must still be learning. + assert head.model.n_updates > 0 + + +def test_an_outcome_does_move_them(): + """The other half: verify()'s channel has to work, or nothing ever learns.""" + _, cols, X, valid = _synthetic() + ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + ens.fit(X, valid, cols) + head = ens.heads[("temperature", 900)] + before = head.weights.copy() + + # Persistence right, everything else wrong, so the weights must move to it. + members = np.array([22.0, 30.0, 30.0, 30.0]) + for k in range(20): + head.observe_outcome(members, 22.0, covered=True, valid_ts=1e9 + k * 900) + + assert head.n_scored == 20 + assert len(head.conformal.scores) == 20 + assert head.weights[MEMBERS.index("persistence")] > before[MEMBERS.index("persistence")] + assert head.weights[MEMBERS.index("climatology")] < before[MEMBERS.index("climatology")] + + +# ------------------------------------------- overlapping outcomes are one fact + +def test_hedge_ignores_outcomes_that_overlap_the_last_one_it_took(): + """Forecasts are issued every retrain tick, so at the 1 day horizon about + two hundred a day resolve against very nearly the same outcome. The + conformal window can absorb that, being an order statistic. Exponentiated + gradient cannot.""" + _, cols, X, valid = _synthetic() + ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + head = ens.heads[("temperature", 86400)] + members = np.array([22.0, 30.0, 30.0, 30.0]) + + t0 = 1.7554e9 + for k in range(50): # 50 forecasts, 10 minutes apart + head.observe_outcome(members, 22.0, covered=True, valid_ts=t0 + k * 600) + + assert head.n_scored == 1, \ + f"Hedge took {head.n_scored} of 50 overlapping outcomes at a 1 d horizon" + # The calibrator still sees all of them: a quantile over duplicated scores + # is over-confident about its sample size, not wrong about its value. + assert len(head.conformal.scores) == 50 + + head.observe_outcome(members, 22.0, covered=True, valid_ts=t0 + 86400) + assert head.n_scored == 2, "a genuinely new outcome must be taken" + + +def test_the_decimation_clock_survives_a_restart(): + """Otherwise every deploy hands the weights a free duplicate.""" + ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + head = ens.heads[("pressure", 3600)] + head.observe_outcome(np.full(len(MEMBERS), 1013.0), 1013.5, + covered=True, valid_ts=1.7554e9) + + back = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + back.load_dict(ens.to_dict()) + assert back.heads[("pressure", 3600)].last_hedge_ts == head.last_hedge_ts + + # A head that has never scored must not refuse its first outcome. + fresh = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + fresh.load_dict(NowcastEnsemble(CONFIG.model.targets, + CONFIG.model.horizons_s, CONFIG.model).to_dict()) + h = fresh.heads[("pressure", 3600)] + h.observe_outcome(np.full(len(MEMBERS), 1013.0), 1013.5, + covered=True, valid_ts=1.7554e9) + assert h.n_scored == 1 + + +# -------------------------------------------------------------- persistence + +def test_forecast_members_round_trip(tmp_path): + """The learned member cannot be recomputed at maturity, the RLS has moved + on, so it has to be written down at issue time.""" + store = Store(str(tmp_path / "f.db")) + members = {"persistence": 21.0, "climatology": 22.0, + "learned": 23.0, "setpoint": 24.0} + store.insert_forecast(1000.0, 900, "temperature", 22.5, 21.0, 24.0, + "ensemble", members=members) + row = store.due_forecasts(now=2000.0)[0] + assert Station._members_of(row).tolist() == [21.0, 22.0, 23.0, 24.0] + assert row["scored"] == 0 + + store.mark_forecasts_scored([(1000.0, 900, "temperature")]) + assert store.due_forecasts(now=2000.0)[0]["scored"] == 1 + + +def test_a_forecast_from_before_the_columns_existed_is_not_fatal(tmp_path): + """It still deserves a conformal score. It just cannot move the weights.""" + store = Store(str(tmp_path / "old.db")) + store.insert_forecast(1000.0, 900, "temperature", 22.5, 21.0, 24.0, "ensemble") + row = store.due_forecasts(now=2000.0)[0] + assert Station._members_of(row) is None + + +def test_migration_adds_the_forecast_columns_to_a_live_table(tmp_path): + """Same hazard as the telemetry columns: CREATE TABLE IF NOT EXISTS is a + no-op against a table that already exists.""" + import sqlite3 + + path = str(tmp_path / "legacy.db") + with sqlite3.connect(path) as c: + c.execute("CREATE TABLE forecasts (" + "issued_ts REAL NOT NULL, valid_ts REAL NOT NULL, " + "horizon_s INTEGER NOT NULL, target TEXT NOT NULL, " + "mu REAL, lo REAL, hi REAL, model TEXT, " + "PRIMARY KEY (issued_ts, horizon_s, target))") + c.execute("INSERT INTO forecasts VALUES (1.0, 901.0, 900, 'temperature'," + " 22.0, 21.0, 23.0, 'ensemble')") + + store = Store(path) + with sqlite3.connect(path) as c: + have = {r[1] for r in c.execute("PRAGMA table_info(forecasts)")} + for col in ("m_persistence", "m_climatology", "m_learned", "m_setpoint", "scored"): + assert col in have, f"migration missed {col}" + + rows = store.due_forecasts(now=2000.0) + assert len(rows) == 1, "the existing row must survive the migration" + assert rows[0]["scored"] == 0 + assert Station._members_of(rows[0]) is None + + +# --------------------------------------------------------------- end to end + +def test_verify_teaches_from_a_matured_forecast_exactly_once(tmp_path): + """A matured forecast stays readable for an hour so the scorecard can + aggregate a rolling window, which means verify() sees it about twelve + times. Aggregating it twelve times is harmless. Teaching from it twelve + times is how the calibrator ends up believing it holds four hundred + independent scores when it holds one.""" + import time + + cfg = load_config() + cfg.storage.db_path = str(tmp_path / "verify.db") + # Without this the station loads the developer's own saved state, whose + # heads already carry six figures of n_scored. + cfg.storage.state_dir = str(tmp_path) + st = Station(cfg) + + now = time.time() + # The record has to run past the forecast's validity time, or verify() + # finds no truth to score it against. + for k in range(61): + ts = now - 1800 + k * 30 + st.store.insert_telemetry({"ts": ts, "temp_raw": 21.0, "temp_smooth": 21.0, + "hum": 50.0, "hum_smooth": 50.0, + "press": 1013.0, "press_slp": 1013.0}) + + issued = now - 1500 + st.store.insert_forecast(issued, 900, "temperature", 21.4, 20.9, 21.9, + "ensemble", + members={"persistence": 21.0, "climatology": 21.6, + "learned": 21.5, "setpoint": 21.0}) + + head = st.nowcast.heads[("temperature", 900)] + first = st.verify() + assert first["learned"] == 1, first + assert head.n_scored == 1 + assert len(head.conformal.scores) == 1 + + for _ in range(5): + again = st.verify() + assert again["learned"] == 0, "the same outcome was taught again" + assert head.n_scored == 1 + assert len(head.conformal.scores) == 1 + # and it is still being aggregated into the scorecard + assert again["scored"] >= 1 + + +def test_state_from_before_the_fix_is_not_trusted(): + """The saturated weights do not decay on their own. + + A pre-fix state file holds weights produced by the same week of weather + read about a thousand times. Hedge needs roughly twenty independent + outcomes to climb back off its 1e-4 floor, and the 1 d head sees one a day, + so leaving them in place would mean three weeks of a forecast pinned to + whichever member the replay happened to favour. + """ + ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + head = ens.heads[("temperature", 86400)] + for k in range(400): + head.observe_outcome(np.array([22.0, 30.0, 30.0, 30.0]), 22.0, + covered=False, valid_ts=1.7554e9 + k * 86400) + saved = ens.to_dict() + assert head.weights.max() > 0.9 and head.conformal.alpha < 0.02 + + for h in saved["heads"]: + h.pop("last_hedge_ts") # a state file from before + + back = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + back.load_dict(saved) + b = back.heads[("temperature", 86400)] + assert np.allclose(b.weights, 1.0 / len(MEMBERS)), "saturated weights were trusted" + assert b.n_scored == 0 + assert b.conformal.alpha == b.conformal.alpha_target, "the ACI integrator kept its wind-up" + # but the residual magnitudes are kept, or the long horizons spend nine + # days on the Gaussian fallback + assert len(b.conformal.scores) == 400 + + # A state file written after the fix must survive untouched. + keep = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model) + keep.load_dict(ens.to_dict()) + assert np.allclose(keep.heads[("temperature", 86400)].weights, head.weights)