buzeraandClaude Opus 5.5
0b82a8e26c
Reproduced: a forecast planning 15 minutes of heat was scored although the relay was on for 12 (the 0.4.1 rule allowed 20 % per 15-minute block), and the same on-time placed elsewhere in a block passed as well. - The relay's on-time per poll interval, integrated at its own switch times, is compared with the plan block by block up to each horizon: the deviation is the time the relay spent in another state than planned. Unobserved relay time counts as possible deviation. - Matched (at most 2 minutes, worst case): Forecast Error figures. Approximate (at most 5 minutes, at most 2 of them unobserved): separate approximate_stats. Counterfactual (more than 5 minutes seen): counted, not scored. Incomplete otherwise, by reason: relay history, assumed inputs at forecast time (an unknown window or electrical heat), no poll near the due time, or a missing reading. - 2 minutes of unplanned boiler heat move a room by roughly half a 0.1 °C sensor step for the learned heat coefficients (0.05-0.12 °C per 5 minutes); 5 minutes by about one step. - Stored ledger statistics are versioned (2); 0.4.1 statistics used the looser rule and are dropped. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
608 lines
32 KiB
Python
608 lines
32 KiB
Python
"""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.
|
|
|
|
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, timezone
|
|
from math import isfinite, sqrt
|
|
from typing import Any
|
|
|
|
from .physics import Learner
|
|
from .windows import ventilation_known, ventilation_vector_for_prediction
|
|
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],
|
|
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
|
|
|
|
|
|
def _dot(theta: tuple[float, ...], x: list[float]) -> float:
|
|
return sum(t * v for t, v in zip(theta, x))
|
|
|
|
|
|
def model_snapshot(model: Any) -> dict[str, Any]:
|
|
"""Frozen parameters of a ZoneModel, so later learning never leaks into a replay."""
|
|
def layer(learner):
|
|
return (list(learner.features), tuple(learner.theta)) if learner is not None else None
|
|
return {"base": tuple(model.base.theta), "uses_extended": model.uses_extended,
|
|
"extended": layer(model.extended), "ventilation": layer(model.ventilation)}
|
|
|
|
|
|
def snapshot_delta(snapshot: dict[str, Any], temp: float, inputs: dict[str, Any]) -> float:
|
|
"""5-minute change predicted by a model snapshot (see zone_model.ZoneModel.delta)."""
|
|
main = None
|
|
if snapshot["uses_extended"] and snapshot["extended"]:
|
|
features, theta = snapshot["extended"]
|
|
x = feature_vector(features, temp, inputs)
|
|
main = _dot(theta, x) if x is not None else None
|
|
if main is None:
|
|
a, b, c = snapshot["base"]
|
|
main = a * (inputs["outside"] - temp) + b * inputs["lag"] + c
|
|
ventilation = inputs.get("ventilation")
|
|
if snapshot["ventilation"] and ventilation:
|
|
features, theta = snapshot["ventilation"]
|
|
x = ventilation_vector_for_prediction(features, temp, inputs["outside"], ventilation)
|
|
if x is not None:
|
|
main += _dot(theta, x)
|
|
return main
|
|
|
|
|
|
def inputs_complete(snapshot: dict[str, Any], temp: float, inputs: dict[str, Any],
|
|
incomplete: set[str]) -> bool:
|
|
"""Whether a step of a zone's model used only observed inputs.
|
|
|
|
Not when a window state, or a weather input an open window's terms need, was
|
|
unknown; nor when the extended layer is in use and its electrical heat or a
|
|
neighbour is unknown, or that neighbour's own replay is already incomplete.
|
|
The replay itself still falls back as forecasts do; the result is just not a
|
|
measurement of the model.
|
|
"""
|
|
vent_features = snapshot["ventilation"][0] if snapshot.get("ventilation") else ()
|
|
if not ventilation_known(vent_features, inputs.get("ventilation")):
|
|
return False
|
|
if snapshot["uses_extended"] and snapshot["extended"]:
|
|
features = snapshot["extended"][0]
|
|
if feature_vector(features, temp, inputs) is None:
|
|
return False
|
|
if any(f.startswith("couple:") and f.split(":", 1)[1] in incomplete for f in features):
|
|
return False
|
|
return True
|
|
|
|
|
|
def replay_zones(origin: dict[str, Any], path: list[dict[str, Any]]) -> tuple[dict[str, float], set[str]]:
|
|
"""All zones advanced together with the origin's models and the actual inputs.
|
|
|
|
Returns the replayed temperatures and the zones whose replay needed an input
|
|
that was not observed at some step (see inputs_complete).
|
|
"""
|
|
state = dict(origin["temps"])
|
|
incomplete: set[str] = set()
|
|
for before, poll in zip(path, path[1:]):
|
|
seconds = (poll["time"] - before["time"]).total_seconds()
|
|
new = {}
|
|
unknown = set(incomplete)
|
|
for zone, temp in state.items():
|
|
actual = poll["inputs"].get(zone, {})
|
|
inputs = {"outside": origin["outside"], "lag": poll["lag"], "temps": state,
|
|
"internal": actual.get("internal"), "ventilation": actual.get("ventilation")}
|
|
snapshot = origin["models"][zone]
|
|
if not inputs_complete(snapshot, temp, inputs, incomplete):
|
|
unknown.add(zone)
|
|
new[zone] = temp + snapshot_delta(snapshot, temp, inputs) * seconds / STEP_SECONDS
|
|
state, incomplete = new, unknown
|
|
return state, incomplete
|
|
|
|
|
|
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|,
|
|
# 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]]]] = {}
|
|
# UTC day -> horizon -> room -> [incomplete spans, skipped spans, sum |base model error|
|
|
# over the incomplete ones]. Neither kind enters the statistics above.
|
|
self._excluded: 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]],
|
|
models: dict[str, dict[str, Any]] | None = None,
|
|
inputs: dict[str, dict[str, Any]] | None = None) -> None:
|
|
"""Record a poll, then score every forecast that has reached its horizon.
|
|
|
|
`models`: model_snapshot per zone (base-only when omitted); `inputs`: this
|
|
poll's actual electrical heat and ventilation inputs per zone.
|
|
"""
|
|
if self._history and not 270 <= (time - self._history[-1]["time"]).total_seconds() <= 330:
|
|
self._history.clear() # irregular poll spacing: replays would be unreliable
|
|
models = models or {}
|
|
snapshots = {zone: models.get(zone) or {"base": tuple(theta[zone]), "uses_extended": False,
|
|
"extended": None, "ventilation": None}
|
|
for zone in temps if zone in theta}
|
|
self._history.append({"time": time, "temps": dict(temps), "outside": outside, "lag": lag,
|
|
"lag_uncertainty": lag_uncertainty, "theta": dict(theta),
|
|
"models": snapshots, "inputs": dict(inputs or {})})
|
|
self._score()
|
|
self._prune(time)
|
|
|
|
def _exclude(self, day: str, horizon: int, room: str, kind: int, base_error: float = 0.0) -> None:
|
|
bucket = self._excluded.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(room, [0, 0, 0.0])
|
|
bucket[kind] += 1
|
|
bucket[2] += abs(base_error)
|
|
|
|
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):]
|
|
origin = path[0]
|
|
# The actual heat input must be known along the way.
|
|
if any(poll["lag_uncertainty"] > self._max_lag_uncertainty for poll in path[1:]):
|
|
for room in origin["temps"]:
|
|
self._exclude(day, horizon, room, SKIPPED)
|
|
continue
|
|
inputs = [(poll["lag"], (poll["time"] - before["time"]).total_seconds())
|
|
for before, poll in zip(path, path[1:])]
|
|
full, incomplete = replay_zones(origin, path)
|
|
for room, start in origin["temps"].items():
|
|
actual, theta = latest["temps"].get(room), origin["theta"].get(room)
|
|
if actual is None or theta is None:
|
|
self._exclude(day, horizon, room, SKIPPED) # no reading at the end
|
|
continue
|
|
error = full[room] - actual
|
|
base_error = replay(start, origin["outside"], theta, inputs) - actual
|
|
if room in incomplete:
|
|
# The base model needs only the relay and outdoor temperature: its
|
|
# error is still a measurement, kept apart from the main figures.
|
|
self._exclude(day, horizon, room, INCOMPLETE, base_error)
|
|
continue
|
|
bucket = self._days.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(
|
|
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()
|
|
for store in (self._days, self._excluded):
|
|
for day in [day for day in store if day < oldest]:
|
|
del store[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, 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(7)]
|
|
excluded: dict[str, list[float]] = {}
|
|
for day in self._excluded.values():
|
|
for room, values in day.get(str(horizon), {}).items():
|
|
if room in rooms:
|
|
totals = excluded.setdefault(room, [0, 0, 0.0])
|
|
for i, value in enumerate(values):
|
|
totals[i] += value
|
|
all_excluded = [sum(values[i] for values in excluded.values()) for i in range(3)]
|
|
result[horizon] = {**_stats(overall, horizon, len(per_room)), **_excluded_stats(all_excluded),
|
|
"rooms": {room: {**_stats(per_room.get(room, [0, 0.0, 0.0, 0.0, 0.0, 0.0, 0]),
|
|
horizon),
|
|
**_excluded_stats(excluded.get(room, [0, 0, 0.0]))}
|
|
for room in rooms if room in per_room or room in excluded}}
|
|
return result
|
|
|
|
def forget(self, room: str) -> None:
|
|
"""Drop a zone's statistics, for example when its temperature is measured differently."""
|
|
for store in (self._days, self._excluded):
|
|
for by_horizon in store.values():
|
|
for by_room in by_horizon.values():
|
|
by_room.pop(room, None)
|
|
for poll in self._history:
|
|
poll["temps"].pop(room, None)
|
|
|
|
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: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()},
|
|
"excluded": {day: {horizon: {room: [v[0], v[1], round(v[2], 6)] for room, v in by_room.items()}
|
|
for horizon, by_room in by_horizon.items()}
|
|
for day, by_horizon in self._excluded.items()}}
|
|
|
|
def load(self, data: Any) -> None:
|
|
"""Restore stored statistics, keeping only well-formed entries."""
|
|
self._days, self._excluded = {}, {}
|
|
excluded = data.get("excluded") if isinstance(data, dict) else None
|
|
for day, horizon, room, values in _entries(excluded):
|
|
if (isinstance(values, list) and len(values) == 3 and all(type(v) is int and v >= 0 for v in values[:2])
|
|
and type(values[2]) in (int, float) and isfinite(values[2]) and values[2] >= 0):
|
|
self._excluded.setdefault(day, {}).setdefault(horizon, {})[room] = [
|
|
values[0], values[1], float(values[2])]
|
|
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):
|
|
# 0.3 buckets had no base column: the model then was the base model.
|
|
# 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),
|
|
*squares]
|
|
|
|
|
|
INCOMPLETE, SKIPPED = 0, 1 # columns of the excluded-span counts
|
|
|
|
|
|
def _entries(days: Any):
|
|
"""(day, horizon, zone, values) of stored per-day statistics with a valid shape."""
|
|
if not isinstance(days, dict):
|
|
return
|
|
horizons = {str(horizon) for horizon in HORIZONS}
|
|
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 in horizons and isinstance(by_zone, dict):
|
|
for zone, values in by_zone.items():
|
|
if isinstance(zone, str):
|
|
yield day, horizon, zone, values
|
|
|
|
|
|
def _excluded_stats(values: list[float]) -> dict[str, Any]:
|
|
incomplete, skipped, base = values
|
|
return {"incomplete": int(incomplete), "skipped": int(skipped),
|
|
"incomplete_base_model_mae": round(base / incomplete, 3) if incomplete else None}
|
|
|
|
|
|
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, 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
|
|
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], 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, "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).
|
|
"skill": round(1 - absolute / persistence, 3) if persistence > 0 else None,
|
|
"base_model_mae": round(base / count, 3),
|
|
# Share of the base model's error removed by electrical heat, neighbours and
|
|
# ventilation (0 = no gain, < 0 = the extra features make it worse).
|
|
"improvement_vs_base": round(1 - absolute / base, 3) if base > 0 else None,
|
|
}
|
|
|
|
|
|
# How far the relay may depart from a forecast's plan (seconds in another state than
|
|
# planned, up to the horizon) for the forecast to be compared with what happened.
|
|
# d seconds of boiler heat that the plan did not assume (or the reverse) move the
|
|
# room by roughly heat x d / 300, where heat is the learned boiler coefficient
|
|
# (typically 0.05-0.12 °C per 5 minutes with warm radiators). Up to 2 minutes keeps
|
|
# that below half a typical 0.1 °C sensor step and absorbs timestamp jitter;
|
|
# matched forecasts make the Forecast Error figures. Up to 5 minutes (one poll) is
|
|
# about one sensor step: kept apart as approximate. Beyond, the relay did
|
|
# something else: the forecast is counterfactual and is not scored. Relay time
|
|
# that was not observed counts as possible deviation, and more than 2 minutes of it
|
|
# leaves the forecast incomplete rather than approximate.
|
|
MATCH_SECONDS = 120
|
|
APPROXIMATE_SECONDS = 300
|
|
DUE_TOLERANCE_SECONDS = 150 # half a poll: the poll closest to the forecast's due time
|
|
VERDICTS = ("matched", "approximate", "counterfactual", "incomplete")
|
|
# Why a forecast (or, for a missing reading, a zone's forecast) could not be compared.
|
|
INCOMPLETE_REASONS = ("relay_history", "uncertain_inputs", "no_poll", "reading")
|
|
LEDGER_VERSION = 2 # 1 (0.4.1) accepted plans within 20 % per 15-minute block
|
|
|
|
|
|
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 compared only when the relay actually followed that plan (see
|
|
MATCH_SECONDS), and not at all when it was made with an input it had to
|
|
assume. 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]]]] = {} # matched forecasts
|
|
self._approximate: dict[str, dict[str, dict[str, list[float]]]] = {}
|
|
# UTC day -> horizon -> verdict (incomplete ones as "incomplete:<reason>") -> count
|
|
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
|
|
verdict = self._classify(forecast, horizon, now)
|
|
self._count(day, horizon, verdict)
|
|
if verdict not in ("matched", "approximate"):
|
|
continue
|
|
store = self._days if verdict == "matched" else self._approximate
|
|
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, "incomplete:reading")
|
|
continue
|
|
error = path[index] - measured
|
|
bucket = store.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._approximate, self._verdicts):
|
|
for old in [d for d in store if d < oldest]:
|
|
del store[old]
|
|
|
|
def _classify(self, forecast: dict[str, Any], horizon: int, now: float) -> str:
|
|
if now > forecast["time"] + horizon * 60 + DUE_TOLERANCE_SECONDS:
|
|
return "incomplete:no_poll"
|
|
if forecast["assumptions"].get("uncertain_inputs"):
|
|
return "incomplete:uncertain_inputs"
|
|
deviation, unobserved = self.deviation(forecast, horizon)
|
|
if deviation > APPROXIMATE_SECONDS:
|
|
return "counterfactual" # whatever happened in the unobserved part
|
|
# Otherwise judged by the worst case: unobserved time may all have deviated.
|
|
if deviation + unobserved <= MATCH_SECONDS:
|
|
return "matched"
|
|
if deviation + unobserved <= APPROXIMATE_SECONDS and unobserved <= MATCH_SECONDS:
|
|
return "approximate"
|
|
return "incomplete:relay_history"
|
|
|
|
def deviation(self, forecast: dict[str, Any], horizon: int) -> tuple[float, float]:
|
|
"""Seconds up to `horizon` the relay was seen in another state than planned,
|
|
and seconds its state is not known.
|
|
|
|
Compared per 15-minute plan block, from the relay's on-time per poll interval,
|
|
itself integrated at the relay's own switch times. So timing counts: the same
|
|
total on-time at other times than planned is a deviation, while a switch a few
|
|
seconds off a block boundary is not.
|
|
"""
|
|
start, end = forecast["time"], forecast["time"] + horizon * 60
|
|
deviation = covered = unobserved = 0.0
|
|
for interval_end, on_seconds, seconds in self._relay:
|
|
lo, hi = max(start, interval_end - seconds), min(end, interval_end)
|
|
if hi <= lo:
|
|
continue
|
|
covered += hi - lo
|
|
if on_seconds is None:
|
|
unobserved += hi - lo
|
|
continue
|
|
share = on_seconds / seconds
|
|
t = lo
|
|
while t < hi: # split at plan block boundaries
|
|
block = min(int((t - start) // 900), len(forecast["plan"]) - 1)
|
|
part_end = min(hi, start + (block + 1) * 900)
|
|
deviation += (part_end - t) * (share if not forecast["plan"][block] else 1.0 - share)
|
|
t = part_end
|
|
return deviation, unobserved + max(0.0, (end - start) - covered)
|
|
|
|
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 store in (self._days, self._approximate):
|
|
for by_horizon in store.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 = _sum_by_zone(self._days, horizon, zones)
|
|
approximate = _sum_by_zone(self._approximate, horizon, zones)
|
|
counts: dict[str, int] = {}
|
|
for by_horizon in self._verdicts.values():
|
|
for verdict, count in by_horizon.get(str(horizon), {}).items():
|
|
counts[verdict] = counts.get(verdict, 0) + count
|
|
reasons = {reason: counts.get(f"incomplete:{reason}", 0) for reason in INCOMPLETE_REASONS}
|
|
overall = [sum(values[i] for values in per_zone.values()) for i in range(5)]
|
|
approx_overall = [sum(values[i] for values in approximate.values()) for i in range(5)]
|
|
result[horizon] = {
|
|
**_operational_stats(overall, horizon, len(per_zone)),
|
|
# Forecasts by how the relay compared with their plan.
|
|
"matched": counts.get("matched", 0), "approximate": counts.get("approximate", 0),
|
|
"counterfactual": counts.get("counterfactual", 0), "incomplete": sum(reasons.values()),
|
|
"incomplete_reasons": reasons,
|
|
"approximate_stats": _operational_stats(approx_overall, horizon, len(approximate)),
|
|
"zones": {zone: _operational_stats(values, horizon) for zone, values in per_zone.items()},
|
|
}
|
|
return result
|
|
|
|
def as_dict(self) -> dict[str, Any]:
|
|
def rounded(store):
|
|
return {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 store.items()}
|
|
return {"version": LEDGER_VERSION, "days": rounded(self._days),
|
|
"approximate": rounded(self._approximate), "verdicts": self._verdicts}
|
|
|
|
def load(self, data: Any) -> None:
|
|
self._days, self._approximate, self._verdicts = {}, {}, {}
|
|
# 0.4.1 statistics used a looser plan criterion: they are not comparable.
|
|
if not isinstance(data, dict) or data.get("version") != LEDGER_VERSION:
|
|
return
|
|
for name, store in (("days", self._days), ("approximate", self._approximate)):
|
|
for day, horizon, zone, values in _entries(data.get(name)):
|
|
if (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):
|
|
store.setdefault(day, {}).setdefault(horizon, {})[zone] = [
|
|
values[0], *(float(v) for v in values[1:])]
|
|
known = {"matched", "approximate", "counterfactual", *(f"incomplete:{r}" for r in INCOMPLETE_REASONS)}
|
|
verdicts = data.get("verdicts") if isinstance(data.get("verdicts"), dict) else {}
|
|
horizons = {str(h) for h in HORIZONS}
|
|
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 known and type(v) is int and v >= 0}
|
|
if clean:
|
|
self._verdicts.setdefault(day, {})[horizon] = clean
|
|
|
|
|
|
def _sum_by_zone(store: dict[str, dict[str, dict[str, list[float]]]], horizon: int,
|
|
zones: list[str]) -> dict[str, list[float]]:
|
|
totals: dict[str, list[float]] = {}
|
|
for by_horizon in store.values():
|
|
for zone, values in by_horizon.get(str(horizon), {}).items():
|
|
if zone in zones:
|
|
bucket = totals.setdefault(zone, [0, 0.0, 0.0, 0.0, 0.0])
|
|
for i, value in enumerate(values):
|
|
bucket[i] += value
|
|
return totals
|
|
|
|
|
|
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,
|
|
}
|