diff --git a/README.md b/README.md index ecc0921..6d3a6da 100644 --- a/README.md +++ b/README.md @@ -68,9 +68,11 @@ systemd service on an always-on box, see [Running wallmonitor](https://github.co - **Automatic derate prevention** — an optional daemon caps the vehicle's charge current through a least-privilege BLE pairing when the forecast firms up, and restores it once the risk clears -- **Degradation watch** — ambient-corrected heat-rise trend with confidence - intervals and a verified-baseline anchor, so "getting worse" is a - statistical claim, not a vibe +- **Degradation watch** — heat rise regressed on time while ambient and + charge current are held, over windows the charger was not regulating, + with confidence intervals and a verified-baseline anchor, so "getting + worse" is a statistical claim about the connector and not about the + weather - **Actionable notifications, local-only** — browser push and LAN webhook/self-hosted ntfy for phones, with a systemd + Docker deploy recipe - Resilience: seamless restarts, downtime recorded as explicit gap events, diff --git a/docs/thermal-model.md b/docs/thermal-model.md index e5eeaff..97e8bc2 100644 --- a/docs/thermal-model.md +++ b/docs/thermal-model.md @@ -58,6 +58,31 @@ itself. The unit of thermal analysis is the **load window** — the stretch where current actually flows. +## Free-running windows, and the ones the charger wrote + +The other way a plateau can be fictitious is that the charger chose it. A +Gen 3 defends its own thermal limit long before alert 40: as the handle +warms it trims charge current back, and it can trim ~10 % without ever +leaving the fitter's steady-current band. The ramp then flattens *because +the current fell*, and the exponential reads that flattening as the +plateau — a lower rise paired with a faster τ, clearing every other gate +with an excellent RMSE. + +That bias is not random. Foldback starts sooner in a hot garage, so the +under-read arrives and leaves with the weather; and a session that runs at +a low enough current never triggers it at all. Measured on one install: +windows whose current sagged fitted a median +33.3 °C rise where the same +charger's steady windows fitted +37.2 °C. + +So every fit records `current_sag_a` — how far current fell from the head +of the window to its tail, compared by quarter-medians — and +`free_plateau`, true when that sag stayed inside 1.5 % of the window's +current. On the install above the two populations did not overlap: steady +windows sagged ≤ 0.6 %, regulated ones ≥ 3.9 %. + +Regulated fits still describe what the handle actually did, so the +forecast keeps them. The degradation watch cannot use them at all. + ## Ambient is a bracket, not a point A single start-of-window ambient silently assumes the weather held still for @@ -120,35 +145,77 @@ per-segment fit, and the drift verdict. ## Degradation watch The same per-segment fits feed a trend. Rising heat at unchanged current -means added resistance — a loose lug, a degrading contact — so when recent -segments' fitted rise climbs past the baseline, the poller raises a monitor -alert and the Alerts page charts the fitted-rise trend. +means added resistance — a loose lug, a degrading contact — so when the +fitted rise climbs over time, the poller raises a monitor alert and the +Alerts page charts the fitted-rise trend. + +### What it estimates, and why not a median split + +The watch **regresses fitted rise on time** across the whole comparable +history, holding ambient and charge current, and reads the *time* +coefficient. The reported Δ is that slope times the observed span. + +It did once compare a recent median against a baseline median, and that +asks the wrong question. "Are the last few fits higher?" is answered for +you by anything that moved with the calendar: a garage that cooled between +the two halves, or a vehicle capped to a lower current whose (48/I)² +normalization then lifts every recent fit at once. On one install the +split reported **+7.2 °C with a 95 % CI of [5.4, 9.1]** — "statistically +confirmed" — for a connector whose rise, regressed on time with ambient and +current held, was moving +0.01 ± 0.04 °C/day. The confidence was real. It +was confidence in the wrong estimand. + +Regression fixes three things at once: + +- **Confounders become covariates.** Ambient and charge current are + adjusted for instead of assumed away, and each one's coefficient is + reported so you can see what it was worth. +- **Every fit counts.** The estimate uses the whole history rather than + three fits against a handful, which is where the precision comes from. +- **Unseparable confounds declare themselves.** When a covariate cannot be + told apart from the calendar — a current cap applied once and kept is + nearly collinear with time — the collinearity inflates the slope's + standard error and the verdict declines to confirm. That is the honest + outcome, and it is reached automatically rather than by a rule someone + had to anticipate. + +A covariate earns a column only when the history actually moved in it +(≥ 3 °C of ambient, ≥ 2 A of current); below that it buys nothing and +spends a degree of freedom. `/api/thermal` reports which columns were used +(`covariates`), each coefficient, the residual scatter, and any covariate +that correlated with time past 0.8 (`collinear_with_time`). ### What counts as drift -The alert needs two things at once: the recent-vs-baseline delta must be -**material** (≥ 2.5 °C — a confirmed 0.3 °C increase is real but not worth -an inspection) and **confirmed** — its 95% confidence interval, built from -this install's own session-to-session scatter, must clear zero. The -effective alert threshold is therefore the larger of the floor and what the -scatter demands, and the dashboard shows which one is binding. A noisy -install (variable ambient, a sensor in a draughty spot) must show more -before the watch alarms; a quiet one, less. A fixed 2.5 °C tripwire sat -near one sigma on a real install and fired on scatter alone. +The alert needs the change to be **material** (≥ 2.5 °C — a confirmed +0.3 °C increase is real but not worth an inspection), **confirmed** — the +slope's 95 % confidence interval, at a small-sample Student-t multiplier, +must clear zero — and **soundly measured** (see the ambient confound +below). The effective threshold is the larger of the floor and what this +install's own scatter demands, and the dashboard shows which one is +binding. A noisy install must show more before the watch alarms; a quiet +one, less. -A delta past the floor whose interval still straddles zero is a **lead**: -shown on the dashboard, pushed once at default priority, no alert row. More -sessions either confirm it or dissolve it. +A change past the floor that fails either of the other two tests is a +**lead**: shown on the dashboard, pushed once at default priority, no alert +row. More sessions either confirm it or dissolve it. ### What it compares +- **Only free-running fits.** Windows the charger was regulating are + excluded outright — their plateau is a setpoint, not an equilibrium. + `regulated_n` says how many sat out, and the Alerts page says so too, + because a verdict resting on four fits should not look like one resting + on twelve. - **Only sessions near the install's recent operating current.** Cap the vehicle at a new amperage and the watch follows, rather than judging forever against a current the install no longer uses. - **Pooled across a wider current band when the fits are clean.** - Ambient-bracketed fits are clean enough under the I² normalization to pool - in, so a baseline recorded at 48 A keeps judging charges after a cap to - 40 A instead of the verdict going dark. + Ambient-bracketed fits join from a wider band, and the regression's own + current term then *adjusts* them: residual error in the I² normalization + lands on that coefficient instead of masquerading as a trend. On the + install above that coefficient read −0.99 °C per amp, which is the whole + of the phantom +7.2 °C. ### How sure it is @@ -168,11 +235,43 @@ healthy", not "vs the first charges the monitor happened to see". ### The confounder the fits can't remove -A **rise-vs-ambient scatter** on the same page separates the one thing -left. Ambient is subtracted per fit, so a healthy install shows a flat cloud -regardless of garage temperature: - -- a cloud **still sloping upward with ambient** exposes an environment - effect the model doesn't carry — multi-day heat soak of cable and - structure in an uninsulated garage; +The premise of `rise_ref` is that subtracting ambient leaves a number that +depends on the connector and not on the weather. The watch checks that +premise instead of assuming it: the regression's **ambient coefficient** +is how much rise still moves per degree of garage air after the +subtraction. A healthy install reads flat. + +An install reading materially non-zero (≥ 0.3 °C/°C, resolved well enough +to be sure of the sign) has something in the measurement that the model +does not carry, and the watch says so and **caps its verdict at a lead** — +it can never raise an alert until the cause is found. Adjusting for a +confounder is not the same as understanding it. + +The **rise-vs-ambient scatter** on the same page shows the shape: + +- a cloud **sloping upward with ambient** exposes an environment effect — + multi-day heat soak of cable and structure in an uninsulated garage; +- a cloud **sloping downward** is the signature of a handle held at a + fixed temperature by something — most often charger regulation that the + free-plateau gate did not catch, since a hotter garage means the trim + starts sooner; - an **elevated-but-flat** cloud is the genuine added-resistance signature. + +### What the watch cannot do, and what to charge to fix it + +Everything above measures a plateau the charger allowed. On an install +where full-rate charging always ends in foldback, the only free-running +windows are the low-current ones, and the watch is judging a handful of +fits at a current the install rarely uses. + +The cheapest fix is a **fixed-condition probe**: once a month, charge at a +current low enough to run unregulated end to end (well under whatever +first triggers foldback), for at least 3 τ. That yields a plateau nobody +imposed, at a repeatable operating point, and comparing those month over +month is a degradation test with no extrapolation in it at all. A single +such charge is worth more to this watch than several at full rate. + +Deliberately varying the current — 32 / 40 / 48 A inside one week, at +similar ambient — is worth doing once for a different reason: it breaks the +collinearity between current and the calendar, which is the one thing that +stops the regression from separating a cap from a trend. diff --git a/tests/test_wallmonitor.py b/tests/test_wallmonitor.py index ba301a6..948b82b 100644 --- a/tests/test_wallmonitor.py +++ b/tests/test_wallmonitor.py @@ -3,6 +3,7 @@ import asyncio import math import time +from statistics import median import aiohttp import pytest @@ -334,7 +335,7 @@ def _seed_idle(db, t_from, t_to, ambient_c, dt=10.0): def _seed_thermal_session(db, start_ts, ambient_c, tau_s=720.0, rise_ref_c=36.0, amps=48.6, charge_s=1500.0, dt=10.0, - ambient_end_c=None, cooldown_s=0.0): + ambient_end_c=None, cooldown_s=0.0, sag_to_a=None): """Idle lead-in plus a charging ramp that follows the first-order model. With ambient_end_c set, ambient drifts linearly across the charge (the @@ -342,24 +343,35 @@ def _seed_thermal_session(db, start_ts, ambient_c, tau_s=720.0, rise_ref_c=36.0, integrating the lag ODE against the moving ambient. cooldown_s appends post-session idle decay samples — the tail the fitter reads the load window's end ambient from. + + With sag_to_a set, the charge current falls linearly to that value + across the window — the charger's own thermal regulation, which flattens + the ramp for a reason that has nothing to do with the connector. Kept + inside the steady-prefix band on purpose: that is what makes it + dangerous, and what the free-plateau gate exists to catch. """ _seed_idle(db, start_ts - 1800, start_ts, ambient_c, dt) sid = db.start_session(start_ts) t0_temp = thermal.idle_handle_c(ambient_c) rise_at = rise_ref_c * (amps / thermal.REF_CURRENT_A) ** 2 + integrate = ambient_end_c is not None or sag_to_a is not None temp = t0_temp ts = start_ts while ts <= start_ts + charge_s: - if ambient_end_c is None: + elapsed = (ts - start_ts) / charge_s + amps_now = amps if sag_to_a is None else amps + (sag_to_a - amps) * elapsed + if not integrate: t_inf = ambient_c + rise_at temp = t_inf - (t_inf - t0_temp) * math.exp(-(ts - start_ts) / tau_s) db.insert_vitals(ts, { - "vehicle_connected": 1, "contactor_closed": 1, "vehicle_current_a": amps, + "vehicle_connected": 1, "contactor_closed": 1, "vehicle_current_a": round(amps_now, 1), "handle_temp_c": round(temp, 3), "pcba_temp_c": 55.0, "mcu_temp_c": 50.0, - }, sid, amps * 233.0) - if ambient_end_c is not None: - ambient_now = ambient_c + (ambient_end_c - ambient_c) * (ts - start_ts) / charge_s - temp += dt * ((ambient_now + rise_at - temp) / tau_s) + }, sid, amps_now * 233.0) + if integrate: + ambient_now = ambient_c if ambient_end_c is None else \ + ambient_c + (ambient_end_c - ambient_c) * elapsed + rise_now = rise_ref_c * (amps_now / thermal.REF_CURRENT_A) ** 2 + temp += dt * ((ambient_now + rise_now - temp) / tau_s) ts += dt db.close_session(sid, start_ts + charge_s, "vehicle_disconnected") ambient_final = ambient_end_c if ambient_end_c is not None else ambient_c @@ -991,7 +1003,11 @@ async def test_thermal_drift_confidence_interval(db): assert drift["drifting"] is True and drift["confident"] is True ci_lo, ci_hi = drift["delta_ci95_c"] assert ci_lo < drift["delta_c"] < ci_hi and ci_lo > 0 - assert drift["baseline_mad_c"] < 1.0 and drift["recent_mad_c"] < 1.0 + # One regression, so one scatter: how far the fits sit from the modelled + # trend, which is what the interval above is actually built from. A step + # change fitted by a straight line leaves residual by construction — the + # scatter here is mostly that, not measurement noise. + assert drift["resid_sd_c"] < 2.0 and drift["n"] == 7 async def test_thermal_drift_wide_scatter_is_not_confident(db): @@ -1019,7 +1035,12 @@ def fits(rises): for i, r in enumerate(rises)] quiet = thermal.detect_drift(fits([36.0, 36.2, 35.9, 36.1, 36.0, 35.8, 39.0, 39.2, 38.9])) noisy = thermal.detect_drift(fits([33.0, 39.0, 34.0, 38.0, 35.0, 37.0, 39.0, 39.2, 38.9])) - assert abs(quiet["delta_c"] - 3.0) < 0.2 and abs(noisy["delta_c"] - 3.0) < 0.6 + # delta_c is the modelled change across the observed span, not a step + # between two medians: a late 3 C jump fitted with a straight line reads + # a little larger than the step itself, and the same jump buried in + # scatter reads larger still. Both stay in the same neighbourhood — what + # separates them is the interval, below. + assert 3.0 < quiet["delta_c"] < 4.5 and 3.0 < noisy["delta_c"] < 5.5 assert quiet["drifting"] is True and quiet["lead"] is False assert abs(quiet["threshold_c"] - thermal.DRIFT_WARN_C) < 0.01 # the floor binds assert noisy["drifting"] is False and noisy["lead"] is True @@ -1056,6 +1077,129 @@ async def test_thermal_drift_pools_bracketed_cross_current_fits(db): assert drift is not None and drift["off_current_n"] >= 1 +async def test_thermal_fit_flags_regulated_windows(db): + # The Gen 3 defends its own thermal limit by trimming charge current as + # the handle warms, and it can trim ~10% without ever leaving the + # steady-prefix band. The ramp then flattens because the current fell, + # and the exponential reads that flattening as the plateau. Each fit has + # to say whether its window was left alone. + now = time.time() + steady = _seed_thermal_session(db, now - 4 * 7200, ambient_c=25.0) + folded = _seed_thermal_session(db, now - 2 * 7200, ambient_c=25.0, sag_to_a=43.5) + fits = {fit["session_id"]: fit for fit in thermal.fit_sessions(db, now)} + assert fits[steady]["free_plateau"] is True + assert abs(fits[steady]["current_sag_a"]) < 0.5 + assert fits[folded]["free_plateau"] is False + assert fits[folded]["current_sag_a"] > 1.0 + # And the regulated window really is biased low, which is why it cannot + # be allowed near a baseline: same seeded connector, smaller answer. + assert fits[folded]["rise_ref_c"] < fits[steady]["rise_ref_c"] - 1.0 + + +async def test_thermal_drift_excludes_regulated_windows(db): + # Regulated windows are biased low, so a history that starts regulated + # and ends free-running manufactures a rise out of nothing. (On a real + # install this is seasonal: foldback starts sooner in a hot garage, so + # the bias arrives and leaves with the weather.) Excluding them is what + # makes the remaining fits comparable to each other at all. + now = time.time() + for i in range(4): + _seed_thermal_session(db, now - (10 - i) * 7200, ambient_c=25.0, sag_to_a=43.5) + for i in range(6): + _seed_thermal_session(db, now - (6 - i) * 7200, ambient_c=25.0) + fits = thermal.fit_sessions(db, now) + assert sum(1 for fit in fits if not fit["free_plateau"]) == 4 + drift = thermal.detect_drift(fits) + assert drift is not None + assert drift["regulated_n"] == 4 and drift["n"] == 6 + assert drift["drifting"] is False and drift["lead"] is False + assert abs(drift["delta_c"]) < thermal.DRIFT_WARN_C + + +async def test_thermal_drift_holds_charge_current(db): + # The install's connector is unchanged; the vehicle simply starts asking + # for 39.6 A instead of 48.6. The (48/I)^2 normalization does not + # describe any real install exactly, so every reduced-current fit + # normalizes ~7 C high — and a recent-vs-baseline median split reports + # that as a confirmed degradation. Holding current in the regression + # reads it as what it is: a level shift with current, not a trend. + now = time.time() + for i, rise in enumerate([36.0, 36.5, 35.8, 36.2, 36.1, 36.4]): + _seed_thermal_session(db, now - (10 - i) * 7200, ambient_c=25.0, rise_ref_c=rise, + cooldown_s=900.0, ambient_end_c=25.0) + for i, rise in enumerate([43.4, 43.7, 43.5]): + _seed_thermal_session(db, now - (4 - i) * 7200, ambient_c=25.0, rise_ref_c=rise, + amps=39.6, cooldown_s=900.0, ambient_end_c=25.0) + fits = thermal.fit_sessions(db, now) + drift = thermal.detect_drift(fits) + assert drift is not None and drift["n"] == 9 + # The raw split the old watch used would have called this an alert. + naive = median(fit["rise_ref_c"] for fit in fits[-3:]) - median(fit["rise_ref_c"] for fit in fits[:-3]) + assert naive > thermal.DRIFT_WARN_C + assert drift["drifting"] is False + assert "current" in drift["covariates"] + # A cap applied once and kept is nearly collinear with the calendar, so + # the slope cannot be pinned down — and the widened interval says so + # instead of the verdict quietly picking one explanation. + assert "current" in drift["collinear_with_time"] + assert drift["confident"] is False + + # Vary the current instead of stepping it, and the confound separates: + # the same nine fits, interleaved, pin the slope near zero. + now2 = now + 400 * 86400 + for i in range(9): + amps, rise = (39.6, 43.5) if i % 2 else (48.6, 36.2) + _seed_thermal_session(db, now2 - (10 - i) * 7200, ambient_c=25.0, rise_ref_c=rise, + amps=amps, cooldown_s=900.0, ambient_end_c=25.0) + mixed = thermal.detect_drift([fit for fit in thermal.fit_sessions(db, now2 + 3600) + if fit["start_ts"] > now2 - 20 * 7200]) + assert mixed is not None and mixed["drifting"] is False + assert not mixed["collinear_with_time"] + assert abs(mixed["delta_c"]) < thermal.DRIFT_WARN_C + + +async def test_thermal_drift_reports_ambient_confound(db): + # The premise of rise_ref is that subtracting ambient leaves a number + # about the connector. On an install where it does not — the fits still + # slope against garage temperature — a cooling autumn walks the rise + # upward all by itself. The regression must attribute that to ambient + # rather than to time, and must say the measurement is compromised + # rather than presenting an adjusted number as clean. + now = time.time() + ambients = [34.0, 32.5, 31.0, 29.5, 28.0, 26.5, 25.0, 24.0] + for i, ambient in enumerate(ambients): + _seed_thermal_session(db, now - (12 - i) * 7200, ambient_c=ambient, + rise_ref_c=36.0 + (29.0 - ambient) * 0.8, + cooldown_s=900.0, ambient_end_c=ambient) + fits = thermal.fit_sessions(db, now) + drift = thermal.detect_drift(fits) + assert drift is not None + naive = median(fit["rise_ref_c"] for fit in fits[-3:]) - median(fit["rise_ref_c"] for fit in fits[:-3]) + assert naive > thermal.DRIFT_WARN_C # the median split would have alerted + assert "ambient" in drift["covariates"] + assert drift["ambient_coef_c_per_c"] < -0.5 + assert "ambient" in drift["compromised_by"] + assert drift["drifting"] is False and abs(drift["delta_c"]) < thermal.DRIFT_WARN_C + + +async def test_thermal_drift_ambient_confound_caps_verdict_at_lead(db): + # Same broken premise, but this time with a real trend on top of it. The + # delta is material and the interval clears zero, yet the number rests on + # fits that are measuring something other than connector resistance — so + # it is a lead to be explained, never an alert to act on. + now = time.time() + ambients = [34.0, 32.5, 31.0, 29.5, 28.0, 26.5, 25.0, 24.0] + for i, ambient in enumerate(ambients): + _seed_thermal_session(db, now - (12 - i) * 7200, ambient_c=ambient, + rise_ref_c=36.0 + (29.0 - ambient) * 0.8 + i * 1.2, + cooldown_s=900.0, ambient_end_c=ambient) + drift = thermal.detect_drift(thermal.fit_sessions(db, now)) + assert drift is not None + assert drift["delta_c"] > thermal.DRIFT_WARN_C and drift["confident"] is True + assert drift["compromised_by"] == ["ambient"] + assert drift["drifting"] is False and drift["lead"] is True + + async def test_thermal_baseline_anchor(db): now = time.time() rises = [36.0, 36.5, 35.8, 36.2, 42.0, 41.5, 42.3] diff --git a/wallmonitor/poller.py b/wallmonitor/poller.py index cde9dce..593cdf9 100644 --- a/wallmonitor/poller.py +++ b/wallmonitor/poller.py @@ -527,9 +527,10 @@ async def recheck_thermal_drift(self, ts: float) -> None: return ci_lo, ci_hi = drift["delta_ci95_c"] body = ( - f"Recent sessions run +{drift['recent_rise_c']:.1f} °C vs a +{drift['baseline_rise_c']:.1f} °C " - f"baseline at the same current (Δ {drift['delta_c']:.1f} °C, 95% CI {ci_lo:.1f}..{ci_hi:.1f}, " - f"n={drift['baseline_n']}+{drift['recent_n']}" + f"Fitted rise is trending +{drift['slope_c_per_day'] * 30.0:.1f} °C/month with ambient and " + f"charge current held: +{drift['baseline_rise_c']:.1f} → +{drift['recent_rise_c']:.1f} °C across " + f"{drift['span_days']:.0f} days (Δ {drift['delta_c']:.1f} °C, 95% CI {ci_lo:.1f}..{ci_hi:.1f}, " + f"n={drift['n']} free-running fits" ) if drift["drifting"]: # Confirmed: the interval clears zero and the delta is material. @@ -550,16 +551,31 @@ async def recheck_thermal_drift(self, ts: float) -> None: if cleared: await self._event(ts, "thermal_drift_cleared", drift) if drift["lead"]: - # Past the floor but inside this install's own scatter: a lead - # for the dashboard and a quiet push, once per episode — no alert. + # Past the floor but not confirmable: either inside this install's + # own scatter, or resting on a measurement the install itself + # undermines. A lead for the dashboard and a quiet push, once per + # episode — no alert. Say which of the two it is, because they + # call for different things: more sessions settle scatter, while a + # rise that still tracks garage temperature needs the confounder + # found before any number here means much. + if "ambient" in drift["compromised_by"]: + why = ( + f"; but this install's rise still moves {drift['ambient_coef_c_per_c']:+.2f} °C per °C of " + "ambient after the subtraction, so the fits are carrying something other than connector " + "resistance) — find that before reading the trend as hardware" + ) + else: + why = ( + f"; needs Δ ≥ {drift['threshold_c']:.1f} °C at this install's scatter to confirm) " + "— worth a look at the handle and pins next time you're there; more sessions will settle it" + ) if not self._drift_lead_active: self._drift_lead_active = True await self._event(ts, "thermal_drift_lead", drift) await self._notify( "thermal_drift_lead", "Heat rise may be climbing — a lead, not yet confirmed", - body + f"; needs Δ ≥ {drift['threshold_c']:.1f} °C at this install's scatter to confirm) " - "— worth a look at the handle and pins next time you're there; more sessions will settle it.", + body + why + ".", drift, ) else: diff --git a/wallmonitor/static/app.js b/wallmonitor/static/app.js index a8d4ee2..9212a9c 100644 --- a/wallmonitor/static/app.js +++ b/wallmonitor/static/app.js @@ -1086,19 +1086,29 @@ async function viewLive(root) { // same current means added resistance somewhere in the current path. const drift = data.drift; let driftLine = null; + // The trend is a modelled slope with ambient and charge current held, + // not a difference of medians — so the number quoted is what the rise + // did over the window at fixed conditions, and the covariate that would + // otherwise have explained it is named when it is doing real work. + const ambientConfounded = drift && (drift.compromised_by || []).includes("ambient"); if (drift && drift.lead) { driftLine = el("div", { class: "note" }, chipFor("warning", "heat rise: lead"), - ` Recent sessions average +${fmtNum(drift.recent_rise_c, 1)} °C at ${fmtNum(model.ref_current_a, 0)} A vs a ` + - `+${fmtNum(drift.baseline_rise_c, 1)} °C baseline (Δ ${fmtNum(drift.delta_c, 1)} °C) — past the ` + - `${fmtNum(drift.floor_c, 1)} °C floor but within this install's session-to-session scatter, which needs ` + - `Δ ≥ ${fmtNum(drift.threshold_c, 1)} °C to confirm. Not an alert; worth a look at the handle and pins.`); + ` Fitted rise trends +${fmtNum(drift.slope_c_per_day * 30, 1)} °C/month at fixed ambient and current ` + + `(+${fmtNum(drift.baseline_rise_c, 1)} → +${fmtNum(drift.recent_rise_c, 1)} °C over ${fmtNum(drift.span_days, 0)} days, ` + + `Δ ${fmtNum(drift.delta_c, 1)} °C) — past the ${fmtNum(drift.floor_c, 1)} °C floor but ` + + (ambientConfounded + ? `this install's rise still moves ${fmtNum(drift.ambient_coef_c_per_c, 2)} °C per °C of ambient after the ` + + "subtraction, so the fits are measuring something besides connector resistance. Not an alert; find that first." + : `within this install's session-to-session scatter, which needs Δ ≥ ${fmtNum(drift.threshold_c, 1)} °C to ` + + "confirm. Not an alert; worth a look at the handle and pins.")); } else if (drift && drift.drifting) { driftLine = el("div", { class: "note" }, chipFor("serious", "heat rise increasing"), - ` Recent sessions average +${fmtNum(drift.recent_rise_c, 1)} °C at ${fmtNum(model.ref_current_a, 0)} A vs a ` + - `+${fmtNum(drift.baseline_rise_c, 1)} °C baseline (Δ ${fmtNum(drift.delta_c, 1)} °C). More heat at the same current ` + - `means added resistance — inspect the handle and charge-port pins, and have the terminal torque checked.` + + ` Fitted rise trends +${fmtNum(drift.slope_c_per_day * 30, 1)} °C/month at fixed ambient and current ` + + `(+${fmtNum(drift.baseline_rise_c, 1)} → +${fmtNum(drift.recent_rise_c, 1)} °C over ${fmtNum(drift.span_days, 0)} days, ` + + `Δ ${fmtNum(drift.delta_c, 1)} °C). More heat at the same current and ambient means added resistance — ` + + "inspect the handle and charge-port pins, and have the terminal torque checked." + (drift.off_current_n ? ` (${drift.off_current_n} session${drift.off_current_n === 1 ? "" : "s"} away from the usual ` + `~${fmtNum(drift.typical_current_a, 0)} A excluded from the comparison.)` : "")); } @@ -1132,7 +1142,7 @@ async function viewLive(root) { const modelNote = `Model: τ ≈ ${fmtNum(model.tau_min, 1)} min, +${fmtNum(model.rise_ref_c, 0)} °C at ${fmtNum(model.ref_current_a, 0)} A — ` + (model.fitted ? `fitted from ${model.tau_fits} recorded session ramp${model.tau_fits === 1 ? "" : "s"}.` + priorNote : "defaults from one verified install, used until this charger has fits of its own; refits automatically as sessions accumulate.") + - (drift && !drift.drifting && !drift.lead ? ` Heat rise stable across the last ${drift.recent_n + drift.baseline_n} fitted sessions` + + (drift && !drift.drifting && !drift.lead ? ` Heat rise stable across ${drift.n} free-running fitted sessions` + `${drift.off_current_n ? ` (${drift.off_current_n} off-current session${drift.off_current_n === 1 ? "" : "s"} excluded)` : ""}.` : "") + idleNote; thermalCard.append(el("div", { class: "chart-card" }, @@ -1746,11 +1756,11 @@ async function viewAlerts(root, rangeKey = "7d") { unit: "°C", digits: 1, height: 180, vlines: calMarks, }); if (drift) { - // The verdict carries its own uncertainty — a delta from a handful of + // The verdict carries its own uncertainty — a slope from a handful of // fits is a lead, not a conviction, and the note must show which. const [ciLo, ciHi] = drift.delta_ci95_c || [null, null]; const sureness = ciLo == null ? "" : - ` · 95% CI ${fmtNum(ciLo, 1)}..${fmtNum(ciHi, 1)} °C from n=${drift.baseline_n}+${drift.recent_n}`; + ` · 95% CI ${fmtNum(ciLo, 1)}..${fmtNum(ciHi, 1)} °C from n=${drift.n}`; const pooled = drift.cross_current_n ? ` (${drift.cross_current_n} ambient-bracketed fit${drift.cross_current_n > 1 ? "s" : ""} pooled from other charge currents)` : ""; @@ -1761,7 +1771,8 @@ async function viewAlerts(root, rangeKey = "7d") { (drift.threshold_c > drift.floor_c + 0.05 ? ` (the ${fmtNum(drift.floor_c, 1)} °C floor, raised to what this install's scatter needs to confirm)` : ` (the ${fmtNum(drift.floor_c, 1)} °C floor)`); - const summary = `recent median +${fmtNum(drift.recent_rise_c, 1)} °C vs baseline +${fmtNum(drift.baseline_rise_c, 1)} °C`; + const summary = `+${fmtNum(drift.baseline_rise_c, 1)} → +${fmtNum(drift.recent_rise_c, 1)} °C over ` + + `${fmtNum(drift.span_days, 0)} days (${fmtNum(drift.slope_c_per_day * 30, 2)} °C/month)`; rise.card.append(el("div", { class: "note" }, (drift.drifting ? `Confirmed: ${summary} (Δ ${fmtNum(drift.delta_c, 1)} °C; ${thresholdNote}) — a monitor alert is active` @@ -1769,6 +1780,43 @@ async function viewAlerts(root, rangeKey = "7d") { ? `Lead: ${summary} (Δ ${fmtNum(drift.delta_c, 1)} °C, past the floor but not yet confirmed; ${thresholdNote}) — ` + "no alert; more sessions will settle it" : `Stable: ${summary} (${thresholdNote})`) + sureness + pooled + ".")); + // What the trend was actually held against, and what it could not be + // held against. A slope is only as meaningful as its controls, so the + // covariates are named rather than left implicit — and an install + // whose rise still tracks ambient is told so plainly, because that is + // the finding, not a footnote to it. + const covs = drift.covariates || []; + rise.card.append(el("div", { class: "note" }, + `Estimated by regressing fitted rise on time` + + (covs.length ? ` while holding ${covs.join(" and ")}` : ", with no covariate varying enough to hold") + + `; residual scatter ±${fmtNum(drift.resid_sd_c, 1)} °C on ${drift.dof} degrees of freedom` + + (drift.current_coef_c_per_a != null + ? `. Charge current carries ${fmtNum(drift.current_coef_c_per_a, 2)} °C per amp here — the residual ` + + "error in the I² normalization, absorbed rather than left to masquerade as a trend" + : "") + ".")); + if (drift.ambient_coef_c_per_c != null) { + const confounded = (drift.compromised_by || []).includes("ambient"); + rise.card.append(el("div", { class: "note" }, + `Rise vs ambient: ${fmtNum(drift.ambient_coef_c_per_c, 2)} ± ${fmtNum(drift.ambient_coef_se, 2)} °C per °C. ` + + (confounded + ? "Subtracting ambient was supposed to leave a number that depends on the connector and not the weather; " + + "here it did not, so something the model does not carry — multi-day heat soak, a charger regulating to a " + + "fixed handle temperature, a badly sited sensor — is inside the measurement. The trend above is adjusted " + + "for it, but this install can raise a lead and never an alert until it is found." + : "Flat enough that the ambient subtraction is doing its job."))); + } + // Regulated windows are the fits the charger wrote itself: current + // trimmed back as the handle warmed, so the plateau is a setpoint and + // its "rise" tracks how close the handle got to the limit. Excluding + // them is what makes the rest comparable; saying how many were + // excluded is what keeps a thin verdict from looking well-fed. + if (drift.regulated_n) { + rise.card.append(el("div", { class: "note" }, + `${drift.regulated_n} further fit${drift.regulated_n === 1 ? " was" : "s were"} excluded: the charger trimmed ` + + "charge current back inside the ramp window, so the plateau it reached was one the charger held, not the " + + "connector's own. A charge at a current low enough to run unregulated end to end is worth more to this watch " + + "than several at full rate.")); + } } // Ambient bracketing: fits that read ambient at both ends of the load // window are de-trended for weather that moved during the charge — the diff --git a/wallmonitor/thermal.py b/wallmonitor/thermal.py index d327419..5ad0498 100644 --- a/wallmonitor/thermal.py +++ b/wallmonitor/thermal.py @@ -208,6 +208,31 @@ def ambient_from_idle_handle(handle_c: float, model: IdleOffset = BUILTIN_IDLE_O PREFIX_SPAN_TAU = 2.5 PREFIX_SPAN_MIN_S = 1800.0 +# The steady-prefix band (10% of the reference current) is wide enough to +# hide the charger's own thermal regulation: a Gen 3 that trims 48.6 A to +# 44.7 A as the handle nears its limit never leaves the band, so the ramp +# keeps collecting samples whose flattening is *caused by the current +# dropping*. The exponential then reads that as the plateau — a lower rise +# paired with a faster tau, passing every gate with a fine RMSE. Measured on +# one install: fits whose current sagged read a median 33.3 C rise against +# 37.2 C for the same charger's steady windows, and because foldback starts +# sooner in a hot garage the bias tracked ambient, manufacturing a 7 C +# "drift" verdict out of the seasons. +# +# So each fit records whether its window was *free-running*: current held +# flat end to end, making the fitted plateau the connector's own equilibrium +# rather than one the charger imposed. Only free-running fits are compared +# by the degradation watch. The other half of "the plateau was real" — +# whether the window ran long enough to observe it — is MIN_SPAN_TAU above, +# already enforced before any fit is emitted. +# +# The threshold separates the two populations with room to spare: on that +# install steady windows sagged <= 0.6% while regulated ones sagged >= 3.9%. +# The absolute floor keeps sensor quantization on a low-current charge from +# reading as regulation. +FREE_CURRENT_SAG_FRAC = 0.015 +FREE_CURRENT_SAG_MIN_A = 0.5 + # Live-forecast gate: a steady-current window must hold this many samples # over this much time before its trajectory is projected. TRAJECTORY_MIN_SAMPLES = 8 @@ -352,6 +377,34 @@ def _steady_current_prefix(samples: list[dict], max_span_s: float = PREFIX_SPAN_ return prefix +def _current_sag_a(prefix: list[dict]) -> float: + """How far the window's current fell from its opening to its close. + + Head and tail quarters are compared by median, so a single dropped + sample or a momentary blip cannot pass for regulation. Signed: a + negative sag means the current *rose* across the window, which breaks + the constant-current premise just as thoroughly. + """ + quarter = max(2, len(prefix) // 4) + head = median(sample["vehicle_current_a"] for sample in prefix[:quarter]) + tail = median(sample["vehicle_current_a"] for sample in prefix[-quarter:]) + return head - tail + + +def _free_plateau(prefix: list[dict], sag_a: float) -> bool: + """Did the charger leave this window alone? + + True when the current held flat across the whole window, so the fitted + plateau is the connector's own equilibrium at that current. False when + the charger was trimming current back — the plateau is then a setpoint + the charger held, and the fitted rise says more about how close the + handle got to the limit than about connector resistance. + """ + i_ref = median(sample["vehicle_current_a"] for sample in prefix) + tolerance = max(FREE_CURRENT_SAG_MIN_A, FREE_CURRENT_SAG_FRAC * i_ref) + return abs(sag_a) <= tolerance + + def _segments(rows: list[dict]) -> list[tuple[float, float]]: """(start, end) timestamps of distinct charging segments in a session. @@ -569,7 +622,14 @@ def fit_sessions(db: Database, now: float, lookback_days: float = 120.0) -> list ambient_source ("measured" for a stationary LAN sensor, "measured_car" for a parked vehicle's sensor, else "pre_idle" or "cooldown_tail"), ambient_c, and — when bracketed — ambient_end_c and ambient_drift_c - (end minus start).""" + (end minus start). + + Every fit also carries current_sag_a (how far current fell across the + window) and free_plateau: False means the charger was trimming current + back as the handle warmed, so the fitted plateau is a setpoint it held + rather than the connector's own equilibrium. Such fits still serve the + forecast — they describe what the handle actually did — but the + degradation watch compares only free-running ones.""" sessions = [ session for session in db.sessions_range(now - lookback_days * 86400, now) @@ -670,6 +730,7 @@ def fit_sessions(db: Database, now: float, lookback_days: float = 120.0) -> list rise = (t_inf - ambient) * (REF_CURRENT_A / i_med) ** 2 if not (RISE_RANGE_C[0] <= rise <= RISE_RANGE_C[1]): rise = None + sag_a = _current_sag_a(prefix) fits.append( { "session_id": sess["id"], @@ -678,6 +739,8 @@ def fit_sessions(db: Database, now: float, lookback_days: float = 120.0) -> list "rise_ref_c": round(rise, 2) if rise is not None else None, "rmse_c": round(rmse, 3), "current_a": round(i_med, 1), + "current_sag_a": round(sag_a, 2), + "free_plateau": _free_plateau(prefix, sag_a), "ambient_source": ambient_source if rise is not None else None, "ambient_c": round(ambient, 2) if rise is not None else None, "ambient_end_c": ( @@ -717,39 +780,139 @@ def fit_history(db: Database, now: float, lookback_days: float = 120.0, # A loose lug or degrading contact shows up as extra resistance: more heat # rise for the same current. Prediction alone hides that (the rolling median -# just follows it), so the drift watch compares recent sessions against the -# earlier baseline and flags a sustained increase. -DRIFT_RECENT_N = 3 -DRIFT_MIN_BASELINE_N = 3 +# just follows it), so the drift watch models rise against time and flags a +# sustained increase. +# +# It models rather than compares, because a recent-vs-baseline median split +# answers the wrong question. The split asks "are the last few fits higher?", +# which any covariate that moved with the calendar answers for it: a garage +# that cooled between the two halves, or a vehicle capped to a lower current +# whose (48/I)^2 normalization then lifts every recent fit. On one install +# that split reported +7.2 C with a 95% CI of [5.4, 9.1] — "statistically +# confirmed" — from a connector whose rise, regressed on time with ambient +# and current held, was moving +0.01 C/day, indistinguishable from flat. The +# confidence was real; it was confidence in the wrong estimand. +# +# So the watch fits rise_ref ~ days + ambient + current over the whole +# comparable history and reads the *days* coefficient. Ambient and current +# stop being confounders and become covariates, the estimate uses every fit +# instead of six, and when a covariate genuinely cannot be separated from +# time the collinearity inflates the slope's standard error and the verdict +# declines to confirm — which is the honest outcome, reached automatically. +DRIFT_MIN_N = 6 +DRIFT_TYPICAL_N = 3 # newest fits defining the install's current operating point DRIFT_WARN_C = 2.5 # materiality floor, not the trigger — see detect_drift DRIFT_ALERT = "Handle heat rise increasing (check connector/wiring)" # Cross-current pooling: fits whose ambient was bracketed at both ends are # trustworthy enough under the I^2 normalization to join the comparison from # a wider current band; start-only fits must still match the typical current. +# The regression carries a current term of its own, so pooled fits are +# adjusted rather than merely admitted — a residual error in the I^2 +# normalization lands on that coefficient instead of on the time slope. DRIFT_POOL_BAND_FRAC = 0.25 +# A covariate earns a column only when the history actually moved in it. +# Regressing on a covariate that barely varies buys nothing and spends a +# degree of freedom the small-sample t-multiplier charges dearly for. +DRIFT_AMBIENT_SPREAD_C = 3.0 +DRIFT_CURRENT_SPREAD_A = 2.0 + +# Reported, not gated: how strongly a covariate moved with the calendar. +# Past this the two cannot be told apart, and the slope's standard error +# will already be showing it — the flag exists so the UI can say *why* an +# apparently large delta refused to confirm. +DRIFT_COLLINEAR_R = 0.8 + +# The premise of rise_ref is that subtracting ambient leaves a number that +# depends on the connector and not on the weather. An install where rise +# still moves this much per degree of ambient — materially, and resolved +# well enough to be sure of the sign — has broken that premise: something +# the model does not carry (multi-day heat soak, a charger regulating to a +# fixed handle temperature, a badly sited sensor) is in the measurement. +# The regression adjusts for it, but adjustment is not understanding, so +# such an install can still raise a lead and never an alert. +DRIFT_AMBIENT_CONFOUND_C = 0.3 + # The settings key holding the baseline anchor: a timestamp before which # fits are excluded from the drift comparison. Set it when the hardware has # been inspected and verified (or fixed) — from then on "baseline" means # "verified healthy", not "the first charges the monitor happened to see". BASELINE_ANCHOR_KEY = "thermal_baseline_anchor_ts" -# Two-sided 95% Student-t multipliers by degrees of freedom, for the delta's -# confidence interval. Small-sample medians are noisy; a plain 1.96 would -# overstate the confidence exactly when the history is thinnest. -_T95 = {1: 12.71, 2: 4.30, 3: 3.18, 4: 2.78, 5: 2.57, 6: 2.45, 7: 2.36, 8: 2.31, 9: 2.26, 10: 2.23} +# Two-sided 95% Student-t multipliers by residual degrees of freedom, for +# the slope's confidence interval. A plain 1.96 would overstate the +# confidence exactly when the history is thinnest. +_T95 = {1: 12.71, 2: 4.30, 3: 3.18, 4: 2.78, 5: 2.57, 6: 2.45, 7: 2.36, 8: 2.31, 9: 2.26, 10: 2.23, + 11: 2.20, 12: 2.18, 13: 2.16, 14: 2.14, 15: 2.13, 16: 2.12, 17: 2.11, 18: 2.10, 19: 2.09, + 20: 2.09, 25: 2.06, 30: 2.04, 40: 2.02, 60: 2.00} + + +def _t95(dof: int) -> float: + """Two-sided 95% t multiplier. An untabulated dof falls back to the + next-lower tabulated one, which is the larger multiplier — rounding + toward caution rather than away from it.""" + if dof < 1: + return _T95[1] + if dof in _T95: + return _T95[dof] + return _T95[max(key for key in _T95 if key <= dof)] + + +def _invert(matrix: list[list[float]]) -> list[list[float]] | None: + """Gauss-Jordan inverse with partial pivoting; None if singular. + + The pivot tolerance is relative to the largest entry, because the design + matrix mixes columns of wildly different scale (a count of fits against + a sum of squared day-offsets).""" + size = len(matrix) + scale = max((abs(value) for row in matrix for value in row), default=0.0) + if scale <= 0.0: + return None + aug = [row[:] + [1.0 if i == j else 0.0 for j in range(size)] for i, row in enumerate(matrix)] + for col in range(size): + pivot = max(range(col, size), key=lambda row: abs(aug[row][col])) + if abs(aug[pivot][col]) < 1e-10 * scale: + return None + aug[col], aug[pivot] = aug[pivot], aug[col] + divisor = aug[col][col] + aug[col] = [value / divisor for value in aug[col]] + for row in range(size): + if row != col and aug[row][col] != 0.0: + factor = aug[row][col] + aug[row] = [value - factor * base for value, base in zip(aug[row], aug[col])] + return [row[size:] for row in aug] + + +def _ols(y: list[float], design: list[list[float]]) -> tuple[list[float], list[float], float, int] | None: + """Ordinary least squares: (coefficients, standard errors, residual sd, + residual dof), or None when the design is singular or leaves too few + degrees of freedom for the standard errors to mean anything.""" + rows, cols = len(y), len(design[0]) + dof = rows - cols + if dof < 2: + return None + xtx = [[sum(design[i][a] * design[i][b] for i in range(rows)) for b in range(cols)] + for a in range(cols)] + inv = _invert(xtx) + if inv is None: + return None + xty = [sum(design[i][a] * y[i] for i in range(rows)) for a in range(cols)] + beta = [sum(inv[a][b] * xty[b] for b in range(cols)) for a in range(cols)] + sse = sum((y[i] - sum(design[i][j] * beta[j] for j in range(cols))) ** 2 for i in range(rows)) + variance = sse / dof + se = [math.sqrt(max(variance * inv[j][j], 0.0)) for j in range(cols)] + return beta, se, math.sqrt(variance), dof -def _median_stats(values: list[float]) -> tuple[float, float, float]: - """(median, MAD, standard error of the median). +def _pearson(xs: list[float], ys: list[float]) -> float: + """Correlation coefficient; 0.0 when either side is constant.""" + x_bar, y_bar = sum(xs) / len(xs), sum(ys) / len(ys) + sxy = sum((x - x_bar) * (y - y_bar) for x, y in zip(xs, ys)) + sxx = sum((x - x_bar) ** 2 for x in xs) + syy = sum((y - y_bar) ** 2 for y in ys) + return sxy / math.sqrt(sxx * syy) if sxx > 0 and syy > 0 else 0.0 - Spread comes from the median absolute deviation — robust to the odd wild - fit — scaled to a normal-equivalent sigma (1.4826) and to the median's - sampling efficiency (1.2533 / sqrt(n)).""" - med = median(values) - mad = median(abs(value - med) for value in values) - return med, mad, 1.2533 * (1.4826 * mad) / math.sqrt(len(values)) # Actionable warning: the live forecast puts the 65 C trip inside this # horizon, so the user still has time to lower the vehicle's charge current @@ -759,50 +922,62 @@ def _median_stats(values: list[float]) -> tuple[float, float, float]: def detect_drift(fits: list[dict], anchor_ts: float | None = None) -> dict | None: - """Compare the last few sessions' fitted rise against the baseline. - - Returns None while there is too little history to judge; otherwise a - verdict dict with the medians compared. Only rise (not tau) is watched: - added contact resistance changes how much heat is made, not how fast the - handle mass warms. - - Only sessions charging near the install's typical current are compared. - rise_ref_c is normalized by (REF/I)^2, and far from the measured current - that normalization amplifies ordinary fit error — a session at 40 A with - an unremarkable raw rise extrapolates to an alarming number at 48 A, and - with only DRIFT_RECENT_N recent sessions a single such point can swing - the median past the threshold and manufacture a drift verdict. Fits with - bracketed ambient are cleaner, so they join from a wider current band — - that is what lets a baseline recorded at 48 A keep judging charges after - the vehicle is capped to 40 A, instead of the verdict going dark. - - "Typical" is the median current of the newest fits, not of all history: - when the user caps the vehicle at a new current, the watch follows the - new operating point. + """Regress fitted rise on time, holding ambient and current, and read the + time coefficient. + + Returns None while there is too little comparable history to judge; + otherwise a verdict dict. Only rise (not tau) is watched: added contact + resistance changes how much heat is made, not how fast the handle mass + warms. + + **Only free-running fits are compared.** A window whose current the + charger was trimming back as the handle warmed has a plateau the charger + chose, and its fitted rise moves with how close the handle got to the + limit — which is to say, with ambient. Those fits still describe what the + handle did, so the forecast keeps them; a degradation comparison cannot + use them at all. + + **Only sessions near the install's recent operating current**, with + ambient-bracketed fits pooled in from a wider band. "Typical" is the + median current of the newest fits, not of all history: when the user caps + the vehicle at a new current, the watch follows the new operating point. anchor_ts, when set, excludes fits from before it: the user has had the hardware inspected and verified, so "baseline" means "verified healthy" from that moment, not "the first charges the monitor happened to see". - The verdict carries its own uncertainty: MAD spread per side and a - Student-t ~95% confidence interval on the delta. "drifting" — the alert - — needs both: the interval clears zero ("confident") *and* the delta is - material (>= DRIFT_WARN_C, the floor below which a real increase isn't - worth an inspection). A delta past the floor whose interval still - straddles zero is a "lead": shown, pushed quietly, but no alert. The - effective threshold ("threshold_c") is therefore the larger of the floor - and what this install's own scatter requires, so a noisy install must - show more before the watch alarms and a quiet one less — the fixed - 2.5 °C tripwire sat near one sigma on a real install and fired on - scatter. + The estimate is the modelled change across the observed span — + slope x days — not a difference of medians, so ambient and charge current + are adjusted for rather than assumed away. Its confidence interval is the + slope's, at a small-sample Student-t multiplier. "drifting" — the alert — + needs the interval to clear zero ("confident"), the change to be material + (>= DRIFT_WARN_C, the floor below which a real increase is not worth an + inspection), and the measurement itself to be sound: an install whose + rise still tracks ambient after the subtraction (see + DRIFT_AMBIENT_CONFOUND_C) is measuring something other than connector + resistance, and can raise a lead but never an alert. A change past the + floor that fails either test is a "lead": shown, pushed quietly, no alert + row. The effective threshold ("threshold_c") is the larger of the floor + and what this install's own scatter requires. """ - usable = [fit for fit in fits if fit["rise_ref_c"] is not None] - if anchor_ts is not None: - usable = [fit for fit in usable if fit["start_ts"] >= anchor_ts] - if len(usable) < DRIFT_RECENT_N + DRIFT_MIN_BASELINE_N: + usable = [ + fit for fit in fits + if fit["rise_ref_c"] is not None + # Fits predating the free-plateau gate carry no verdict either way; + # nothing better to assume than that the window was steady. + and fit.get("free_plateau", True) + and (anchor_ts is None or fit["start_ts"] >= anchor_ts) + ] + regulated_n = sum( + 1 for fit in fits + if fit["rise_ref_c"] is not None + and not fit.get("free_plateau", True) + and (anchor_ts is None or fit["start_ts"] >= anchor_ts) + ) + if len(usable) < DRIFT_MIN_N: return None usable.sort(key=lambda fit: fit["start_ts"]) - typical_a = median(fit["current_a"] for fit in usable[-DRIFT_RECENT_N:]) + typical_a = median(fit["current_a"] for fit in usable[-DRIFT_TYPICAL_N:]) band = max(2.0, 0.1 * typical_a) pool_band = DRIFT_POOL_BAND_FRAC * typical_a comparable = [ @@ -814,34 +989,100 @@ def detect_drift(fits: list[dict], anchor_ts: float | None = None) -> dict | Non and abs(fit["current_a"] - typical_a) <= pool_band ) ] - rises = [(fit["start_ts"], fit["rise_ref_c"]) for fit in comparable] - rises.sort(key=lambda entry: entry[0]) - if len(rises) < DRIFT_RECENT_N + DRIFT_MIN_BASELINE_N: + if len(comparable) < DRIFT_MIN_N: + return None + + origin = comparable[0]["start_ts"] + days = [(fit["start_ts"] - origin) / 86400.0 for fit in comparable] + span_days = days[-1] - days[0] + if span_days <= 0.0: + return None + rises = [fit["rise_ref_c"] for fit in comparable] + currents = [fit["current_a"] for fit in comparable] + # Ambient is missing on fits that read it from neither end; centre what + # is there on its own mean so the intercept keeps its meaning. + ambients = [fit.get("ambient_c") for fit in comparable] + have_ambient = all(value is not None for value in ambients) + + # Each covariate column is earned by variation. Centring them makes the + # intercept the predicted rise at the first fit under average conditions, + # which is what the UI reports as the baseline. + columns: list[tuple[str, list[float]]] = [] + if have_ambient and max(ambients) - min(ambients) >= DRIFT_AMBIENT_SPREAD_C: + mean_ambient = sum(ambients) / len(ambients) + columns.append(("ambient", [value - mean_ambient for value in ambients])) + if max(currents) - min(currents) >= DRIFT_CURRENT_SPREAD_A: + mean_current = sum(currents) / len(currents) + columns.append(("current", [value - mean_current for value in currents])) + + # A thin history cannot afford every column. Drop them back to front — + # current first, ambient last — until the design is estimable, because + # ambient is the confounder this watch exists to survive and the one + # most likely to move with the calendar on its own. + fit_result = None + while True: + design = [[1.0, day] + [column[i] for _, column in columns] for i, day in enumerate(days)] + fit_result = _ols(rises, design) + if fit_result is not None or not columns: + break + columns.pop() + if fit_result is None: return None - recent = [rise for _, rise in rises[-DRIFT_RECENT_N:]] - baseline = [rise for _, rise in rises[:-DRIFT_RECENT_N]] - recent_med, recent_mad, recent_se = _median_stats(recent) - baseline_med, baseline_mad, baseline_se = _median_stats(baseline) - delta = recent_med - baseline_med - delta_se = math.sqrt(recent_se**2 + baseline_se**2) - t_mult = _T95.get(len(recent) + len(baseline) - 2, 2.0) + beta, se, resid_sd, dof = fit_result + names = [name for name, _ in columns] + coef = {name: (beta[2 + i], se[2 + i]) for i, name in enumerate(names)} + + slope, slope_se = beta[1], se[1] + delta = slope * span_days + delta_se = slope_se * span_days + t_mult = _t95(dof) ci_lo, ci_hi = delta - t_mult * delta_se, delta + t_mult * delta_se confident = ci_lo > 0.0 threshold = max(DRIFT_WARN_C, t_mult * delta_se) - drifting = confident and delta >= DRIFT_WARN_C + + # What could not be told apart from the calendar, and what the ambient + # subtraction failed to remove. Neither is an error; both are reasons a + # delta that looks large is not yet an alert. + compromised: list[str] = [] + collinear = { + name: _pearson(days, column) + for name, column in columns + if abs(_pearson(days, column)) > DRIFT_COLLINEAR_R + } + ambient_coef, ambient_se = coef.get("ambient", (None, None)) + ambient_confounded = ( + ambient_coef is not None + and abs(ambient_coef) >= DRIFT_AMBIENT_CONFOUND_C + and abs(ambient_coef) > 2.0 * ambient_se + ) + if ambient_confounded: + compromised.append("ambient") + + drifting = confident and delta >= DRIFT_WARN_C and not compromised cross_current = sum(1 for fit in comparable if abs(fit["current_a"] - typical_a) > band) return { "drifting": drifting, "lead": delta >= DRIFT_WARN_C and not drifting, "confident": confident, - "recent_rise_c": round(recent_med, 2), - "baseline_rise_c": round(baseline_med, 2), + "baseline_rise_c": round(beta[0], 2), + "recent_rise_c": round(beta[0] + slope * span_days, 2), "delta_c": round(delta, 2), "delta_ci95_c": [round(ci_lo, 2), round(ci_hi, 2)], - "recent_mad_c": round(recent_mad, 2), - "baseline_mad_c": round(baseline_mad, 2), - "recent_n": len(recent), - "baseline_n": len(baseline), + "slope_c_per_day": round(slope, 4), + "slope_se_c_per_day": round(slope_se, 4), + "span_days": round(span_days, 1), + "n": len(comparable), + "resid_sd_c": round(resid_sd, 2), + "dof": dof, + "covariates": names, + "ambient_coef_c_per_c": round(ambient_coef, 3) if ambient_coef is not None else None, + "ambient_coef_se": round(ambient_se, 3) if ambient_se is not None else None, + "current_coef_c_per_a": ( + round(coef["current"][0], 3) if "current" in coef else None + ), + "collinear_with_time": {name: round(value, 2) for name, value in collinear.items()}, + "compromised_by": compromised, + "regulated_n": regulated_n, "typical_current_a": round(typical_a, 1), "off_current_n": len(usable) - len(comparable), "cross_current_n": cross_current,