From b6ba626e5694ef3be186b02f096ab9aed45ddb90 Mon Sep 17 00:00:00 2001 From: Awesome-Ace22 Date: Tue, 3 Nov 2026 23:26:13 -0500 Subject: [PATCH] feat(data_service): add new linear calibration methods and corresponding tests --- app/models.py | 3 ++ app/services/data_service.py | 73 +++++++++++++++++++++++++++- tests/test_robust_calibration.py | 82 ++++++++++++++++++++++++++++++++ 3 files changed, 156 insertions(+), 2 deletions(-) create mode 100644 tests/test_robust_calibration.py diff --git a/app/models.py b/app/models.py index 7846ffe..a0ad2a9 100644 --- a/app/models.py +++ b/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): diff --git a/app/services/data_service.py b/app/services/data_service.py index 11e6248..7b348a4 100644 --- a/app/services/data_service.py +++ b/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( diff --git a/tests/test_robust_calibration.py b/tests/test_robust_calibration.py new file mode 100644 index 0000000..c588b2d --- /dev/null +++ b/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"