feat(data_service): add new linear calibration methods

app/models.py

@@ -52,6 +52,9 @@ class ConvertMethod(str, Enum):
     ZERO_OFFSET = "zero_offset"
     LINEAR = "linear"
     POLYNOMIAL = "polynomial"
+    # Robust straight-line fits: a bad calibration point doesn't pull the line.
+    THEIL_SEN = "theil_sen"
+    REPEATED_MEDIAN = "repeated_median"
 
 
 class ConvertSpec(BaseModel):

app/services/data_service.py

@@ -2,11 +2,12 @@
 
 from collections import defaultdict, deque
 from dataclasses import dataclass
+from functools import lru_cache
 from typing import Iterable
 
 import numpy as np
 
-from app.models import ConvertMethod, FasSensorBinding, GcsSensorBinding, SensorEntry, SourceTarget, SystemConfig
+from app.models import ConvertMethod, ConvertSpec, FasSensorBinding, GcsSensorBinding, SensorEntry, SourceTarget, SystemConfig
 
 
 def linear_interpolate(raw_value: float, points: list[tuple[float, float]] | None, degree: int = 1) -> float:
@@ -22,6 +23,74 @@ def linear_interpolate(raw_value: float, points: list[tuple[float, float]] | Non
     return m * raw_value + b
 
 
+def _pairwise_slopes(xs: np.ndarray, ys: np.ndarray) -> np.ndarray:
+    """slopes[i, j] is the slope from point i to point j; NaN where i == j or the x values match."""
+    dx = xs[None, :] - xs[:, None]
+    dy = ys[None, :] - ys[:, None]
+    with np.errstate(divide="ignore", invalid="ignore"):
+        slopes = dy / dx
+    slopes[dx == 0] = np.nan
+    return slopes
+
+
+def theil_sen_fit(xs: np.ndarray, ys: np.ndarray) -> tuple[float, float] | None:
+    """Slope is the median of the slopes between every pair of points. Ignores up to ~29% bad points."""
+    slopes = _pairwise_slopes(xs, ys)[np.triu_indices(len(xs), k=1)]
+    slopes = slopes[~np.isnan(slopes)]
+    if slopes.size == 0:
+        return None
+    m = float(np.median(slopes))
+    return m, float(np.median(ys - m * xs))
+
+
+def repeated_median_fit(xs: np.ndarray, ys: np.ndarray) -> tuple[float, float] | None:
+    """Siegel's repeated median: each point's median slope to the others, then the median of those.
+
+    Ignores up to half the points being bad.
+    """
+    slopes = _pairwise_slopes(xs, ys)
+    has_slope = ~np.all(np.isnan(slopes), axis=1)
+    if not has_slope.any():
+        return None
+    m = float(np.median(np.nanmedian(slopes[has_slope], axis=1)))
+    return m, float(np.median(ys - m * xs))
+
+
+_ROBUST_FITS = {
+    ConvertMethod.THEIL_SEN: theil_sen_fit,
+    ConvertMethod.REPEATED_MEDIAN: repeated_median_fit,
+}
+
+
+@lru_cache(maxsize=512)
+def _robust_line(method: ConvertMethod, points: tuple[tuple[float, float], ...]) -> tuple[float, float] | None:
+    # Cached so a sensor's table is fitted once, not on every sample.
+    calibration_points = np.asarray(points, dtype=float)
+    if calibration_points.ndim != 2 or calibration_points.shape[1] != 2:
+        raise ValueError("Calibration must be a list of [voltage, reading] pairs.")
+    return _ROBUST_FITS[method](calibration_points[:, 0], calibration_points[:, 1])
+
+
+def robust_interpolate(raw_value: float, points: list[tuple[float, float]] | None, method: ConvertMethod) -> float:
+    """Apply a robust straight-line fit. Like `linear_interpolate`, the line carries on past the
+    outermost points, and fewer than two points with different raw values pass the raw value through."""
+    if not points or len(points) < 2:
+        return raw_value
+    line = _robust_line(method, tuple((float(x), float(y)) for x, y in points))
+    if line is None:
+        return raw_value
+    m, b = line
+    return m * raw_value + b
+
+
+def convert_value(raw_value: float, convert: ConvertSpec) -> float:
+    """Calibrate one raw reading with the sensor's conversion method (not called for `none`)."""
+    if convert.method in _ROBUST_FITS:
+        return robust_interpolate(raw_value, convert.calibration, convert.method)
+    # linear, and zero_offset / polynomial, which aren't implemented separately yet.
+    return linear_interpolate(raw_value, convert.calibration)
+
+
 @dataclass
 class ParsedSensor:
     name: str
@@ -102,7 +171,7 @@ def parse(self, source: str, raw_sensors: list[dict] | dict[str, dict], config:
 
             value = raw_value
             if calibration_enabled and sensor_cfg.convert.method != ConvertMethod.NONE:
-                value = linear_interpolate(raw_value, sensor_cfg.convert.calibration)
+                value = convert_value(raw_value, sensor_cfg.convert)
             averaged_value = self._rolling_store.add(sensor_cfg.name, value)
 
             parsed.append(

tests/test_robust_calibration.py

@@ -0,0 +1,82 @@
+import pytest
+
+from app.models import ConvertMethod, SystemConfig
+from app.services.data_service import (
+    RollingAverageStore,
+    SensorParser,
+    linear_interpolate,
+    robust_interpolate,
+)
+
+ROBUST = [ConvertMethod.THEIL_SEN, ConvertMethod.REPEATED_MEDIAN]
+
+# y = 10x, except the last point, which is badly off.
+ONE_BAD_POINT = [(0.0, 0.0), (1.0, 10.0), (2.0, 20.0), (3.0, 30.0), (4.0, 400.0)]
+# y = 10x for 5 of 7 points: more bad points than Theil-Sen can ignore, fewer than half.
+TWO_BAD_OF_SEVEN = [(0.0, 0.0), (1.0, 10.0), (2.0, 20.0), (3.0, 30.0), (4.0, 40.0), (5.0, 500.0), (6.0, 900.0)]
+
+
+@pytest.mark.parametrize("method", ROBUST)
+def test_clean_points_match_least_squares(method: ConvertMethod) -> None:
+    points = [(0.0, 0.0), (2.0, 100.0)]
+    assert robust_interpolate(1.0, points, method) == pytest.approx(linear_interpolate(1.0, points))
+
+
+@pytest.mark.parametrize("method", ROBUST)
+def test_one_bad_point_does_not_skew_the_line(method: ConvertMethod) -> None:
+    # Least squares is pulled well away from y = 10x; the robust fits aren't.
+    assert linear_interpolate(2.0, ONE_BAD_POINT) != pytest.approx(20.0, abs=1.0)
+    assert robust_interpolate(2.0, ONE_BAD_POINT, method) == pytest.approx(20.0)
+
+
+@pytest.mark.parametrize("method", ROBUST)
+def test_extrapolates_past_the_outermost_points(method: ConvertMethod) -> None:
+    assert robust_interpolate(-1.0, ONE_BAD_POINT, method) == pytest.approx(-10.0)
+    assert robust_interpolate(10.0, ONE_BAD_POINT, method) == pytest.approx(100.0)
+
+
+def test_repeated_median_ignores_one_bad_point_in_four() -> None:
+    points = [(0.0, 0.0), (1.0, 10.0), (2.0, 20.0), (3.0, 330.0)]
+    assert robust_interpolate(2.0, points, ConvertMethod.REPEATED_MEDIAN) == pytest.approx(20.0)
+
+
+def test_repeated_median_survives_more_bad_points_than_theil_sen() -> None:
+    assert robust_interpolate(2.0, TWO_BAD_OF_SEVEN, ConvertMethod.REPEATED_MEDIAN) == pytest.approx(20.0)
+    assert robust_interpolate(2.0, TWO_BAD_OF_SEVEN, ConvertMethod.THEIL_SEN) != pytest.approx(20.0, abs=1.0)
+
+
+@pytest.mark.parametrize("method", ROBUST)
+def test_repeated_raw_values_are_ignored_when_computing_slopes(method: ConvertMethod) -> None:
+    points = [(1.0, 10.0), (1.0, 12.0), (2.0, 20.0), (3.0, 30.0)]
+    assert robust_interpolate(4.0, points, method) == pytest.approx(40.0)
+
+
+@pytest.mark.parametrize("method", ROBUST)
+@pytest.mark.parametrize("points", [None, [], [(1.0, 1.0)], [(1.0, 1.0), (1.0, 5.0)]])
+def test_too_few_distinct_points_pass_the_raw_value_through(method: ConvertMethod, points) -> None:
+    assert robust_interpolate(7.3, points, method) == 7.3
+
+
+@pytest.mark.parametrize("method", ["theil_sen", "repeated_median"])
+def test_parser_applies_robust_methods(method: str) -> None:
+    config = SystemConfig.model_validate(
+        {
+            "Sensors": [
+                {
+                    "name": "PGSO",
+                    "type": "PT",
+                    "unit": "psi",
+                    "binding": {"source": "GCS", "hat_id": 0, "channel_id": 0},
+                    "convert": {"method": method, "calibration": [list(p) for p in ONE_BAD_POINT]},
+                }
+            ]
+        }
+    )
+    parsed = SensorParser(RollingAverageStore(window_size=5)).parse(
+        source="GCS",
+        raw_sensors=[{"hat_id": 0, "channel_id": 0, "value": 2.0, "timestamp": 1}],
+        config=config,
+        calibration_enabled=True,
+    )
+    assert parsed[0].value == pytest.approx(20.0)
+    assert parsed[0].unit == "psi"