auto select knots using aicc

Author: CamDavidsonPilon

core/pioreactor/calibrations/protocols/od_reference_standard.py

@@ -256,7 +256,7 @@ def get_reference_standard_step(
     return get_session_step(_REFERENCE_STANDARD_STEPS, session, executor)
 
 
-def get_valid_od_devices_for_this_unit() -> list[str]:
+def get_valid_od_devices_for_this_unit() -> list[pt.ODCalibrationDevices]:
 
     pd_channels = config["od_config.photodiode_channel"]
     valid_devices: list[pt.ODCalibrationDevices] = []

core/pioreactor/calibrations/protocols/od_single_vial.py

@@ -484,7 +484,7 @@ def run_od_calibration(target_device: pt.ODCalibrationDevices) -> structs.OD600C
 
 
 class SingleVialODProtocol(CalibrationProtocol[pt.ODCalibrationDevices]):
-    target_device = pt.OD_DEVICES
+    target_device = pt.OD_DEVICES  # type: ignore
     protocol_name = "single_vial"
     description = "Calibrate OD using a single vial"
 

core/pioreactor/calibrations/protocols/od_standards.py

@@ -187,8 +187,7 @@ def _calculate_curve_data(
     if len(od600_values) >= 3:
         from pioreactor.utils.splines import spline_fit
 
-        knots = min(4, len(od600_values))
-        return "spline", spline_fit(od600_values, voltages, knots=knots, weights=weights)
+        return "spline", spline_fit(od600_values, voltages, knots="auto", weights=weights)
 
     degree = min(3, max(1, len(od600_values) - 1))
     return "poly", utils.calculate_poly_curve_of_best_fit(od600_values, voltages, degree, weights)
@@ -598,7 +597,7 @@ def get_standards_step(
     return get_session_step(_OD_STANDARDS_STEPS, session, executor)
 
 
-def get_valid_od_devices_for_this_unit() -> list[str]:
+def get_valid_od_devices_for_this_unit() -> list[pt.ODCalibrationDevices]:
     pd_channels = config["od_config.photodiode_channel"]
     valid_devices: list[pt.ODCalibrationDevices] = []
 

core/pioreactor/calibrations/protocols/pump_duration_based.py

@@ -580,7 +580,7 @@ def run_pump_calibration(
 
 
 class DurationBasedPumpProtocol(CalibrationProtocol[pt.PumpCalibrationDevices]):
-    target_device = pt.PUMP_DEVICES
+    target_device = cast(list[pt.PumpCalibrationDevices], pt.PUMP_DEVICES)
     protocol_name = "duration_based"
     title = "Duration-based pump calibration"
     description = "Build a duration-to-volume curve for the {device} pump using a simple multi-step flow."

core/pioreactor/calibrations/protocols/stirring_dc_based.py

@@ -140,14 +140,13 @@ def _build_stirring_calibration_from_measurements(
 
     from pioreactor.utils.splines import spline_fit
 
-    knots = min(3, len(dcs))
     return SimpleStirringCalibration(
         pwm_hz=config.getfloat("stirring.config", "pwm_hz"),
         voltage=voltage,
         calibration_name=f"stirring-calibration-{current_utc_datetime().strftime('%Y-%m-%d_%H-%M')}",
         calibrated_on_pioreactor_unit=unit,
         created_at=current_utc_datetime(),
-        curve_data_=spline_fit(dcs, rpms, knots=knots),
+        curve_data_=spline_fit(dcs, rpms, knots="auto"),
         curve_type="spline",
         recorded_data={"x": dcs, "y": rpms},
     )

core/pioreactor/calibrations/registry.py

@@ -16,7 +16,7 @@
 
 class CalibrationProtocol(Generic[Device]):
     protocol_name: ClassVar[ProtocolName]
-    target_device: ClassVar[str | list[str]]
+    target_device: ClassVar[str | list[Device]]
     title: ClassVar[str] = ""
     description: ClassVar[str] = ""
     requirements: ClassVar[tuple[str, ...]] = ()

core/pioreactor/utils/polys.py

@@ -6,6 +6,11 @@
 import numpy as np
 
 
+def _to_pyfloat(seq: list[float]) -> list[float]:
+    # we have trouble serializing numpy floats
+    return [float(_) for _ in seq]
+
+
 def poly_fit(
     x: Sequence[float],
     y: Sequence[float],
@@ -33,7 +38,7 @@ def poly_fit(
             raise ValueError("weights must be non-negative.")
 
     coefs = np.polyfit(x_values, y_values, deg=degree, w=weight_values)
-    return coefs.tolist()
+    return _to_pyfloat(coefs.tolist())
 
 
 def poly_eval(poly_data: list[float], x: float) -> float:

core/pioreactor/utils/splines.py

@@ -7,10 +7,15 @@
 import numpy as np
 
 
+def _to_pyfloat(seq: list[float]) -> list[float]:
+    # we have trouble serializing numpy floats
+    return [float(_) for _ in seq]
+
+
 def spline_fit(
     x: Sequence[float],
     y: Sequence[float],
-    knots: int | Sequence[float],
+    knots: int | Sequence[float] | str | None = "auto",
     weights: Sequence[float] | None = None,
 ) -> list:
     """
@@ -21,7 +26,8 @@ def spline_fit(
     x, y
         Observations.
     knots
-        Either the number of knots to use (including boundaries) or explicit knot positions.
+        Either the number of knots to use (including boundaries), explicit knot positions, or "auto".
+        When "auto" (default), knot count is selected by AICc over a small candidate range.
     weights
         Optional weights for each observation.
 
@@ -41,10 +47,6 @@ def spline_fit(
     if np.allclose(x_values, x_values[0]):
         raise ValueError("x values must not all be the same.")
 
-    knot_positions = _normalize_knots(x_values, knots)
-    if len(knot_positions) < 2:
-        raise ValueError("At least two knots are required.")
-
     if weights is None:
         weight_values = np.ones_like(x_values)
     else:
@@ -54,14 +56,19 @@ def spline_fit(
         if np.any(weight_values < 0):
             raise ValueError("weights must be non-negative.")
 
-    design_matrix = _build_spline_design_matrix(knot_positions, x_values)
-    weighted_design = design_matrix * np.sqrt(weight_values)[:, None]
-    weighted_y = y_values * np.sqrt(weight_values)
+    if knots is None or knots == "auto":
+        knot_positions = _auto_select_knots(x_values, y_values, weight_values)
+    elif isinstance(knots, str):
+        raise ValueError('knots must be an int, a sequence of floats, or "auto".')
+    else:
+        knot_positions = _normalize_knots(x_values, knots)
+    if len(knot_positions) < 2:
+        raise ValueError("At least two knots are required.")
 
-    knot_values, *_ = np.linalg.lstsq(weighted_design, weighted_y, rcond=None)
+    knot_values, _ = _fit_knot_values(knot_positions, x_values, y_values, weight_values)
     coefficients = _natural_cubic_spline_coefficients(knot_positions, knot_values)
 
-    return [knot_positions.tolist(), [coeff.tolist() for coeff in coefficients]]
+    return [_to_pyfloat(knot_positions.tolist()), [_to_pyfloat(coeff.tolist()) for coeff in coefficients]]
 
 
 def spline_eval(spline_data: list, x: float) -> float:
@@ -133,6 +140,61 @@ def _normalize_knots(x_values: np.ndarray, knots: int | Sequence[float]) -> np.n
     return knot_positions
 
 
+def _fit_knot_values(
+    knot_positions: np.ndarray,
+    x_values: np.ndarray,
+    y_values: np.ndarray,
+    weight_values: np.ndarray,
+) -> tuple[np.ndarray, np.ndarray]:
+    design_matrix = _build_spline_design_matrix(knot_positions, x_values)
+    sqrt_weights = np.sqrt(weight_values)
+    weighted_design = design_matrix * sqrt_weights[:, None]
+    weighted_y = y_values * sqrt_weights
+    knot_values, *_ = np.linalg.lstsq(weighted_design, weighted_y, rcond=None)
+    return knot_values, design_matrix
+
+
+def _aicc_score(weighted_sse: float, n_obs: int, n_params: int) -> float:
+    if n_obs <= n_params + 1:
+        return float("inf")
+    sse = max(weighted_sse, np.finfo(float).tiny)
+    correction = (2 * n_params * (n_params + 1)) / (n_obs - n_params - 1)
+    return n_obs * np.log(sse / n_obs) + 2 * n_params + correction
+
+
+def _auto_select_knots(
+    x_values: np.ndarray,
+    y_values: np.ndarray,
+    weight_values: np.ndarray,
+    *,
+    max_knots: int | None = None,
+) -> np.ndarray:
+    n_obs = x_values.size
+    unique_x = np.unique(x_values).size
+    if max_knots is None:
+        max_knots = min(6, n_obs)
+    max_knots = min(max_knots, unique_x)
+    max_knots = max(2, max_knots)
+
+    best_score = float("inf")
+    best_knots: np.ndarray | None = None
+
+    for count in range(2, max_knots + 1):
+        knot_positions = _normalize_knots(x_values, count)
+        knot_values, design_matrix = _fit_knot_values(knot_positions, x_values, y_values, weight_values)
+        y_pred = design_matrix @ knot_values
+        residual = y_values - y_pred
+        weighted_sse = float(np.sum(weight_values * residual**2))
+        score = _aicc_score(weighted_sse, n_obs, knot_positions.size)
+        if score < best_score:
+            best_score = score
+            best_knots = knot_positions
+
+    if best_knots is None:
+        return _normalize_knots(x_values, 2)
+    return best_knots
+
+
 def _build_spline_design_matrix(knots: np.ndarray, x_values: np.ndarray) -> np.ndarray:
     n = x_values.size
     m = knots.size

core/tests/test_splines.py

@@ -15,6 +15,16 @@ def test_spline_fit_and_eval_linear() -> None:
     assert spline_eval(spline_data, -1.0) == pytest.approx(-1.0, rel=1e-6)
 
 
+def test_spline_fit_auto_selects_knots() -> None:
+    x = [0.0, 1.0, 2.0, 3.0]
+    y = [1.0, 3.0, 5.0, 7.0]
+    spline_data = spline_fit(x, y, knots="auto")
+
+    assert spline_eval(spline_data, 2.5) == pytest.approx(6.0, rel=1e-6)
+    assert spline_data[0][0] == pytest.approx(min(x), rel=1e-6)
+    assert spline_data[0][-1] == pytest.approx(max(x), rel=1e-6)
+
+
 def test_spline_fit_explicit_knots_interpolate_at_knots() -> None:
     x = [0.0, 1.0, 2.0]
     y = [0.0, 1.0, 0.0]
@@ -133,6 +143,19 @@ def test_spline_fit_reduces_residuals_on_noisy_linear_data() -> None:
     assert mse < 0.2
 
 
+def test_spline_fit_auto_selects_linear_curve_for_noisy_linear_data() -> None:
+    rng = np.random.default_rng(123)
+    x = np.linspace(0.0, 10.0, 30)
+    y = 2.5 * x - 1.0 + rng.normal(0.0, 0.1, size=x.size)
+    spline_data = spline_fit(x.tolist(), y.tolist(), knots="auto")
+
+    knots, coefficients = spline_data
+    assert len(knots) == 2
+    assert len(coefficients) == 1
+    assert coefficients[0][2] == pytest.approx(0.0, abs=1e-12)
+    assert coefficients[0][3] == pytest.approx(0.0, abs=1e-12)
+
+
 def test_spline_fit_respects_sorted_or_unsorted_input() -> None:
     x = [0.0, 1.0, 2.0, 3.0]
     y = [1.0, 2.0, 0.0, 3.0]