heating-mpc-ha/custom_components/heating_mpc/evaluation.py

"""Walk-forward evaluation of the thermal model against what actually happened.

Pure Python; no Home Assistant imports.

At every poll the readings, the outdoor temperature, the radiator-lag estimate
and each room's model parameters are recorded. When 15/30/60/120 minutes have
passed, each room is replayed from that earlier poll with the parameters it had
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. A persistence forecast (no
change) is scored alongside as the baseline to beat.
"""
from __future__ import annotations

from collections import deque
from datetime import datetime, timedelta
from math import isfinite
from typing import Any

from .physics import Learner

HORIZONS = (15, 30, 60, 120)  # minutes
WINDOW_DAYS = 7
STEP_SECONDS = 300


def replay(temp: float, outside: float, theta: tuple[float, float, float],
           steps: list[tuple[float, float]]) -> float:
    """Advance one room through (end-of-step lag, step seconds) with fixed parameters."""
    model = Learner(theta=list(theta))
    for lag, seconds in steps:
        # Same 5-minute model and jitter normalization as learning.
        temp += model.predict_delta(temp, outside, lag) * seconds / STEP_SECONDS
    return temp


class ForecastEvaluator:
    """Score replayed forecasts per room and horizon over the last WINDOW_DAYS days."""

    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|]
        self._days: dict[str, dict[str, dict[str, list[float]]]] = {}

    def reset_window(self) -> None:
        """Inputs were missing: no forecast may be scored across this poll."""
        self._history.clear()

    def observe(self, time: datetime, temps: dict[str, float], outside: float, lag: float,
                lag_uncertainty: float, theta: dict[str, tuple[float, float, float]]) -> None:
        """Record a poll, then score every forecast that has reached its horizon."""
        if self._history and not 270 <= (time - self._history[-1]["time"]).total_seconds() <= 330:
            self._history.clear()  # irregular poll spacing: replays would be unreliable
        self._history.append({"time": time, "temps": dict(temps), "outside": outside, "lag": lag,
                              "lag_uncertainty": lag_uncertainty, "theta": dict(theta)})
        self._score()
        self._prune(time)

    def _score(self) -> None:
        history = list(self._history)
        latest = history[-1]
        day = latest["time"].date().isoformat()
        for horizon in HORIZONS:
            steps = horizon * 60 // STEP_SECONDS
            if len(history) < steps + 1:
                continue
            path = history[-(steps + 1):]
            # The actual heat input must be known along the way.
            if any(poll["lag_uncertainty"] > self._max_lag_uncertainty for poll in path[1:]):
                continue
            origin = path[0]
            inputs = [(poll["lag"], (poll["time"] - before["time"]).total_seconds())
                      for before, poll in zip(path, path[1:])]
            for room, actual in latest["temps"].items():
                start, theta = origin["temps"].get(room), origin["theta"].get(room)
                if start is None or theta is None:
                    continue
                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])
                bucket[0] += 1
                bucket[1] += error
                bucket[2] += abs(error)
                bucket[3] += abs(start - actual)

    def _prune(self, time: datetime) -> None:
        oldest = (time.date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
        for day in [day for day in self._days if day < oldest]:
            del self._days[day]

    def summary(self, rooms: list[str]) -> dict[int, dict[str, Any]]:
        """Error statistics per horizon, overall and per configured room."""
        result = {}
        for horizon in HORIZONS:
            per_room: dict[str, list[float]] = {}
            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])
                        for i, value in enumerate(values):
                            totals[i] += value
            overall = [sum(values[i] for values in per_room.values()) for i in range(4)]
            result[horizon] = {**_stats(overall),
                               "rooms": {room: _stats(values) for room, values in per_room.items()}}
        return result

    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:])]
                                         for room, v in by_room.items()}
                               for horizon, by_room in by_horizon.items()}
                         for day, by_horizon in self._days.items()}}

    def load(self, data: Any) -> None:
        """Restore stored statistics, keeping only well-formed entries."""
        self._days = {}
        days = data.get("days") if isinstance(data, dict) else None
        if not isinstance(days, dict):
            return
        horizons = {str(horizon) for horizon in HORIZONS}
        for day, by_horizon in days.items():
            try:
                datetime.strptime(day, "%Y-%m-%d")
            except (TypeError, ValueError):
                continue
            if not isinstance(by_horizon, dict):
                continue
            for horizon, by_room in by_horizon.items():
                if horizon not in horizons or not isinstance(by_room, dict):
                    continue
                for room, values in by_room.items():
                    if _valid_bucket(room, values):
                        self._days.setdefault(day, {}).setdefault(horizon, {})[room] = [
                            values[0], float(values[1]), float(values[2]), float(values[3])]


def _valid_bucket(room: Any, values: Any) -> bool:
    if not isinstance(room, str) or not isinstance(values, list) or len(values) != 4:
        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 sums[1] >= 0 and sums[2] >= 0)


def _stats(values: list[float]) -> dict[str, Any]:
    count, error, absolute, persistence = values
    if not count:
        return {"samples": 0, "mae": None, "bias": None, "persistence_mae": None, "skill": None}
    return {
        "samples": int(count),
        "mae": round(absolute / count, 3),
        "bias": round(error / count, 3),
        "persistence_mae": round(persistence / count, 3),
        # Share of the persistence error removed by the model (1 = perfect, < 0 = worse).
        "skill": round(1 - absolute / persistence, 3) if persistence > 0 else None,
    }