From 8faf76d3925a78bb0e5e9d3a0ea197be446cd2e1 Mon Sep 17 00:00:00 2001 From: CamDavidsonPilon Date: Thu, 28 Sep 2023 13:33:54 -0400 Subject: [PATCH] fix calibration for signals out of range --- pioreactor/actions/od_calibration.py | 30 +++--- pioreactor/automations/dosing/base.py | 4 +- pioreactor/background_jobs/dosing_control.py | 3 + pioreactor/background_jobs/led_control.py | 4 + pioreactor/background_jobs/monitor.py | 17 ++-- pioreactor/background_jobs/od_reading.py | 45 ++++++--- .../background_jobs/temperature_control.py | 4 + pioreactor/structs.py | 12 +-- pioreactor/tests/test_od_reading.py | 91 ++++++++++++++++++- pioreactor/types.py | 1 + requirements/requirements.txt | 2 +- 11 files changed, 164 insertions(+), 49 deletions(-) diff --git a/pioreactor/actions/od_calibration.py b/pioreactor/actions/od_calibration.py index e8ab9f11..994d6fe3 100644 --- a/pioreactor/actions/od_calibration.py +++ b/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, ) diff --git a/pioreactor/automations/dosing/base.py b/pioreactor/automations/dosing/base.py index f2ed87f3..3f9d849f 100644 --- a/pioreactor/automations/dosing/base.py +++ b/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: diff --git a/pioreactor/background_jobs/dosing_control.py b/pioreactor/background_jobs/dosing_control.py index 6f2cbcd6..b3415a51 100644 --- a/pioreactor/background_jobs/dosing_control.py +++ b/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())}" diff --git a/pioreactor/background_jobs/led_control.py b/pioreactor/background_jobs/led_control.py index c7fbf287..096241d4 100644 --- a/pioreactor/background_jobs/led_control.py +++ b/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())}" ) diff --git a/pioreactor/background_jobs/monitor.py b/pioreactor/background_jobs/monitor.py index d0080eae..9fddf365 100644 --- a/pioreactor/background_jobs/monitor.py +++ b/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() diff --git a/pioreactor/background_jobs/od_reading.py b/pioreactor/background_jobs/od_reading.py index b6069f34..e705a8de 100644 --- a/pioreactor/background_jobs/od_reading.py +++ b/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 diff --git a/pioreactor/background_jobs/temperature_control.py b/pioreactor/background_jobs/temperature_control.py index 6b34d9ba..c1a48113 100644 --- a/pioreactor/background_jobs/temperature_control.py +++ b/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())}" ) diff --git a/pioreactor/structs.py b/pioreactor/structs.py index 9d3c157c..9e3d77e3 100644 --- a/pioreactor/structs.py +++ b/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 diff --git a/pioreactor/tests/test_od_reading.py b/pioreactor/tests/test_od_reading.py index 0744c5d5..ad146314 100644 --- a/pioreactor/tests/test_od_reading.py +++ b/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"] diff --git a/pioreactor/types.py b/pioreactor/types.py index af334a4f..75fe02ff 100644 --- a/pioreactor/types.py +++ b/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] diff --git a/requirements/requirements.txt b/requirements/requirements.txt index 3151ba72..d082e92d 100644 --- a/requirements/requirements.txt +++ b/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