Separate model accuracy from forecast accuracy

The Forecast Error sensors replayed the model with the relay, electrical heat
and window history that actually followed, so they measured the model, not a
forecast anyone could have made. They are renamed Model Error (unique IDs
unchanged) and gain RMSE, effective sample count and a quality label; their
state stays unknown until there are at least 5 roughly independent spans.

A new Forecast Error sensor scores forecasts as published, with the plan and
assumptions of the moment they were made. A forecast is scored only when the
relay followed its plan; otherwise it is counted as counterfactual, not as
error. Statistics are saved with the model and dropped with a zone's
recalibration.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

custom_components/heating_mpc/evaluation.py

@@ -9,18 +9,28 @@ then (so later data is never used), the outdoor temperature held at its starting
 value, and the radiator lag that actually followed from the relay's real ON/OFF
 history. The replayed temperature is compared with the reading now.
 
-This scores the thermal model, not the MPC's proposed plan: Generic Thermostat
-controls the relay, so the plan is never executed. Two baselines are scored on
-the same spans: persistence (no change) and the base 3-parameter model. The full
-model replays all zones together with the electrical heat, neighbour exchange and
-window openings that actually happened, so its improvement over the base model
-shows whether those features earn their keep.
+Two kinds of accuracy are kept apart:
+
+* ForecastEvaluator, **model accuracy** (retrospective): the model is replayed
+  with the relay, electrical heat and window history that actually happened. It
+  answers "how well does the model reproduce the rooms under known inputs?" and
+  compares the full model with the base model and with persistence. It is not a
+  forecast accuracy: a real forecast does not know the future inputs.
+* OperationalForecastLedger, **forecast accuracy**: forecasts as published, with
+  the boiler plan, electrical heat and windows they assumed, compared with what
+  happened once their horizon has passed. A forecast assuming the MPC's plan can
+  only be checked when the relay actually followed that plan; otherwise it is
+  counted as counterfactual and not scored.
+
+Forecasts start every 5 minutes, so consecutive ones overlap: n forecasts at a
+horizon of h minutes carry roughly n x 5 / h independent samples. Metrics report
+that number and a quality label, and sensors show no value while it is too small.
 """
 from __future__ import annotations
 
 from collections import deque
-from datetime import datetime, timedelta
-from math import isfinite
+from datetime import datetime, timedelta, timezone
+from math import isfinite, sqrt
 from typing import Any
 
 from .physics import Learner
@@ -30,6 +40,20 @@ from .zone_model import feature_vector
 HORIZONS = (15, 30, 60, 120)  # minutes
 WINDOW_DAYS = 7
 STEP_SECONDS = 300
+# Independent samples (overlap-adjusted) needed before a metric is shown / trusted.
+MIN_EFFECTIVE_SAMPLES = 5
+GOOD_EFFECTIVE_SAMPLES = 30
+
+
+def effective_samples(count: int, horizon: int, zones: int = 1) -> float:
+    """Roughly independent samples among `count` forecasts started 5 minutes apart."""
+    return round(count / max(1, zones) * min(1.0, STEP_SECONDS / 60 / horizon), 1)
+
+
+def quality(effective: float) -> str:
+    if effective < MIN_EFFECTIVE_SAMPLES:
+        return "insufficient"
+    return "indicative" if effective < GOOD_EFFECTIVE_SAMPLES else "ok"
 
 
 def replay(temp: float, outside: float, theta: tuple[float, float, float],
@@ -94,8 +118,8 @@ class ForecastEvaluator:
     def __init__(self, max_lag_uncertainty: float) -> None:
         self._max_lag_uncertainty = max_lag_uncertainty
         self._history: deque[dict[str, Any]] = deque(maxlen=HORIZONS[-1] * 60 // STEP_SECONDS + 1)
-        # UTC day -> horizon -> room ->
-        #   [count, sum error, sum |error|, sum |persistence error|, sum |base model error|]
+        # UTC day -> horizon -> room -> [count, sum error, sum |error|, sum |persistence error|,
+        #   sum |base model error|, sum error^2, count with squares (0 for 0.3/0.4 data)]
         self._days: dict[str, dict[str, dict[str, list[float]]]] = {}
 
     def reset_window(self) -> None:
@@ -146,12 +170,14 @@ class ForecastEvaluator:
                 error = full[room] - actual
                 base_error = replay(start, origin["outside"], theta, inputs) - actual
                 bucket = self._days.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(
-                    room, [0, 0.0, 0.0, 0.0, 0.0])
+                    room, [0, 0.0, 0.0, 0.0, 0.0, 0.0, 0])
                 bucket[0] += 1
                 bucket[1] += error
                 bucket[2] += abs(error)
                 bucket[3] += abs(start - actual)
                 bucket[4] += abs(base_error)
+                bucket[5] += error * error
+                bucket[6] += 1
 
     def _prune(self, time: datetime) -> None:
         oldest = (time.date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
@@ -166,12 +192,12 @@ class ForecastEvaluator:
             for day in self._days.values():
                 for room, values in day.get(str(horizon), {}).items():
                     if room in rooms:
-                        totals = per_room.setdefault(room, [0, 0.0, 0.0, 0.0, 0.0])
+                        totals = per_room.setdefault(room, [0, 0.0, 0.0, 0.0, 0.0, 0.0, 0])
                         for i, value in enumerate(values):
                             totals[i] += value
-            overall = [sum(values[i] for values in per_room.values()) for i in range(5)]
-            result[horizon] = {**_stats(overall),
-                               "rooms": {room: _stats(values) for room, values in per_room.items()}}
+            overall = [sum(values[i] for values in per_room.values()) for i in range(7)]
+            result[horizon] = {**_stats(overall, horizon, len(per_room)),
+                               "rooms": {room: _stats(values, horizon) for room, values in per_room.items()}}
         return result
 
     def forget(self, room: str) -> None:
@@ -184,7 +210,7 @@ class ForecastEvaluator:
 
     def as_dict(self) -> dict[str, Any]:
         # Rounded copy: shorter JSON on every save; in-memory sums stay exact.
-        return {"days": {day: {horizon: {room: [v[0], *(round(x, 6) for x in v[1:])]
+        return {"days": {day: {horizon: {room: [v[0], *(round(x, 6) for x in v[1:6]), v[6]]
                                          for room, v in by_room.items()}
                                for horizon, by_room in by_horizon.items()}
                          for day, by_horizon in self._days.items()}}
@@ -209,28 +235,38 @@ class ForecastEvaluator:
                 for room, values in by_room.items():
                     if _valid_bucket(room, values):
                         # 0.3 buckets had no base column: the model then was the base model.
-                        base = values[4] if len(values) == 5 else values[2]
+                        # Older buckets have no squared errors: they do not count for RMSE.
+                        base = values[4] if len(values) >= 5 else values[2]
+                        squares = (float(values[5]), values[6]) if len(values) == 7 else (0.0, 0)
                         self._days.setdefault(day, {}).setdefault(horizon, {})[room] = [
-                            values[0], float(values[1]), float(values[2]), float(values[3]), float(base)]
+                            values[0], float(values[1]), float(values[2]), float(values[3]), float(base),
+                            *squares]
 
 
 def _valid_bucket(room: Any, values: Any) -> bool:
-    if not isinstance(room, str) or not isinstance(values, list) or len(values) not in (4, 5):
+    if not isinstance(room, str) or not isinstance(values, list) or len(values) not in (4, 5, 7):
+        return False
+    count, *sums = values[:6]
+    if len(values) == 7 and (type(values[6]) is not int or not 0 <= values[6] <= count):
         return False
-    count, *sums = values
     return (type(count) is int and count >= 0
             and all(type(v) in (int, float) and isfinite(v) for v in sums)
             and all(v >= 0 for v in sums[1:]))
 
 
-def _stats(values: list[float]) -> dict[str, Any]:
-    count, error, absolute, persistence, base = values
+def _stats(values: list[float], horizon: int, zones: int = 1) -> dict[str, Any]:
+    count, error, absolute, persistence, base, squares, square_count = values
+    effective = effective_samples(int(count), horizon, zones)
     if not count:
-        return {"samples": 0, "mae": None, "bias": None, "persistence_mae": None, "skill": None,
+        return {"samples": 0, "effective_samples": 0.0, "quality": "insufficient", "mae": None,
+                "rmse": None, "bias": None, "persistence_mae": None, "skill": None,
                 "base_model_mae": None, "improvement_vs_base": None}
     return {
         "samples": int(count),
+        "effective_samples": effective,
+        "quality": quality(effective),
         "mae": round(absolute / count, 3),
+        "rmse": round(sqrt(squares / square_count), 3) if square_count else None,
         "bias": round(error / count, 3),
         "persistence_mae": round(persistence / count, 3),
         # Share of the persistence error removed by the model (1 = perfect, < 0 = worse).
@@ -240,3 +276,191 @@ def _stats(values: list[float]) -> dict[str, Any]:
         # ventilation (0 = no gain, < 0 = the extra features make it worse).
         "improvement_vs_base": round(1 - absolute / base, 3) if base > 0 else None,
     }
+
+
+# A forecast's plan counts as executed when, in every 15-minute block, the share
+# of time the relay was on is within this of what the plan assumed (0 or 1).
+PLAN_TOLERANCE = 0.2
+DUE_TOLERANCE_SECONDS = 150  # half a poll: the poll closest to the forecast's due time
+VERDICTS = ("counterfactual", "unverifiable", "missing")
+
+
+class OperationalForecastLedger:
+    """Forecasts as published, scored against what happened once their horizon passed.
+
+    Nothing is replaced with hindsight: the start temperatures, trajectory, boiler
+    plan and assumptions are those of the moment the forecast was made. Because the
+    plan is the MPC's proposal while Generic Thermostat controls the relay, a forecast
+    is only scored when the relay actually followed the plan up to that horizon;
+    otherwise it is counted as counterfactual. Pending forecasts live in memory;
+    the statistics are saved.
+    """
+
+    def __init__(self, saved: Any = None) -> None:
+        self._pending: deque[dict[str, Any]] = deque(maxlen=60)
+        # (interval end, seconds the relay was on or None if unknown, interval seconds)
+        self._relay: deque[tuple[float, float | None, float]] = deque(maxlen=48)
+        # UTC day -> horizon -> zone -> [count, sum error, sum |error|, sum error^2, sum |persistence error|]
+        self._days: dict[str, dict[str, dict[str, list[float]]]] = {}
+        # UTC day -> horizon -> verdict -> count of forecasts not scored, and why
+        self._verdicts: dict[str, dict[str, dict[str, int]]] = {}
+        self.load(saved)
+
+    def record(self, time: float, start: dict[str, float], trajectory: dict[str, list[float]],
+               plan: list[int], revision: str, assumptions: dict[str, Any]) -> None:
+        self._pending.append({"time": time, "start": dict(start), "trajectory": dict(trajectory),
+                              "plan": list(plan), "revision": revision, "assumptions": assumptions,
+                              "done": set()})
+
+    def observe(self, now: float, actual: dict[str, float | None], on_seconds: float | None,
+                seconds: float | None) -> None:
+        """One poll: the relay's last interval, then score forecasts that are due."""
+        if seconds:
+            self._relay.append((now, on_seconds, seconds))
+        day = datetime.fromtimestamp(now, timezone.utc).date().isoformat()
+        for forecast in list(self._pending):
+            for horizon in HORIZONS:
+                index = horizon // 15 - 1
+                due = forecast["time"] + horizon * 60
+                if horizon in forecast["done"] or now < due - DUE_TOLERANCE_SECONDS:
+                    continue
+                forecast["done"].add(horizon)
+                if index >= len(forecast["plan"]):
+                    continue  # beyond this forecast's horizon
+                if now > due + DUE_TOLERANCE_SECONDS:
+                    self._count(day, horizon, "missing")  # no poll near the due time
+                    continue
+                verdict = self._verify(forecast, horizon)
+                if verdict:
+                    self._count(day, horizon, verdict)
+                    continue
+                for zone, path in forecast["trajectory"].items():
+                    measured, start = actual.get(zone), forecast["start"].get(zone)
+                    if measured is None or start is None or len(path) <= index:
+                        self._count(day, horizon, "missing")
+                        continue
+                    error = path[index] - measured
+                    bucket = self._days.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(
+                        zone, [0, 0.0, 0.0, 0.0, 0.0])
+                    bucket[0] += 1
+                    bucket[1] += error
+                    bucket[2] += abs(error)
+                    bucket[3] += error * error
+                    bucket[4] += abs(start - measured)
+            if len(forecast["done"]) == len(HORIZONS) or now > forecast["time"] + (HORIZONS[-1] + 5) * 60:
+                self._pending.remove(forecast)
+        oldest = (datetime.fromtimestamp(now, timezone.utc).date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
+        for store in (self._days, self._verdicts):
+            for old in [d for d in store if d < oldest]:
+                del store[old]
+
+    def _verify(self, forecast: dict[str, Any], horizon: int) -> str | None:
+        """None when the relay followed the plan up to `horizon`, else the reason."""
+        for block in range(horizon // 15):
+            start = forecast["time"] + block * 900
+            end = start + 900
+            on = covered = 0.0
+            for interval_end, on_seconds, seconds in self._relay:
+                overlap = min(end, interval_end) - max(start, interval_end - seconds)
+                if overlap <= 0:
+                    continue
+                if on_seconds is None:
+                    return "unverifiable"
+                on += on_seconds * overlap / seconds
+                covered += overlap
+            if covered < 0.8 * 900:
+                return "unverifiable"
+            if abs(on / covered - forecast["plan"][block]) > PLAN_TOLERANCE:
+                return "counterfactual"
+        return None
+
+    def _count(self, day: str, horizon: int, verdict: str) -> None:
+        counts = self._verdicts.setdefault(day, {}).setdefault(str(horizon), {})
+        counts[verdict] = counts.get(verdict, 0) + 1
+
+    def forget(self, zone: str) -> None:
+        for by_horizon in self._days.values():
+            for by_zone in by_horizon.values():
+                by_zone.pop(zone, None)
+        for forecast in self._pending:
+            forecast["trajectory"].pop(zone, None)
+
+    @property
+    def pending(self) -> int:
+        return len(self._pending)
+
+    def summary(self, zones: list[str]) -> dict[int, dict[str, Any]]:
+        result = {}
+        for horizon in HORIZONS:
+            per_zone: dict[str, list[float]] = {}
+            verdicts = {verdict: 0 for verdict in VERDICTS}
+            for by_horizon in self._days.values():
+                for zone, values in by_horizon.get(str(horizon), {}).items():
+                    if zone in zones:
+                        totals = per_zone.setdefault(zone, [0, 0.0, 0.0, 0.0, 0.0])
+                        for i, value in enumerate(values):
+                            totals[i] += value
+            for by_horizon in self._verdicts.values():
+                for verdict, count in by_horizon.get(str(horizon), {}).items():
+                    verdicts[verdict] = verdicts.get(verdict, 0) + count
+            overall = [sum(values[i] for values in per_zone.values()) for i in range(5)]
+            result[horizon] = {**_operational_stats(overall, horizon, len(per_zone)), **verdicts,
+                               "zones": {zone: _operational_stats(values, horizon)
+                                         for zone, values in per_zone.items()}}
+        return result
+
+    def as_dict(self) -> dict[str, Any]:
+        return {"days": {day: {h: {zone: [v[0], *(round(x, 6) for x in v[1:])] for zone, v in by_zone.items()}
+                               for h, by_zone in by_horizon.items()} for day, by_horizon in self._days.items()},
+                "verdicts": self._verdicts}
+
+    def load(self, data: Any) -> None:
+        self._days, self._verdicts = {}, {}
+        if not isinstance(data, dict):
+            return
+        horizons = {str(h) for h in HORIZONS}
+        days = data.get("days") if isinstance(data.get("days"), dict) else {}
+        for day, by_horizon in days.items():
+            if not _valid_day(day) or not isinstance(by_horizon, dict):
+                continue
+            for horizon, by_zone in by_horizon.items():
+                if horizon not in horizons or not isinstance(by_zone, dict):
+                    continue
+                for zone, values in by_zone.items():
+                    if (isinstance(zone, str) and isinstance(values, list) and len(values) == 5
+                            and type(values[0]) is int and values[0] >= 0
+                            and all(type(v) in (int, float) and isfinite(v) for v in values[1:])
+                            and min(values[2:]) >= 0):
+                        self._days.setdefault(day, {}).setdefault(horizon, {})[zone] = [
+                            values[0], *(float(v) for v in values[1:])]
+        verdicts = data.get("verdicts") if isinstance(data.get("verdicts"), dict) else {}
+        for day, by_horizon in verdicts.items():
+            if not _valid_day(day) or not isinstance(by_horizon, dict):
+                continue
+            for horizon, counts in by_horizon.items():
+                if horizon in horizons and isinstance(counts, dict):
+                    clean = {k: v for k, v in counts.items() if k in VERDICTS and type(v) is int and v >= 0}
+                    if clean:
+                        self._verdicts.setdefault(day, {})[horizon] = clean
+
+
+def _valid_day(day: Any) -> bool:
+    try:
+        datetime.strptime(day, "%Y-%m-%d")
+        return True
+    except (TypeError, ValueError):
+        return False
+
+
+def _operational_stats(values: list[float], horizon: int, zones: int = 1) -> dict[str, Any]:
+    count, error, absolute, squares, persistence = values
+    effective = effective_samples(int(count), horizon, zones)
+    if not count:
+        return {"samples": 0, "effective_samples": 0.0, "quality": "insufficient", "mae": None, "rmse": None,
+                "bias": None, "persistence_mae": None, "improvement_vs_persistence": None}
+    return {
+        "samples": int(count), "effective_samples": effective, "quality": quality(effective),
+        "mae": round(absolute / count, 3), "rmse": round(sqrt(squares / count), 3),
+        "bias": round(error / count, 3), "persistence_mae": round(persistence / count, 3),
+        "improvement_vs_persistence": round(1 - absolute / persistence, 3) if persistence > 0 else None,
+    }