Score the full model against the base model; compile the search step

Validation: the walk-forward evaluator now replays all zones together
with the full model (electrical heat, neighbour exchange, ventilation
and recovery) using the parameters frozen at the forecast start and the
electrical heat and window openings that actually happened, and scores
the base 3-parameter model on the same spans. Forecast Error sensors
add base_model_mae and improvement_vs_base (share of the base model's
error the new features remove; 0 = no gain, negative = worse). Saved
0.3 statistics load with their own error as the base column.

On synthetic apartments, 60-minute error full vs base: gaming office
0.031 vs 0.182 °C (83 % better), bedroom with ventilation 0.026 vs
0.218 °C (88 %), two coupled rooms 0.024 vs 0.037 °C (34 %; gone by
120 min, where the constant-outdoor assumption dominates). The same
ventilation week without window sensors shows no improvement, as it
should. These check the method, not real-world accuracy.

Performance: every term of a zone's 5-minute change is affine in the
zone temperatures and the radiator lag, and which layer applies depends
only on which inputs exist. compile_step_model reads the coefficients
off the real model per zone and step, so the plan search does a few
multiply-adds per step. It matches the direct model to 1e-12 and makes
a 180-minute, 5-zone search 4.5x faster (0.73 s here, close to the old
base-only search).

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

custom_components/heating_mpc/coordinator.py

@@ -21,7 +21,7 @@ from .const import (
     CONF_VENT_WIND, CONF_WINDOW_OPEN_MINUTES, DEFAULTS, READY_OFF_SAMPLES, READY_ON_SAMPLES, READY_SAMPLES,
 )
 from .data_quality import ReportAgeHistory, freshness
-from .evaluation import ForecastEvaluator
+from .evaluation import ForecastEvaluator, model_snapshot
 from .heat_sources import SUBENTRY_HEAT_SOURCE, HeatSource, HeatSourceTracker, build_sources, normalize_power
 from .learning import FeatureLearner
 from .physics import (
@@ -30,7 +30,7 @@ from .physics import (
 )
 from .weather_sources import choose_weather_sources, sources_for, normalize
 from .windows import SUBENTRY_WINDOW, VentilationTracker, Window, build_windows, facing_cos, ventilation_vector
-from .zone_model import ZoneModel, extended_features, feature_vector
+from .zone_model import ZoneModel, compile_step_model, extended_features, feature_vector
 from .zones import SUBENTRY_ZONE, Zone, build_zones, legacy_zone_data, zone_temperature
 
 _LOGGER = logging.getLogger(__name__)
@@ -289,13 +289,8 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             return None
         half_life = float(self.config.get(CONF_INTERNAL_HALF_LIFE, 30))
         forecast = {key: self._heat_sources.zone_forecast(key, steps, half_life) for key in temps}
-
-        def step(state: dict[str, float], lag: float, k: int) -> dict[str, float]:
-            k = min(k, steps - 1)
-            return {key: value + models[key].delta(value, {
-                "outside": outside, "lag": lag, "internal": forecast[key][k], "temps": state,
-                "ventilation": ventilation[key][k]}) for key, value in state.items()}
-        return step
+        return compile_step_model(models, list(temps), steps, lambda key, k: {
+            "outside": outside, "internal": forecast[key][k], "ventilation": ventilation[key][k]})
 
     def _zone_details(self, temps: dict[str, float], outside: float, internal: dict[str, float | None],
                       vent_inputs: dict[str, Any]) -> dict[str, Any]:
@@ -557,8 +552,12 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             "outside_source": active_sources["temperature"], "compositions": compositions,
         }
         # Parameters after this poll's learning: forecasts from here use only data so far.
-        self._evaluator.observe(now, vals, outside, self._lag["lag"], self._lag["uncertainty"],
-                                {room: tuple(self.learners[room].theta) for room in vals})
+        self._evaluator.observe(
+            now, vals, outside, self._lag["lag"], self._lag["uncertainty"],
+            {room: tuple(self.learners[room].theta) for room in vals},
+            models={room: model_snapshot(ZoneModel(self.learners[room], self.extended.get(room),
+                                                   self.ventilation.get(room))) for room in vals},
+            inputs={room: {"internal": internal[room], "ventilation": vent_inputs.get(room)} for room in vals})
         status["forecast_accuracy"] = self._evaluator.summary(self.rooms)
         if learned:
             # Persist at each sample so restart does not erase hours of training.

custom_components/heating_mpc/evaluation.py

@@ -10,8 +10,11 @@ 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.
+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.
 """
 from __future__ import annotations
 
@@ -21,6 +24,8 @@ from math import isfinite
 from typing import Any
 
 from .physics import Learner
+from .windows import ventilation_vector
+from .zone_model import feature_vector
 
 HORIZONS = (15, 30, 60, 120)  # minutes
 WINDOW_DAYS = 7
@@ -37,13 +42,60 @@ def replay(temp: float, outside: float, theta: tuple[float, float, float],
     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(features, temp, inputs["outside"], ventilation)
+        if x is not None:
+            main += _dot(theta, x)
+    return main
+
+
+def replay_zones(origin: dict[str, Any], path: list[dict[str, Any]]) -> dict[str, float]:
+    """All zones advanced together with the origin's models and the actual inputs."""
+    state = dict(origin["temps"])
+    for before, poll in zip(path, path[1:]):
+        seconds = (poll["time"] - before["time"]).total_seconds()
+        new = {}
+        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")}
+            new[zone] = temp + snapshot_delta(origin["models"][zone], temp, inputs) * seconds / STEP_SECONDS
+        state = new
+    return state
+
+
 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|]
+        # UTC day -> horizon -> room ->
+        #   [count, sum error, sum |error|, sum |persistence error|, sum |base model error|]
         self._days: dict[str, dict[str, dict[str, list[float]]]] = {}
 
     def reset_window(self) -> None:
@@ -51,12 +103,23 @@ class ForecastEvaluator:
         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."""
+                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)})
+                              "lag_uncertainty": lag_uncertainty, "theta": dict(theta),
+                              "models": snapshots, "inputs": dict(inputs or {})})
         self._score()
         self._prune(time)
 
@@ -75,17 +138,20 @@ class ForecastEvaluator:
             origin = path[0]
             inputs = [(poll["lag"], (poll["time"] - before["time"]).total_seconds())
                       for before, poll in zip(path, path[1:])]
+            full = replay_zones(origin, path)
             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
+                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])
+                    room, [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)
 
     def _prune(self, time: datetime) -> None:
         oldest = (time.date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
@@ -100,10 +166,10 @@ 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])
+                        totals = per_room.setdefault(room, [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(4)]
+            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()}}
         return result
@@ -134,23 +200,26 @@ class ForecastEvaluator:
                     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.
+                        base = values[4] if len(values) == 5 else values[2]
                         self._days.setdefault(day, {}).setdefault(horizon, {})[room] = [
-                            values[0], float(values[1]), float(values[2]), float(values[3])]
+                            values[0], float(values[1]), float(values[2]), float(values[3]), float(base)]
 
 
 def _valid_bucket(room: Any, values: Any) -> bool:
-    if not isinstance(room, str) or not isinstance(values, list) or len(values) != 4:
+    if not isinstance(room, str) or not isinstance(values, list) or len(values) not in (4, 5):
         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)
+            and all(v >= 0 for v in sums[1:]))
 
 
 def _stats(values: list[float]) -> dict[str, Any]:
-    count, error, absolute, persistence = values
+    count, error, absolute, persistence, base = values
     if not count:
-        return {"samples": 0, "mae": None, "bias": None, "persistence_mae": None, "skill": None}
+        return {"samples": 0, "mae": None, "bias": None, "persistence_mae": None, "skill": None,
+                "base_model_mae": None, "improvement_vs_base": None}
     return {
         "samples": int(count),
         "mae": round(absolute / count, 3),
@@ -158,4 +227,8 @@ def _stats(values: list[float]) -> dict[str, Any]:
         "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,
     }

custom_components/heating_mpc/sensor.py

@@ -164,13 +164,15 @@ class MPCForecastError(CoordinatorEntity, SensorEntity):
             "bias": accuracy.get("bias"),
             "persistence_mae": accuracy.get("persistence_mae"),
             "skill": accuracy.get("skill"),
+            "base_model_mae": accuracy.get("base_model_mae"),
+            "improvement_vs_base": accuracy.get("improvement_vs_base"),
             # Per zone, by name.
             "zones": {self.coordinator.zones[key].name if key in self.coordinator.zones else key: stats
                       for key, stats in accuracy.get("rooms", {}).items()},
             "window_days": WINDOW_DAYS,
             "unit_of_attributes": "°C",
-            "basis": ("thermal model replayed with the actual relay history; outdoor "
-                      "temperature held at the start value"),
+            "basis": ("all zones replayed together with the actual relay, electrical heat and "
+                      "window history; outdoor temperature held at the start value"),
         }
 
 

custom_components/heating_mpc/zone_model.py

@@ -13,7 +13,7 @@ samples where all of its inputs are known.
 """
 from __future__ import annotations
 
-from typing import Any, Mapping
+from typing import Any, Callable, Mapping
 
 from .const import READY_OFF_SAMPLES, READY_ON_SAMPLES, READY_SAMPLES
 from .learning import FeatureLearner
@@ -117,3 +117,44 @@ class ZoneModel:
         a, b, c = self.base.theta
         return {"loss": round(12 * a * (inputs["outside"] - temp), 3),
                 "heat": round(12 * b * inputs["lag"], 3), "bias": round(12 * c, 3)}
+
+
+def _probe(model: ZoneModel, key: str, inputs: Mapping[str, Any], state: Mapping[str, float], lag: float) -> float:
+    return model.delta(state[key], {**inputs, "temps": state, "lag": lag})
+
+
+def compile_step_model(models: Mapping[str, ZoneModel], keys: list[str], steps: int,
+                       inputs_at: Callable[[str, int], Mapping[str, Any]]):
+    """A fast, exact multi-zone step for the optimizer's plan search.
+
+    Every term of ZoneModel.delta is affine in the zone temperatures and the
+    radiator lag (loss, neighbours, ventilation and recovery all scale with
+    temperature differences), and which layer applies depends only on which
+    inputs exist, not on their values. So per zone and step the change is
+    a + sum_j b_j T_j + h lag, with the coefficients read off the real model at
+    a few probe points. The search then does a handful of multiply-adds per step
+    instead of rebuilding feature vectors for every candidate plan.
+    """
+    zero = {key: 0.0 for key in keys}
+    terms: dict[str, list[tuple[float, list[tuple[str, float]], float]]] = {}
+    for key in keys:
+        model = models[key]
+        related = [key] + ([name.split(":", 1)[1] for name in model.extended.features
+                            if name.startswith("couple:") and name.split(":", 1)[1] in zero]
+                           if model.extended is not None else [])
+        per_step = []
+        for k in range(steps):
+            inputs = dict(inputs_at(key, k))
+            a = _probe(model, key, inputs, zero, 0.0)
+            per_step.append((a, [(j, _probe(model, key, inputs, {**zero, j: 1.0}, 0.0) - a) for j in related],
+                             _probe(model, key, inputs, zero, 1.0) - a))
+        terms[key] = per_step
+
+    def step(state: dict[str, float], lag: float, k: int) -> dict[str, float]:
+        k = min(k, steps - 1)
+        new = {}
+        for key, value in state.items():
+            a, coefficients, h = terms[key][k]
+            new[key] = value + a + h * lag + sum(b * state[j] for j, b in coefficients)
+        return new
+    return step

tests/test_coordinator.py

@@ -784,6 +784,10 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
         self.assertEqual(details['model_in_use'], 'extended')
         self.assertGreater(details['contributions_c_per_h']['internal'], 0.2)
         self.assertTrue(details['coefficients']['extended']['internal']['identified'])
+        # Replayed against what happened, the electrical term removes most of the base error.
+        accuracy = data['forecast_accuracy'][60]
+        self.assertLess(accuracy['mae'], 0.06)
+        self.assertGreater(accuracy['improvement_vs_base'], 0.6)
 
     async def test_unknown_power_pauses_extended_learning_only(self):
         self.office_with_ups()
@@ -942,6 +946,11 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
         mixed = without.learners['bed'].theta
         self.assertLess(mixed[2], -0.008)
         self.assertGreater(mixed[1], 1.05 * b)
+        # Forecasts replayed with the actual window history are far better than the base model's.
+        configured = with_window._evaluator.summary(with_window.rooms)[60]
+        self.assertGreater(configured['improvement_vs_base'], 0.6)
+        self.assertLess(configured['mae'], 0.06)
+        self.assertEqual(without._evaluator.summary(without.rooms)[60]['improvement_vs_base'], 0)
 
     async def test_forecast_assumes_window_closes_then_room_recovers(self):
         self.c = self.bedroom_with_window()
@@ -1016,6 +1025,8 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
                 # ...but a steady difference between rooms is collinear with each room's
                 # bias, so the split stays uncertain and is not reported as identified.
                 self.assertFalse(info['identified'])
+        # Joint multi-zone forecasts still beat the base model at an hour.
+        self.assertGreater(c._evaluator.summary(c.rooms)[60]['improvement_vs_base'], 0.15)
 
     async def test_inter_zone_transfer_stays_near_zero_when_absent(self):
         c = await self.coupled_rooms(0.0, bed=(0.004, 0.11, 0.0))

tests/test_evaluation.py

@@ -76,6 +76,25 @@ class EvaluationTests(unittest.TestCase):
         self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 7)
         self.assertEqual(before, 7)
 
+    def test_full_model_replays_actual_inputs_and_is_compared_with_base(self):
+        base = (0.004, 0.12, 0.0005)
+        full = {'base': base, 'uses_extended': True,
+                'extended': (['loss', 'heat', 'bias', 'internal'], base + (0.15,)),
+                'ventilation': (['vent_natural'], (0.01,))}
+        evaluator, temp = ev.ForecastEvaluator(0.1), 20.0
+        for k in range(40):
+            internal = 0.4 if (k // 8) % 2 else 0.0           # gaming on and off
+            opening = 1.0 if k % 10 in (3, 4) else 0.0        # short ventilation
+            if k:
+                temp += 0.004 * (5 - temp) + 0.12 * 0.5 + 0.0005 + 0.15 * internal + 0.01 * opening * (5 - temp)
+            evaluator.observe(START + timedelta(seconds=300 * k), {'r': temp}, 5.0, 0.5, 0.0, {'r': base},
+                              models={'r': full},
+                              inputs={'r': {'internal': internal, 'ventilation': {'opening': opening}}})
+        stats = evaluator.summary(['r'])[60]
+        self.assertAlmostEqual(stats['mae'], 0, places=6)
+        self.assertGreater(stats['base_model_mae'], 0.05)
+        self.assertEqual(stats['improvement_vs_base'], 1.0)
+
     def test_saved_statistics_round_trip_and_reject_garbage(self):
         evaluator = plant(10)
         restored = ev.ForecastEvaluator(0.1)
@@ -86,8 +105,11 @@ class EvaluationTests(unittest.TestCase):
             'not-a-day': {'15': {'sensor.room': [1, 0.0, 0.1, 0.1]}},
             day: {'45': {'sensor.room': [1, 0.0, 0.1, 0.1]},
                   '15': {'sensor.a': [1, 0.0, -0.1, 0.1], 'sensor.b': [True, 0, 0, 0],
-                         'sensor.c': [1, float('nan'), 0, 0], 'sensor.d': [2, 0.1, 0.2, 0.4]}}}})
-        self.assertEqual(restored.as_dict(), {'days': {day: {'15': {'sensor.d': [2, 0.1, 0.2, 0.4]}}}})
+                         'sensor.c': [1, float('nan'), 0, 0], 'sensor.d': [2, 0.1, 0.2, 0.4],
+                         'sensor.e': [2, 0.1, 0.2, 0.4, 0.3], 'sensor.f': [2, 0.1, 0.2, 0.4, -1]}}}})
+        # A 0.3 bucket has no base-model column; its model was the base model then.
+        self.assertEqual(restored.as_dict(), {'days': {day: {'15': {
+            'sensor.d': [2, 0.1, 0.2, 0.4, 0.2], 'sensor.e': [2, 0.1, 0.2, 0.4, 0.3]}}}})
         for garbage in (None, [], {'days': []}):
             restored.load(garbage)
             self.assertEqual(restored.as_dict(), {'days': {}})

tests/test_zone_model.py

@@ -0,0 +1,87 @@
+"""Layered zone model: fallbacks, contributions and the compiled optimizer step."""
+import random
+import sys
+import types
+import unittest
+from pathlib import Path
+
+COMPONENTS = Path(__file__).resolve().parents[1] / 'custom_components' / 'heating_mpc'
+sys.modules.setdefault('heating_mpc', types.ModuleType('heating_mpc'))
+sys.modules['heating_mpc'].__path__ = [str(COMPONENTS)]
+from heating_mpc import learning, physics, zone_model as zm  # noqa: E402
+
+VENT = {'opening': 1.0, 'facing': -0.5, 'other_opening': 0.5, 'recovery': 0.0, 't_ref': 21.0,
+        'wind': 4.0, 'gust': 7.0}
+
+
+def trained(features, theta):
+    learner = learning.FeatureLearner(features)
+    learner.theta = list(theta)
+    learner.count, learner.on_count, learner.off_count = 300, 100, 100
+    return learner
+
+
+def models():
+    living = zm.ZoneModel(physics.Learner(), trained(zm.extended_features(True, ('bed',)),
+                                                     (0.004, 0.12, 0.0005, 0.15, 0.01)))
+    bed = zm.ZoneModel(physics.Learner(), trained(zm.extended_features(False, ('living',)),
+                                                  (0.006, 0.09, 0.0, 0.01)),
+                       trained(['vent_natural', 'recovery', 'vent_wind', 'vent_wind_dir', 'vent_cross'],
+                               (0.01, 0.15, 0.001, 0.0006, 0.002)))
+    return {'living': living, 'bed': bed}
+
+
+class ZoneModelTests(unittest.TestCase):
+    def test_layers_and_fallbacks(self):
+        self.assertEqual(zm.extended_features(False, ()), [])
+        self.assertEqual(zm.extended_features(True, ('b',)), ['loss', 'heat', 'bias', 'internal', 'couple:b'])
+        living = models()['living']
+        inputs = {'outside': 5.0, 'lag': 0.5, 'internal': 0.4, 'temps': {'living': 21.0, 'bed': 20.0}}
+        self.assertAlmostEqual(living.delta(21.0, inputs),
+                               0.004 * -16 + 0.06 + 0.0005 + 0.06 + 0.01 * -1)
+        # A missing input (electrical heat or a neighbour) falls back to the base model.
+        base = living.base.predict_delta(21.0, 5.0, 0.5)
+        self.assertAlmostEqual(living.delta(21.0, {**inputs, 'internal': None}), base)
+        self.assertAlmostEqual(living.delta(21.0, {**inputs, 'temps': {'living': 21.0}}), base)
+        # Until the extended model is as trained as the base needs to be, the base is used.
+        living.extended.count = 10
+        self.assertFalse(living.uses_extended)
+        self.assertAlmostEqual(living.delta(21.0, inputs), base)
+
+    def test_contributions_add_up_to_the_change(self):
+        bed = models()['bed']
+        inputs = {'outside': 5.0, 'lag': 0.5, 'temps': {'living': 21.0, 'bed': 19.0}, 'ventilation': VENT}
+        terms = bed.contributions(19.0, inputs)
+        self.assertEqual(set(terms), {'loss', 'heat', 'bias', 'neighbors', 'ventilation', 'recovery'})
+        self.assertAlmostEqual(sum(terms.values()), 12 * bed.delta(19.0, inputs), places=2)
+        self.assertLess(terms['ventilation'], 0)
+        # Unknown window state: no ventilation term rather than a guess.
+        self.assertAlmostEqual(bed.ventilation_delta(19.0, {**inputs, 'ventilation': {**VENT, 'opening': None}}), 0)
+
+    def test_compiled_step_matches_direct_model_exactly(self):
+        zone_models = models()
+        steps = 24
+        ventilation = [{**VENT, 'opening': 1.0 if k < 3 else 0.0, 'recovery': 0.0 if k < 3 else 0.6}
+                       for k in range(steps)]
+
+        def inputs_at(key, k):
+            return {'outside': 5.0, 'internal': 0.4 if key == 'living' else None,
+                    'ventilation': ventilation[k] if key == 'bed' else {}}
+        fast = zm.compile_step_model(zone_models, ['living', 'bed'], steps, inputs_at)
+        rng = random.Random(5)
+        for _ in range(300):
+            state = {'living': rng.uniform(15, 26), 'bed': rng.uniform(15, 26)}
+            lag, k = rng.random(), rng.randrange(steps + 3)  # beyond the horizon: last step repeats
+            direct = {key: value + zone_models[key].delta(value, {
+                **inputs_at(key, min(k, steps - 1)), 'lag': lag, 'temps': state}) for key, value in state.items()}
+            compiled = fast(state, lag, k)
+            for key in state:
+                self.assertAlmostEqual(compiled[key], direct[key], places=12)
+        # With a zone missing, its neighbour falls back to its base model, also when compiled.
+        solo = zm.compile_step_model(zone_models, ['living'], steps, inputs_at)
+        base = zone_models['living'].base.predict_delta(20.0, 5.0, 0.3)
+        self.assertAlmostEqual(solo({'living': 20.0}, 0.3, 0)['living'], 20.0 + base, places=12)
+
+
+if __name__ == '__main__':
+    unittest.main()