Skip to content

Commit

Permalink
fix calibibrations near 0
Browse files Browse the repository at this point in the history
  • Loading branch information
CamDavidsonPilon committed Jan 15, 2025
1 parent 200b9d6 commit ae6625e
Show file tree
Hide file tree
Showing 5 changed files with 43 additions and 3 deletions.
1 change: 1 addition & 0 deletions pioreactor/background_jobs/od_reading.py
Original file line number Diff line number Diff line change
Expand Up @@ -1260,3 +1260,4 @@ def click_od_reading(
calibration=possible_calibration,
)
od.block_until_disconnected()

6 changes: 4 additions & 2 deletions pioreactor/calibrations/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,14 @@ def bold(string: str) -> str:
def calculate_poly_curve_of_best_fit(x: list[float], y: list[float], degree: int) -> list[float]:
import numpy as np

# weigh the last point, the "blank measurement", more.
# weigh the smallest point, the "blank measurement", more.
# 1. It's far away from the other points
# 2. We have prior knowledge that OD~0 when V~0.
n = len(x)
weights = np.ones_like(x)
weights[-1] = n / 2
weights[0] = n / 2

x, y = zip(*sorted(zip(x, y), key=lambda t: t[0]))

try:
coefs = np.polyfit(x, y, deg=degree, w=weights)
Expand Down
2 changes: 1 addition & 1 deletion pioreactor/structs.py
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ def ipredict(self, y: Y, enforce_bounds=False) -> X:
coef_shift = zeros_like(poly)
coef_shift[-1] = y
solve_for_poly = poly - coef_shift
roots_ = roots(solve_for_poly)
roots_ = roots(solve_for_poly).tolist()
plausible_sols_: list[X] = sorted([real(r) for r in roots_ if (abs(imag(r)) < 1e-10)])

if len(plausible_sols_) == 0:
Expand Down
31 changes: 31 additions & 0 deletions pioreactor/tests/test_od_calibration.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@

from pioreactor.calibrations.utils import calculate_poly_curve_of_best_fit
from pioreactor.calibrations.utils import curve_to_callable
from pioreactor.structs import ODCalibration
from pioreactor.utils.timing import current_utc_datetime


def test_linear_data_produces_linear_curve_in_range_even_if_high_degree() -> None:
Expand All @@ -24,3 +26,32 @@ def test_linear_data_produces_linear_curve_in_range_even_if_high_degree() -> Non
curve_callable = curve_to_callable("poly", curve_data_)
for od_, v_ in zip(od, curve_callable(od)):
assert (v_ - od_ * 0.5) < 0.035


def test_mandys_data_for_pathological_poly() -> None:

od = [0.0, 0.139, 0.155, 0.378, 0.671, 0.993, 1.82, 4.061]
v = [0.0, 0.0158, 0.0322, 0.0589, 0.1002, 0.1648, 0.4045, 0.5463]

curve_data_ = calculate_poly_curve_of_best_fit(od, v, degree=3) # type: ignore
curve_callable = curve_to_callable("poly", curve_data_)
assert abs(curve_callable(0.002) - 0.002) < 0.1

mcal = ODCalibration(
calibration_name='mandy',
calibrated_on_pioreactor_unit='pio1',
created_at=current_utc_datetime(),
curve_data_=curve_data_,
curve_type='poly',
recorded_data={'x': od, 'y': v},
ir_led_intensity=70.0,
angle='90',
pd_channel='2')

assert abs(mcal.predict(0.002) - curve_callable(0.002)) < 1e-10
assert abs(mcal.ipredict(0.002) - 0.002) < 0.1





6 changes: 6 additions & 0 deletions pioreactor/tests/test_od_reading.py
Original file line number Diff line number Diff line change
Expand Up @@ -1120,3 +1120,9 @@ def test_missing_calibration_data():
cal_transformer = CachedCalibrationTransformer()
cal_transformer.models = {"1": lambda v: v * 2}
assert cal_transformer({"1": 1.0, "2": 2.0}) == {"1": 2.0, "2": 2.0}


def test_mandys_calibration():
from pioreactor.calibrations import load_calibration
mcal = structs.ODCalibration(calibration_name='mandy', calibrated_on_pioreactor_unit='pio1', created_at=current_utc_datetime(), curve_data_=[-0.03112259838616315, 0.14606367297714123, 0.05224678328234911, 0.009665339167023364], curve_type='poly', x='voltage', y='od600s', recorded_data={'x': [0.0, 0.139, 0.155, 0.378, 0.671, 0.993, 1.82, 4.061], 'y': [0.0, 0.0158, 0.0322, 0.0589, 0.1002, 0.1648, 0.4045, 0.5463]}, ir_led_intensity=70.0, angle='90', pd_channel='2')
assert 0.0 < mcal.ipredict(0.002, enforce_bounds=True) < 1.0

0 comments on commit ae6625e

Please sign in to comment.