fix calibration for signals out of range

pioreactor/actions/od_calibration.py

@@ -172,8 +172,8 @@ def plot_data(
 
 def start_recording_and_diluting(
     st: Stirrer,
-    initial_od600: float,
-    minimum_od600: float,
+    initial_od600: pt.OD,
+    minimum_od600: pt.OD,
     dilution_amount: float,
     signal_channel,
 ):
@@ -194,7 +194,7 @@ def start_recording_and_diluting(
         use_calibration=False,
     ) as od_reader:
 
-        def get_voltage_from_adc() -> float:
+        def get_voltage_from_adc() -> pt.Voltage:
             od_readings1 = od_reader.record_from_adc()
             od_readings2 = od_reader.record_from_adc()
             return 0.5 * (od_readings1.ods[signal_channel].od + od_readings2.ods[signal_channel].od)
@@ -285,19 +285,19 @@ def get_voltage_from_adc() -> float:
             x_max=initial_od600,
         )
         click.echo("Empty the vial and replace with 10 mL of the media you used.")
-        inferred_od600 = click.prompt("What is the OD600 of your blank?", type=float)
+        od600_of_blank = click.prompt("What is the OD600 of your blank?", type=float)
         click.echo("Confirm vial outside is dry and clean. Place back into Pioreactor.")
         while not click.confirm("Continue?", default=True):
             pass
 
         voltages.append(get_voltage_from_adc())
-        inferred_od600s.append(inferred_od600)
+        inferred_od600s.append(od600_of_blank)
 
         return inferred_od600s, voltages
 
 
 def calculate_curve_of_best_fit(
-    voltages: list[float], inferred_od600s: list[float], degree: int
+    voltages: list[pt.Voltage], inferred_od600s: list[pt.OD], degree: int
 ) -> tuple[list[float], str]:
     import numpy as np
 
@@ -320,8 +320,8 @@ def calculate_curve_of_best_fit(
 def show_results_and_confirm_with_user(
     curve_data: list[float],
     curve_type: str,
-    voltages: list[float],
-    inferred_od600s: list[float],
+    voltages: list[pt.Voltage],
+    inferred_od600s: list[pt.OD],
 ) -> tuple[bool, int]:
     click.clear()
 
@@ -361,12 +361,10 @@ def show_results_and_confirm_with_user(
 def save_results(
     curve_data_: list[float],
     curve_type: str,
-    voltages: list[float],
-    od600s: list[float],
+    voltages: list[pt.Voltage],
+    od600s: list[pt.OD],
     angle,
     name: str,
-    maximum_od600: float,
-    minimum_od600: float,
     signal_channel: pt.PdChannel,
     unit: str,
 ) -> structs.ODCalibration:
@@ -386,8 +384,8 @@ def save_results(
         pioreactor_unit=unit,
         name=name,
         angle=angle,
-        maximum_od600=maximum_od600,
-        minimum_od600=0,
+        maximum_od600=max(od600s),
+        minimum_od600=min(od600s),
         minimum_voltage=min(voltages),
         maximum_voltage=max(voltages),
         curve_data_=curve_data_,
@@ -431,7 +429,7 @@ def od_calibration() -> None:
                 st, initial_od600, minimum_od600, dilution_amount, signal_channel
             )
 
-        degree = 4
+        degree = 5 if len(voltages) > 5 else 3
         while True:
             curve_data_, curve_type = calculate_curve_of_best_fit(voltages, inferred_od600s, degree)
             okay_with_result, degree = show_results_and_confirm_with_user(
@@ -447,8 +445,6 @@ def od_calibration() -> None:
             inferred_od600s,
             angle,
             name,
-            initial_od600,
-            minimum_od600,
             signal_channel,
             unit,
         )

pioreactor/automations/dosing/base.py

@@ -295,6 +295,9 @@ def run(self, timeout: float = 60.0) -> Optional[events.AutomationEvent]:
 
         """
         event: Optional[events.AutomationEvent]
+
+        self._latest_run_at = current_utc_datetime()
+
         if self.state == self.DISCONNECTED:
             # NOOP
             # we ended early.
@@ -326,7 +329,6 @@ def run(self, timeout: float = 60.0) -> Optional[events.AutomationEvent]:
             self.logger.info(str(event))
 
         self.latest_event = event
-        self._latest_run_at = current_utc_datetime()
         return event
 
     def block_until_not_sleeping(self) -> bool:

pioreactor/background_jobs/dosing_control.py

@@ -65,6 +65,9 @@ def __init__(self, unit: str, experiment: str, automation_name: str, **kwargs) -
         try:
             automation_class = self.available_automations[automation_name]
         except KeyError:
+            self.logger.error(
+                f"Unable to find automation {automation_name}. Available automations are {list(self.available_automations.keys())}"
+            )
             self.clean_up()
             raise KeyError(
                 f"Unable to find automation {automation_name}. Available automations are {list(self.available_automations.keys())}"

pioreactor/background_jobs/led_control.py

@@ -43,6 +43,10 @@ def __init__(self, unit: str, experiment: str, automation_name: str, **kwargs) -
         try:
             automation_class = self.available_automations[automation_name]
         except KeyError:
+            self.logger.error(
+                f"Unable to find automation {automation_name}. Available automations are {list(self.available_automations.keys())}"
+            )
+            self.clean_up()
             raise KeyError(
                 f"Unable to find automation {automation_name}. Available automations are {list(self.available_automations.keys())}"
             )

pioreactor/background_jobs/monitor.py

@@ -448,17 +448,14 @@ def rpi_is_having_power_problems(self) -> tuple[bool, float]:
         from pioreactor.hardware import voltage_in_aux
 
         voltage_read = voltage_in_aux(precision=0.05)
-        if voltage_read <= 4.80:
-            return (True, voltage_read)
-
-        under_voltage = new_under_voltage()
-        if under_voltage is None:
-            # not supported on system
-            return (False, voltage_read)
-        elif under_voltage.get():
-            return (True, voltage_read)
+        under_voltage_flag = new_under_voltage()
+
+        under_voltage_status = under_voltage_flag.get() if under_voltage_flag else None
+
+        if voltage_read <= 4.80 and (under_voltage_status is None or under_voltage_status):
+            return True, voltage_read
         else:
-            return (False, voltage_read)
+            return False, voltage_read
 
     def check_for_power_problems(self) -> None:
         is_rpi_having_power_probems, voltage = self.rpi_is_having_power_problems()

pioreactor/background_jobs/od_reading.py

@@ -655,7 +655,7 @@ def __init__(self) -> None:
 class CachedCalibrationTransformer(CalibrationTransformer):
     def __init__(self, channel_angle_map: dict[pt.PdChannel, pt.PdAngle]) -> None:
         super().__init__()
-
+        self.has_logged_warning = False
         try:
             self.models = self.get_models_from_disk(channel_angle_map)
         except Exception as e:
@@ -711,17 +711,22 @@ def _hydrate_model(self, calibration_data: structs.ODCalibration) -> Callable[[f
             this procedure effectively ignores it.
 
             """
-            from numpy import roots, zeros_like, iscomplex, real
+            from numpy import roots, zeros_like, real
 
-            def calibration(x: float) -> float:
+            def calibration(observed_voltage: pt.Voltage) -> pt.OD:
                 poly = calibration_data.curve_data_
+                min_OD, max_OD = calibration_data.minimum_od600, calibration_data.maximum_od600
+                min_voltage, max_voltage = (
+                    calibration_data.minimum_voltage,
+                    calibration_data.maximum_voltage,
+                )
 
-                coef_shift = zeros_like(calibration_data.curve_data_)
-                coef_shift[-1] = x
-                roots_ = roots(poly - coef_shift)
-                min_OD, max_OD = 0, float(calibration_data.maximum_od600)
+                coef_shift = zeros_like(poly)
+                coef_shift[-1] = observed_voltage
+                solve_for_poly = poly - coef_shift
+                roots_ = roots(solve_for_poly)
                 plausible_roots_ = sorted(
-                    [real(r) for r in roots_ if not iscomplex(r) and (min_OD <= r <= max_OD)]
+                    [real(r) for r in roots_ if (r.imag == 0) and (min_OD <= real(r) <= max_OD)]
                 )
 
                 try:
@@ -730,13 +735,29 @@ def calibration(x: float) -> float:
                     return ideal_root
 
                 except IndexError:
-                    self.logger.debug(f"Outside suggested calibration range [{min_OD}, {max_OD}].")
-                    return max_OD
+                    if observed_voltage <= min_voltage:
+                        # voltage less than the blank recorded during the calibration and the calibration curve doesn't have solutions (ex even-deg poly)
+                        # this isn't great, as there is nil noise in the signal.
+
+                        if not self.has_logged_warning:
+                            self.logger.warning(
+                                f"Signal below suggested calibration range. Calibrated for OD=[{min_OD:0.3g}, {max_OD:0.3g}], V=[{min_voltage:0.3g}, {max_voltage:0.3g}]. Observed {observed_voltage:0.3f}V."
+                            )
+                            self.has_logged_warning = True
+                        return min_OD
+
+                    else:
+                        if not self.has_logged_warning:
+                            self.logger.warning(
+                                f"Signal outside suggested calibration range. Calibrated for OD=[{min_OD:0.3g}, {max_OD:0.3g}], V=[{min_voltage:0.3g}, {max_voltage:0.3g}]. Observed {observed_voltage:0.3f}V."
+                            )
+                            self.has_logged_warning = True
+                        return max_OD
 
         else:
 
-            def calibration(x: float) -> float:
-                return x
+            def calibration(observed_voltage: pt.Voltage) -> pt.OD:
+                return observed_voltage
 
         return calibration
 

pioreactor/background_jobs/temperature_control.py

@@ -135,6 +135,10 @@ def __init__(
         try:
             automation_class = self.available_automations[automation_name]
         except KeyError:
+            self.logger.error(
+                f"Unable to find automation {automation_name}. Available automations are {list(self.available_automations.keys())}"
+            )
+            self.clean_up()
             raise KeyError(
                 f"Unable to find automation {automation_name}. Available automations are {list(self.available_automations.keys())}"
             )

pioreactor/structs.py

@@ -227,14 +227,14 @@ class WastePumpCalibration(PumpCalibration, tag="waste_pump"):
 
 class ODCalibration(Calibration):
     angle: pt.PdAngle
-    maximum_od600: float
-    minimum_od600: float
-    minimum_voltage: float
-    maximum_voltage: float
+    maximum_od600: pt.OD
+    minimum_od600: pt.OD
+    minimum_voltage: pt.Voltage
+    maximum_voltage: pt.Voltage
     curve_type: str
     curve_data_: list[float]
-    voltages: list[float]
-    od600s: list[float]
+    voltages: list[pt.Voltage]
+    od600s: list[pt.OD]
     ir_led_intensity: float
     pd_channel: pt.PdChannel
 

pioreactor/tests/test_od_reading.py

@@ -597,7 +597,7 @@ def test_calibration_simple_linear_calibration():
         pause()
         pause()
         pause()
-        with collect_all_logs_of_level("debug", unit=get_unit_name(), experiment="+") as bucket:
+        with collect_all_logs_of_level("warning", unit=get_unit_name(), experiment="+") as bucket:
             voltage = 10.0
             pause()
             pause()
@@ -646,7 +646,7 @@ def test_calibration_simple_linear_calibration_negative_slope():
         voltage = 0.5
         assert od.calibration_transformer.models["2"](voltage) == (voltage - 2) / (-0.1)
 
-        with collect_all_logs_of_level("debug", unit=get_unit_name(), experiment="+") as bucket:
+        with collect_all_logs_of_level("warning", unit=get_unit_name(), experiment="+") as bucket:
             voltage = 12.0
             assert voltage > 2.0
 
@@ -839,3 +839,90 @@ def test_ODReader_with_multiple_angles_and_a_ref():
             print(signal)
             if i == 3:
                 break
+
+
+def test_calibration_data_from_user1():
+    # the problem is that the 4th degree polynomial doesn't always have a solution to the inverse problem.
+    experiment = "test_calibration_data_from_user1"
+    poly = [2.583, -3.447, 1.531, 0.223, 0.017]  # email correspondence
+
+    with local_persistant_storage("current_od_calibration") as c:
+        c["90"] = encode(
+            structs.OD90Calibration(
+                created_at=current_utc_datetime(),
+                curve_type="poly",
+                curve_data_=poly,
+                name="multi_test",
+                maximum_od600=1.0,
+                minimum_od600=0.01,
+                ir_led_intensity=90.0,
+                angle="90",
+                minimum_voltage=0.018,
+                maximum_voltage=1.0,
+                voltages=[],
+                od600s=[],
+                pd_channel="2",
+                pioreactor_unit=get_unit_name(),
+            )
+        )
+
+    with start_od_reading("REF", "90", interval=None, fake_data=True, experiment=experiment) as od:
+        assert isinstance(od.calibration_transformer, CachedCalibrationTransformer)
+        infer = od.calibration_transformer.models["2"]
+
+        # try varying voltage up over and across the lower bound, and assert we are always non-decreasing.
+        od_0 = 0
+        for i in range(10):
+            voltage = i / 5 * 0.018
+            od_1 = infer(voltage)
+            assert od_0 <= od_1
+            od_0 = od_1
+
+    with local_persistant_storage("current_od_calibration") as c:
+        del c["90"]
+
+
+def test_calibration_data_from_user2():
+    # the difference here is that the 3 degree polynomial always has a solution to the inverse problem.
+    experiment = "test_calibration_data_from_user2"
+    poly = [
+        1.71900012,
+        -1.77900665,
+        0.95000656,
+        -0.01770485,
+    ]  # looks like the degree 4 above: https://chat.openai.com/share/2ef30900-22ef-4a7f-8f34-14a88ffc65a8
+
+    with local_persistant_storage("current_od_calibration") as c:
+        c["90"] = encode(
+            structs.OD90Calibration(
+                created_at=current_utc_datetime(),
+                curve_type="poly",
+                curve_data_=poly,
+                name="multi_test",
+                maximum_od600=1.0,
+                minimum_od600=0.01,
+                ir_led_intensity=90.0,
+                angle="90",
+                minimum_voltage=0.018,
+                maximum_voltage=1.0,
+                voltages=[],
+                od600s=[],
+                pd_channel="2",
+                pioreactor_unit=get_unit_name(),
+            )
+        )
+
+    with start_od_reading("REF", "90", interval=None, fake_data=True, experiment=experiment) as od:
+        assert isinstance(od.calibration_transformer, CachedCalibrationTransformer)
+        infer = od.calibration_transformer.models["2"]
+
+        # try varying voltage up over and across the lower bound, and assert we are always non-decreasing.
+        od_0 = 0
+        for i in range(10):
+            voltage = i / 5 * 0.018
+            od_1 = infer(voltage)
+            assert od_0 <= od_1
+            od_0 = od_1
+
+    with local_persistant_storage("current_od_calibration") as c:
+        del c["90"]

pioreactor/types.py

@@ -90,6 +90,7 @@ class PublishableSetting(t.TypedDict, total=False):
 # hardware level stuff
 AnalogValue = t.Union[int, float]
 Voltage = float  # maybe should be non-negative?
+OD = float  # maybe should be non-negative?
 
 AdcChannel = t.Literal[0, 1, 2, 3]
 

requirements/requirements.txt

@@ -5,5 +5,5 @@ sh==1.14.2
 JSON-log-formatter==0.4.0
 rpi_hardware_pwm==0.1.3
 colorlog==6.6.0
-msgspec==0.18.1
+msgspec==0.18.2
 diskcache==5.6.1