Reject outlier readings in the battery estimate

A single noisy ADC/regulator glitch (see firmware/main/battery.c)
survives the existing recharge filter: whichever way it reads, one of
the two steps around it (into a dip, or out of a spike) still looks
like an ordinary drop and got averaged straight into the remaining-
time estimate, letting one bad reading swing it dramatically.

Added a MAD-based modified z-score outlier check on top of the
existing recency-weighted average. Had to special-case the standard
MAD degenerating to exactly 0, which happens whenever more than half
the steps share the same value -- the norm for battery data (most
wakes cost the same small integer percent), and exactly the shape a
single spliced-in glitch among a steady discharge rate has, so the
naive case would have let the outlier this is for sail straight
through. Falls back to mean absolute deviation there instead.

Verified against constructed glitch scenarios in both directions
(spurious dip and spurious spike): estimate now comes out identical
to the same series with the glitch removed entirely.

server/app/routers/common.py

@@ -5,6 +5,7 @@ from __future__ import annotations
 import io
 import logging
 import os
+import statistics
 import time
 from datetime import datetime, timedelta
 from urllib.parse import urlparse
@@ -40,6 +41,11 @@ RECHARGE_JUMP_PCT = 5  # a report this much above the recent baseline = battery
 RECHARGE_LOOKBACK = 3
 BATTERY_ESTIMATE_SAMPLE_COUNT = 100  # most recent battery_log rows considered
 MIN_ESTIMATE_SAMPLES = 5  # discharge steps needed before trusting the average
+# Modified z-score cutoff (Iglewicz & Hoaglin's standard figure) for
+# _reject_outlier_drops -- see that function's docstring for why a
+# single noisy reading needs rejecting at the per-wake-drop level, not
+# just at the recharge-detection level.
+OUTLIER_MODIFIED_Z_THRESHOLD = 3.5
 
 # "Overdue" threshold multiplier: the device should check in roughly every
 # refresh_interval_s; give it half again as long before flagging it.
@@ -128,6 +134,50 @@ def _avg_wake_interval_s(frame: Frame) -> float:
     return 86400 / wakes_per_day
 
 
+def _reject_outlier_drops(steps: list[tuple[int, float]]) -> list[tuple[int, float]]:
+    """Drops (weight, drop_pct) pairs whose drop is a wild outlier
+    relative to the rest of the recent steps. A single noisy ADC/
+    regulator glitch (see firmware/main/battery.c) corrupts one of the
+    two steps around it, whichever way it reads: a glitch that dips low
+    then recovers makes the step INTO it a spurious huge drop (the step
+    back out is an increase, already excluded above as a "recharge");
+    one that spikes high then settles makes the step OUT OF it the
+    spurious one instead (the step into it is the excluded "recharge").
+    Either way, one bad reading survives the recharge filter looking
+    like an ordinary, legitimately huge drop and swings the whole
+    remaining-time estimate on its own.
+
+    Uses a MAD-based modified z-score (robust to a small number of
+    extreme values in a way a plain mean/stdev z-score isn't -- a single
+    huge outlier inflates the stdev itself, which just hides the outlier
+    from a stdev-based test) rather than a fixed percent-point cutoff, so
+    it adapts to how noisy a given frame's own sensor actually is
+    instead of guessing one global threshold for every install."""
+    drops = [drop for _, drop in steps]
+    median = statistics.median(drops)
+    abs_devs = [abs(d - median) for d in drops]
+    mad = statistics.median(abs_devs)
+    if mad == 0:
+        # The standard median-based MAD degenerates to exactly 0 as soon
+        # as more than half the steps share the median exactly -- and
+        # real battery data is small integer percents, so "most wakes
+        # cost exactly 1%" ties are the norm, not an edge case. That's
+        # precisely the shape a single spliced-in glitch among a steady
+        # discharge rate has (18 steps at "1", one at "26"), so treating
+        # MAD==0 as "no spread, nothing to reject" would let exactly the
+        # outlier this function exists for sail straight through. Fall
+        # back to mean absolute deviation instead, which only reaches 0
+        # when every single step is identical.
+        mad = statistics.mean(abs_devs)
+        if mad == 0:
+            return steps  # every step really is identical -- nothing to reject
+    kept = [
+        (weight, drop) for weight, drop in steps
+        if abs(0.6745 * (drop - median) / mad) <= OUTLIER_MODIFIED_Z_THRESHOLD
+    ]
+    return kept or steps  # never filter down to nothing
+
+
 def battery_estimate_s(frame: Frame, db: Session) -> int | None:
     """Remaining-time estimate from a recency-weighted average of the
     *per-wake* percent drop, over the last BATTERY_ESTIMATE_SAMPLE_COUNT
@@ -143,10 +193,13 @@ def battery_estimate_s(frame: Frame, db: Session) -> int | None:
     entirely rather than folded in as a weird outlier; a flat step
     (0% change) still counts as a real, cheap wake -- excluding those
     would systematically overstate the per-wake cost by only counting
-    the wakes that happened to tick the percentage down. Steps are
-    weighted linearly by recency (step i of n gets weight i, 1-indexed)
-    so a recent change in usage pattern shows up quickly instead of
-    being washed out by a long flat history.
+    the wakes that happened to tick the percentage down. The remaining
+    steps then get one more pass, _reject_outlier_drops, to catch the
+    single-noisy-reading case that "percent went up" alone can't (see
+    that function's docstring). Steps are weighted linearly by recency
+    (step i of n gets weight i, 1-indexed) so a recent change in usage
+    pattern shows up quickly instead of being washed out by a long flat
+    history.
 
     The resulting %/wake rate is then converted to wall-clock time using
     the frame's *current* refresh_interval_s and quiet-hours settings
@@ -167,21 +220,21 @@ def battery_estimate_s(frame: Frame, db: Session) -> int | None:
         return None
     percents = list(reversed(rows))  # chronological order
 
-    weighted_drop_total = 0.0
-    weight_total = 0.0
-    valid_steps = 0
+    steps: list[tuple[int, float]] = []  # (recency_weight, drop_pct)
     for i in range(1, len(percents)):
         prev_pct, next_pct = percents[i - 1], percents[i]
         if next_pct > prev_pct:
             continue  # recharge (or a swap) -- not a discharge sample
-        weight = i  # later steps (larger i) count more
-        weighted_drop_total += weight * (prev_pct - next_pct)
-        weight_total += weight
-        valid_steps += 1
+        steps.append((i, prev_pct - next_pct))  # later steps (larger i) weigh more
 
-    if valid_steps < MIN_ESTIMATE_SAMPLES or weight_total <= 0:
+    if len(steps) < MIN_ESTIMATE_SAMPLES:
         return None
-    avg_drop_per_wake = weighted_drop_total / weight_total
+    steps = _reject_outlier_drops(steps)
+
+    weight_total = sum(weight for weight, _ in steps)
+    if weight_total <= 0:
+        return None
+    avg_drop_per_wake = sum(weight * drop for weight, drop in steps) / weight_total
     if avg_drop_per_wake <= 0:
         return None  # flat -- no honest rate to extrapolate