calibration: expose diode calibration

changelog.md

@@ -13,6 +13,7 @@
 * Added `applied_at` property to a [calibration item](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html) obtained from a force slice. This property returns the timestamp in nanoseconds at which the force calibration was applied.
 * Added improved printing of calibration items under `channel.calibration` providing a more convenient overview of the items associated with a `Slice`.
 * Added improved printing of calibrations performed with `Pylake`.
+* Added property `diode_calibration` to access diode calibration model, and `trap_power` to access the used trap power in [calibration item](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html).
 * Added parameter `titles` to customize title of each subplot in [`Kymo.plot_with_channels()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.kymo.Kymo.html#lumicks.pylake.kymo.Kymo.plot_with_channels).
 * Added [`KymoTrack.sample_from_channel()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.kymotracker.kymotrack.KymoTrack.html#lumicks.pylake.kymotracker.kymotrack.KymoTrack.sample_from_channel) to downsample channel data to the time points of a kymotrack.
 

docs/api.rst

@@ -41,6 +41,7 @@ Force calibration
     ActiveCalibrationModel
     force_calibration.power_spectrum.PowerSpectrum
     force_calibration.power_spectrum_calibration.CalibrationResults
+    force_calibration.calibration_models.DiodeCalibrationModel
 
     :template: function.rst
 
@@ -157,4 +158,4 @@ Simulation
     :toctree: _api
     :template: function.rst
 
-    simulation.simulate_diffusive_tracks
+    simulation.simulate_diffusive_tracks

lumicks/pylake/force_calibration/calibration_item.py

@@ -2,6 +2,7 @@
 from functools import wraps
 from collections import UserDict
 
+from lumicks.pylake.force_calibration.calibration_models import DiodeCalibrationModel
 from lumicks.pylake.force_calibration.detail.calibration_properties import (
     CalibrationPropertiesMixin,
 )
@@ -37,6 +38,43 @@ def _fitted_diode(self):
         """Diode parameters were fitted"""
         return "f_diode (Hz)" in self or "alpha" in self
 
+    @property
+    def diode_calibration(self) -> DiodeCalibrationModel:
+        """Diode calibration model
+
+        The detector used to measure forces has a limited bandwidth. This bandwidth is typically
+        characterized by two parameters, a diode frequency and relaxation factor (which gives the
+        fraction of light that is instantaneously transmitted).
+
+        Some force detection diodes have been pre-calibrated in the factory. For these sensors,
+        this property will return a model that will allow you to calculate the diode parameters
+        at a particular trap sum power.
+
+        Examples
+        --------
+        ::
+
+            import lumicks.pylake as lk
+
+            f = lk.File("passive_calibration.h5")
+            item = f.force1x.calibration[0]  # Grab a calibration item for force 1x
+            diode_model = item.diode_calibration  # Grab diode model
+            diode_parameters = diode_model(item.trap_power)  # Grab diode parameters
+
+            # Verify that the calibration parameters at this trap power are the same as we found
+            # in the item in the first place.
+            assert diode_parameters["fixed_diode"] == item.diode_frequency
+            assert diode_parameters["fixed_alpha"] == item.diode_relaxation_factor
+        """
+        return DiodeCalibrationModel.from_calibration_dict(self.data)
+
+    @property
+    def trap_power(self):
+        """Average trap sum power in volts during calibration
+
+        Note that this property is only available for full calibrations with calibrated diodes."""
+        return self.data.get("Trap sum power (V)")
+
     @property
     def _sensor_type(self):
         if self.fast_sensor:

lumicks/pylake/force_calibration/calibration_models.py

@@ -1,7 +1,7 @@
 from copy import copy
 from typing import Callable
 from functools import partial
-from dataclasses import dataclass
+from dataclasses import field, dataclass
 
 import numpy as np
 
@@ -76,6 +76,53 @@ def diode_params_from_voltage(
     )
 
 
+@dataclass
+class DiodeCalibrationModel:
+    """Diode calibration model
+
+    This model takes a trap voltage and returns the diode parameters at that voltage. These
+    parameters can be passed directly to :func:`~lumicks.pylake.calibrate_force`.
+
+    Attributes
+    ----------
+    params : dict
+        Dictionary of diode parameters
+    """
+
+    _model_fun: Callable = field(repr=False)
+    params: dict
+
+    @staticmethod
+    def from_calibration_dict(calibration_dict):
+        try:
+            params = {
+                "delta_f_diode": calibration_dict["Diode frequency delta"],
+                "rate_f_diode": calibration_dict["Diode frequency rate"],
+                "max_f_diode": calibration_dict["Diode frequency max"],
+                "delta_alpha": calibration_dict["Diode alpha delta"],
+                "rate_alpha": calibration_dict["Diode alpha rate"],
+                "max_alpha": calibration_dict["Diode alpha max"],
+            }
+        except KeyError:
+            return None
+
+        def diode_fun(trap_voltage, params):
+            f_diode, alpha, _ = diode_params_from_voltage(trap_voltage, **params)
+            return {"fixed_diode": f_diode, "fixed_alpha": alpha}
+
+        return DiodeCalibrationModel(diode_fun, params)
+
+    def __call__(self, trap_voltage):
+        """Function to look up the diode parameters at a given trap power.
+
+        Parameters
+        ----------
+        trap_voltage : array_like
+            Array of trap voltages [V].
+        """
+        return self._model_fun(trap_voltage, self.params)
+
+
 def density_of_water(temperature, molarity, pressure=0.101325):
     """Determine the density of water with NaCl.
 

lumicks/pylake/force_calibration/tests/test_calibration_item.py

@@ -1,3 +1,5 @@
+import pickle
+
 import pytest
 
 from lumicks.pylake.calibration import ForceCalibrationList
@@ -143,8 +145,8 @@ def test_passive_item(compare_to_reference_dict, reference_data, calibration_dat
     assert item.number_of_samples == 781250
     assert not item.active_calibration
     assert not item.fast_sensor
-    assert item.start is 1696171376701856700
-    assert item.stop is 1696171386701856700
+    assert item.start == 1696171376701856700
+    assert item.stop == 1696171386701856700
     assert item.stiffness is ref_passive_fixed_diode_with_height["kappa (pN/nm)"]
     assert item.force_sensitivity is ref_passive_fixed_diode_with_height["Rf (pN/V)"]
     assert item.displacement_sensitivity is ref_passive_fixed_diode_with_height["Rd (um/V)"]
@@ -170,7 +172,7 @@ def test_passive_item(compare_to_reference_dict, reference_data, calibration_dat
     assert item.diffusion_volts_std_err == ref_passive_fixed_diode_with_height["err_D (V^2/s)"]
     assert not item.diode_frequency_std_err
     assert not item.diode_relaxation_factor_std_err
-    assert item.applied_at is 1696171386701856700
+    assert item.applied_at == 1696171386701856700
 
     compare_to_reference_dict(item.power_spectrum_params(), test_name="power")
     compare_to_reference_dict(item._model_params(), test_name="model")
@@ -194,8 +196,8 @@ def test_active_item_fixed_diode(compare_to_reference_dict, calibration_data):
     assert item.sample_rate == 78125
     assert item.active_calibration  # It is an active item!
     assert not item.fast_sensor
-    assert item.start is 1713785826919398000
-    assert item.stop is 1713785836900152600
+    assert item.start == 1713785826919398000
+    assert item.stop == 1713785836900152600
     assert item.stiffness is ref_active["kappa (pN/nm)"]
     assert item.force_sensitivity is ref_active["Rf (pN/V)"]
     assert item.displacement_sensitivity is ref_active["Rd (um/V)"]
@@ -264,9 +266,9 @@ def test_non_full(compare_to_reference_dict):
     )
     assert not item.stiffness
     assert not item.displacement_sensitivity
-    assert item.force_sensitivity is 1.0
-    assert item.start is 1714391268938540100
-    assert item.stop is 1714391268938540200
+    assert item.force_sensitivity == 1.0
+    assert item.start == 1714391268938540100
+    assert item.stop == 1714391268938540200
 
     for func in ("calibration_params", "power_spectrum_params"):
         with pytest.raises(
@@ -327,3 +329,39 @@ def create_item(t_start, t_stop, **kwargs):
         assert it == fc
 
     assert list(items) == items._src
+
+
+def test_item_pickle():
+    pickled_str = pickle.dumps(ForceCalibrationItem(ref_passive_fixed_diode_with_height))
+    loaded_object = pickle.loads(pickled_str)
+    assert loaded_object == ForceCalibrationItem(ref_passive_fixed_diode_with_height)
+
+
+@pytest.mark.parametrize(
+    "trap_power, reference_values",
+    [
+        (0, {"fixed_alpha": 0.45, "fixed_diode": 10000.0}),
+        (1000000, {"fixed_alpha": 0.5, "fixed_diode": 14000.0}),
+        (1.0, {"fixed_alpha": 0.48160602794142787, "fixed_diode": 12853.980812559239}),
+        ([0.5, 1.5], {"fixed_alpha": 0.48160602794142787, "fixed_diode": 12853.980812559239}),
+    ],
+)
+def test_diode_model(trap_power, reference_values):
+    item = ForceCalibrationItem(ref_passive_fixed_diode_with_height)
+    diode_calibration = item.diode_calibration
+    assert item.trap_power == ref_passive_fixed_diode_with_height["Trap sum power (V)"]
+    assert diode_calibration(trap_power) == reference_values
+    assert diode_calibration.params == {
+        "delta_f_diode": ref_passive_fixed_diode_with_height["Diode frequency delta"],
+        "rate_f_diode": ref_passive_fixed_diode_with_height["Diode frequency rate"],
+        "max_f_diode": ref_passive_fixed_diode_with_height["Diode frequency max"],
+        "delta_alpha": ref_passive_fixed_diode_with_height["Diode alpha delta"],
+        "rate_alpha": ref_passive_fixed_diode_with_height["Diode alpha rate"],
+        "max_alpha": ref_passive_fixed_diode_with_height["Diode alpha max"],
+    }
+
+
+def test_no_diode_model():
+    item = ForceCalibrationItem(ref_active)
+    assert item.trap_power is None
+    assert item.diode_calibration is None