Fix the runaway forecasts: refits accumulated, and annual terms fitted too early

Reported from a real station after 1.5 days: a six hour temperature forecast of
53 C in a 24 C room, and 9 C at one day, both carrying a plus or minus of 0.43.
Confidently wrong is the one failure this project is supposed to refuse.

Root cause. fit() replayed history into the live RLS on every retrain tick and
never reset, so 453 grid rows had produced 64,676 updates in a day and a half.
RLS with forgetting reads every update as fresh evidence, so the model believed
it had a hundred times the data it had: P collapsed, in-sample error looked
excellent, and the weights drifted without bound in directions the data never
excited. Measured: cond(P) 3.1e9 and ||theta|| 1680 against a median |theta| of
1.67. A refit now starts from the prior, which makes retraining idempotent.
Across 25 refits on the real data ||theta|| holds at 11.35, drifting 0.03, where
before it grew without limit.

The two largest weights were sin_doy and cos_doy at +1174 and +1191. Annual
harmonics were in the design matrix from the first sample, where they are
near-constant, near-collinear with each other and with the bias, and a
rank-deficient regressor is what RLS answers with enormous cancelling weights.
They are now held at zero until the record spans the same 120 days the
climatology fit already requires, because a day and a half of data says nothing
whatsoever about the season.

Also raised the standardiser's variance floor from 1e-8, which only caught a
bit-exactly constant column, to 1e-3. A feature that merely barely moves was
being divided by its own noise.

The conformal calibrators and Hedge weights are deliberately not reset by a
refit: those are earned from scored forecasts, not from this regression.

Backtest unchanged within noise, coverage still 89 to 91 across all 18 heads.
Four regression tests added, including that refitting the same history twice
must give the same model.

ashvale/features.py

@@ -82,7 +82,8 @@ def _rolling(a: np.ndarray, win: int, fn) -> np.ndarray:
 
 def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray,
                    press_slp: np.ndarray, lux: np.ndarray,
-                   grid_s: int, latitude: float, longitude: float
+                   grid_s: int, latitude: float, longitude: float,
+                   min_days_annual: float = 120.0
                    ) -> Tuple[np.ndarray, np.ndarray]:
     """Return (X of shape (n, N_FEATURES), valid mask of shape (n,))."""
     n = grid_ts.size
@@ -127,6 +128,20 @@ def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray,
     hour = (grid_ts % 86400.0) / 86400.0
     doy = (grid_ts % 31557600.0) / 31557600.0
 
+    # Annual harmonics are held at zero until the record spans enough of a year
+    # to excite them, exactly as the climatology fit already gates its annual
+    # terms. Left on from day one they are near-constant, near-collinear with
+    # each other and with the bias, and RLS answers that rank-deficient system
+    # with enormous cancelling weights. Measured on a real station after 1.5
+    # days: cos_doy +1191, sin_doy +1174, ||theta|| 1680 against a median |theta|
+    # of 1.67, cond(P) 3.1e9, and a six hour forecast of 53 C in a 24 C room.
+    # Zero is the honest value: with a day and a half of data the station knows
+    # nothing whatsoever about the season.
+    span_days = float(grid_ts[-1] - grid_ts[0]) / 86400.0 if n > 1 else 0.0
+    annual_on = 1.0 if span_days >= min_days_annual else 0.0
+    sin_doy = np.sin(2 * np.pi * doy) * annual_on
+    cos_doy = np.cos(2 * np.pi * doy) * annual_on
+
     X = np.column_stack([
         np.ones(n),
         temp, temp_rate_1h, temp_rate_3h, temp_std_3h, temp_dev_24h,
@@ -136,7 +151,7 @@ def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray,
         log_lux, cloud, elev, np.clip(elev, 0.0, None), (elev > 0.0).astype(float),
         np.sin(2 * np.pi * hour), np.cos(2 * np.pi * hour),
         np.sin(4 * np.pi * hour), np.cos(4 * np.pi * hour),
-        np.sin(2 * np.pi * doy), np.cos(2 * np.pi * doy),
+        sin_doy, cos_doy,
         press_anom * (hum - 70.0) / 100.0,
         press_tend_3h * dep,
     ])
@@ -172,7 +187,12 @@ class Standardiser:
         if self.n < 2:
             return np.atleast_2d(X)
         std = np.sqrt(self.m2 / max(self.n - 1, 1))
-        std = np.where(std < 1e-8, 1.0, std)
+        # 1e-8 was a token guard: it only catches a bit-exactly constant column.
+        # A feature that merely barely moves sails through and gets divided by
+        # its own noise, which manufactures a large z-score out of nothing. A
+        # feature with this little spread carries no information, so scale it by
+        # one and let it stay near zero rather than amplifying it.
+        std = np.where(std < 1e-3, 1.0, std)
         out = (np.atleast_2d(X) - self.mean) / std
         out[:, 0] = 1.0            # keep the bias column intact
         return out

ashvale/models/nowcast.py

@@ -179,6 +179,18 @@ class NowcastEnsemble:
         """Batch pass over history. Called on startup and every retrain tick."""
         if X.shape[0] < 10:
             return {"rows": 0}
+        # A refit starts from the prior. Without this, every retrain tick replays
+        # the same history into a live filter, and RLS with forgetting reads that
+        # as new evidence each time: measured on a real station after 1.5 days,
+        # 453 grid rows had produced 64,676 updates, cond(P) of 3.1e9 and a
+        # weight vector of norm 1680 whose two largest entries were the annual
+        # 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.
+        for head in self.heads.values():
+            head.model.reset()
         self.scaler.partial_fit(X[valid][:: max(1, X.shape[0] // 2000)])
         Xs = self.scaler.transform(X)
 

ashvale/models/rls.py

@@ -47,8 +47,9 @@ class RecursiveLeastSquares:
         self.d = int(n_features)
         self.lam = float(forgetting)
         self.p_max = float(p_max)
+        self.delta = float(delta)          # kept so a refit can return to the prior
         self.theta = np.zeros(self.d)
-        self.P = np.eye(self.d) * float(delta)
+        self.P = np.eye(self.d) * self.delta
         self.n_updates = 0
         self.ewma_sq_error = 0.0
 
@@ -102,14 +103,32 @@ class RecursiveLeastSquares:
     def noise_var(self) -> float:
         return float(max(self.ewma_sq_error, 1e-9))
 
+    def reset(self) -> None:
+        """Return to the prior, keeping the configuration.
+
+        A batch refit has to start from here rather than continuing, because
+        replaying the same history into a live filter is not the same as seeing
+        new data. RLS with forgetting treats every update as fresh evidence, so
+        feeding it the same rows on each retrain tick makes it believe it has
+        many times the data it has: P collapses, and the weights in directions
+        the data never excites drift without anything to pull them back.
+        """
+        self.theta = np.zeros(self.d)
+        self.P = np.eye(self.d) * self.delta
+        self.n_updates = 0
+        self.ewma_sq_error = 0.0
+
     def to_dict(self) -> Dict:
         return {"d": self.d, "lam": self.lam, "p_max": self.p_max,
+                "delta": self.delta,
                 "theta": self.theta.tolist(), "P": self.P.tolist(),
                 "n": self.n_updates, "ewma": self.ewma_sq_error}
 
     @classmethod
     def from_dict(cls, s: Dict) -> "RecursiveLeastSquares":
-        m = cls(s["d"], s["lam"], 1.0, s.get("p_max", 1e6))
+        # delta must survive the round trip or a refit after a restart would
+        # return to the wrong prior.
+        m = cls(s["d"], s["lam"], s.get("delta", 100.0), s.get("p_max", 1e6))
         m.theta = np.array(s["theta"], dtype=float)
         m.P = np.array(s["P"], dtype=float)
         m.n_updates = s.get("n", 0)

ashvale/station.py

@@ -505,6 +505,7 @@ class Station:
             grid_ts, cols["temperature"], cols["humidity"], cols["pressure"],
             cols["lux"], self.cfg.model.grid_s,
             self.cfg.site.latitude, self.cfg.site.longitude,
+            self.cfg.model.climatology_min_days_annual,
         )
         return grid_ts, cols, X, valid
 

scripts/evaluate.py

@@ -86,7 +86,8 @@ def main() -> None:
     )
     X, valid = build_features(grid_ts, cols["temperature"], cols["humidity"],
                               cols["pressure"], cols["lux"], cfg.model.grid_s,
-                              cfg.site.latitude, cfg.site.longitude)
+                              cfg.site.latitude, cfg.site.longitude,
+                              cfg.model.climatology_min_days_annual)
 
     n = grid_ts.size
     split = int(n * args.train_frac)

tests/test_models.py

@@ -149,3 +149,94 @@ def test_zambretti_rain_prior_rises_with_z():
     settled = zambretti(1035.0, 1.5, 6)
     stormy = zambretti(960.0, -2.5, 6)
     assert stormy["prior_rain_prob"] > settled["prior_rain_prob"]
+
+
+# ---------------------------------------------------------------- refit safety
+
+def test_rls_reset_returns_to_the_prior():
+    m = RecursiveLeastSquares(n_features=5, forgetting=0.999, delta=100.0)
+    rng = np.random.default_rng(9)
+    for _ in range(500):
+        m.update(rng.normal(size=5), float(rng.normal()))
+    assert m.n_updates == 500
+    m.reset()
+    assert m.n_updates == 0
+    assert np.allclose(m.theta, 0.0)
+    assert np.allclose(m.P, np.eye(5) * 100.0)
+
+
+def test_rls_delta_survives_serialisation():
+    """A refit after a restart must return to the same prior it started from."""
+    m = RecursiveLeastSquares(n_features=4, forgetting=0.99, delta=100.0)
+    m.update(np.ones(4), 1.0)
+    back = RecursiveLeastSquares.from_dict(m.to_dict())
+    back.reset()
+    assert np.allclose(back.P, np.eye(4) * 100.0), "reload lost the prior"
+
+
+def test_repeated_refits_do_not_accumulate():
+    """Refitting the same history must be idempotent, not cumulative.
+
+    This is the bug that put a 53 C six-hour forecast on a real station in a
+    24 C room. fit() replayed history into a live filter on every retrain tick
+    and never reset, so 453 grid rows had produced 64,676 updates in a day and a
+    half. RLS with forgetting reads each update as fresh evidence, so P
+    collapsed and the weights drifted without bound in the directions the data
+    never excited.
+    """
+    from ashvale.config import CONFIG
+    from ashvale.models.nowcast import NowcastEnsemble
+
+    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),
+        "pressure": 1013 + np.sin(np.arange(n) / 77.0),
+        "lux": np.clip(300 * np.sin(np.arange(n) / 120.0), 0, None),
+    }
+    X = rng.normal(size=(n, 33))
+    X[:, 0] = 1.0
+    valid = np.ones(n, dtype=bool)
+
+    ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model)
+    head = ens.heads[("temperature", 21600)]
+
+    ens.fit(X, valid, cols, None, ts)
+    first_norm = float(np.linalg.norm(head.model.theta))
+    first_updates = head.model.n_updates
+
+    for _ in range(15):
+        ens.fit(X, valid, cols, None, ts)
+
+    assert head.model.n_updates == first_updates, "updates accumulated across refits"
+    assert float(np.linalg.norm(head.model.theta)) == pytest.approx(first_norm, rel=0.05)
+
+
+def test_annual_harmonics_are_zero_until_the_record_spans_a_season():
+    """Two near-constant, near-collinear columns are a rank-deficient regressor.
+
+    Left on from day one, sin_doy and cos_doy carried +1174 and +1191 on a real
+    station whose median weight was 1.67. Zero is the honest value: a day and a
+    half of data says nothing whatsoever about the season.
+    """
+    from ashvale.features import FEATURE_NAMES, build_features
+
+    n = 450
+    ts = np.arange(n) * 300.0 + 1.7554e9          # about 1.5 days
+    t = 22 + 2 * np.sin(np.arange(n) / 40.0)
+    h = 50 + 5 * np.cos(np.arange(n) / 33.0)
+    p = 1013 + np.sin(np.arange(n) / 77.0)
+    lux = np.clip(300 * np.sin(np.arange(n) / 120.0), 0, None)
+
+    si, ci = FEATURE_NAMES.index("sin_doy"), FEATURE_NAMES.index("cos_doy")
+
+    X, _ = build_features(ts, t, h, p, lux, 300, 52.2, 0.12, min_days_annual=120.0)
+    assert np.all(X[:, si] == 0.0) and np.all(X[:, ci] == 0.0)
+
+    # A record that does span the year keeps them.
+    ts_long = np.arange(n) * (200 * 86400.0 / n) + 1.7554e9
+    X2, _ = build_features(ts_long, t, h, p, lux, 300, 52.2, 0.12, min_days_annual=120.0)
+    assert X2[:, si].std() > 0.1, "annual terms should return once the record is long enough"