channel: support recalibrating the force

changelog.md

@@ -6,6 +6,7 @@
 
 * Added model to correct for bead-bead coupling when using active calibration deep in bulk with two beads. See [`tutorial`](https://lumicks-pylake.readthedocs.io/en/latest/tutorial/force_calibration.html#active-calibration-with-two-beads-far-away-from-the-surface) and [`theory`](https://lumicks-pylake.readthedocs.io/en/latest/theory/force_calibration/active.html#bead-bead-coupling) for more information.
 * Added option to highlight a region on a time plot using [`Slice.highlight_time_range()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.channel.Slice.html#lumicks.pylake.channel.Slice.highlight_time_range). For more information see the [`tutorial`](https://lumicks-pylake.readthedocs.io/en/latest/tutorial/file.html#highlight-time-range).
+* Support recalibrating force data using [`Slice.recalibrate_force()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.channel.Slice.html#lumicks.pylake.channel.Slice.recalibrate_force).
 * Added `__array__` to [`Slice`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.channel.Slice.html). This allows passing slices directly to `numpy` functions such as `np.mean()`or `np.sum()`.
 * Added parameter `allow_overwrite` to [`lk.download_from_doi()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.download_from_doi.html#lumicks.pylake.download_from_doi) to allow re-downloading only those files where the checksum does not match.
 * Added force calibration information to channels accessed directly via the square bracket notation (e.g. `file["Force HF"]["Force 1x"].calibration`).

lumicks/pylake/channel.py

@@ -215,6 +215,66 @@ def calibration(self) -> list:
         else:
             return []
 
+    def recalibrate_force(self, calibration):
+        """Recalibrate the data in this slice using a new force calibration item
+
+        Parameters
+        ----------
+        calibration : ForceCalibrationItem
+            The calibration to apply to the data
+
+        Examples
+        --------
+        ::
+
+            import lumicks.pylake as lk
+
+            f = lk.File("passive_calibration.h5")
+            calibration = f.force1x.calibration[1]  # Grab a calibration item for force 1x
+
+            # Slice the data corresponding to the calibration we want to reproduce.
+            calib_slice = f.force1x[calibration]
+
+            # De-calibrate to volts using the calibration that was active before this slice.
+            previous_calibration = calib_slice.calibration[0]
+            calib_slice = calib_slice / previous_calibration.force_sensitivity
+
+            calibration_params = previous_calibration.calibration_params()
+            new_calibration = lk.calibrate_force(calib_slice.data, **calibration_params)
+            new_calibration.plot()
+
+            # Make a new calibration, but change the amount of blocking
+            less_blocking_params = calibration_params | {"num_points_per_block": 200}
+            less_blocking = lk.calibrate_force(calib_slice.data, **less_blocking_params)
+            less_blocking.plot()
+
+            # Recalibrate the force channels
+            recalibrated_force1x = f.force1x.recalibrate_force(less_blocking)
+        """
+        from lumicks.pylake.calibration import ForceCalibrationList
+
+        if len(self.calibration) > 1:
+            raise NotImplementedError(
+                "Slice contains multiple calibrations. Recalibration is only implemented for slices"
+                "that do not contain calibration events."
+            )
+
+        if not self.calibration:
+            raise RuntimeError(
+                "Slice does not contain any calibration items. Is this an unprocessed force "
+                "channel?"
+            )
+
+        return Slice(
+            self._src._with_data(
+                self._unpack_other(self)
+                * calibration.force_sensitivity
+                / self.calibration[0].force_sensitivity
+            ),
+            labels=self.labels,
+            calibration=ForceCalibrationList(items=[calibration._with_timestamp(self.start)]),
+        )
+
     @property
     def sample_rate(self) -> Union[float, None]:
         """The data frequency for `Continuous` and `TimeSeries` data sources or `None` if it is not

lumicks/pylake/conftest.py

@@ -2,6 +2,7 @@
 import hashlib
 import warnings
 import importlib
+import contextlib
 
 import numpy as np
 import pytest
@@ -70,6 +71,27 @@ def report():
     return reporter
 
 
+@pytest.fixture()
+def mock_datetime(monkeypatch):
+    class FakeDateTime:
+        @classmethod
+        def fromtimestamp(cls, timestamp, *args, **kwargs):
+            """Inject a timezone"""
+            return FakeDateTime()
+
+        def strftime(self, str_format):
+            return str_format
+
+    @contextlib.contextmanager
+    def fake_datetime(function_location):
+        # Representation of time is timezone and locale dependent, hence we monkeypatch it
+        with monkeypatch.context() as m:
+            m.setattr(function_location, FakeDateTime)
+            yield m
+
+    return fake_datetime
+
+
 @pytest.fixture(autouse=True)
 def configure_warnings():
     # importing scipy submodules on some version of Python

lumicks/pylake/force_calibration/calibration_item.py

@@ -186,8 +186,7 @@ def calibration_params(self):
             less_blocking.plot()
 
             # Recalibrate the force channels
-            rf_ratio = less_blocking.force_sensitivity / previous_calibration.force_sensitivity
-            recalibrated_force1x = f.force1x * rf_ratio
+            recalibrated_force1x = f.force1x.recalibrate_force(less_blocking)
         """
         return (
             self.power_spectrum_params()

lumicks/pylake/force_calibration/calibration_results.py

@@ -52,9 +52,13 @@ def __call__(self, frequency):
 
     def _with_timestamp(self, applied_timestamp):
         """Return a copy of this item with a timestamp of when it was applied"""
+        from lumicks.pylake.force_calibration.power_spectrum_calibration import CalibrationParameter
+
         item = copy.copy(self)
         item.params = copy.deepcopy(self.params)
-        item.params["Timestamp (ns)"] = applied_timestamp
+        item.params["Timestamp"] = CalibrationParameter(
+            "Timestamp when item was applied", applied_timestamp, "nanoseconds"
+        )
         return item
 
     def __contains__(self, key):

lumicks/pylake/tests/test_channels/test_channels.py

@@ -1,4 +1,5 @@
 import re
+import textwrap
 from collections import namedtuple
 from dataclasses import dataclass
 
@@ -10,6 +11,7 @@
 from lumicks.pylake import channel
 from lumicks.pylake.low_level import make_continuous_slice
 from lumicks.pylake.calibration import ForceCalibrationList
+from lumicks.pylake.force_calibration import power_spectrum_calibration as psc
 
 
 def with_offset(t, start_time=1592916040906356300):
@@ -1002,3 +1004,68 @@ def test_low_level_construction():
     slc = make_continuous_slice(data, start, int(1e9 / 78125), name="hi", y_label="there")
     assert slc.labels["title"] == "hi"
     assert slc.labels["y"] == "there"
+
+
+def test_recalibrate_force_wrong_number_of_calibrations():
+    too_many = ForceCalibrationList._from_dicts(
+        items=[
+            {"Calibration Data": 50, "Stop time (ns)": with_offset(50)},
+            {"Calibration Data": 80, "Stop time (ns)": with_offset(80)},
+        ],
+    )
+
+    cc = channel.Slice(
+        channel.Continuous(np.arange(100), int(with_offset(40)), 10), calibration=too_many
+    )
+    with pytest.raises(NotImplementedError, match="Slice contains multiple calibrations"):
+        cc.recalibrate_force(None)
+
+    cc = channel.Slice(
+        channel.Continuous(np.arange(100), int(with_offset(40)), 10),
+        calibration=ForceCalibrationList([]),
+    )
+    with pytest.raises(RuntimeError, match="Slice does not contain any calibration items"):
+        cc.recalibrate_force(None)
+
+
+def test_recalibrate_force(mock_datetime):
+    slc = channel.Slice(
+        channel.Continuous(np.arange(100), with_offset(40), 10),
+        calibration=ForceCalibrationList._from_dicts(
+            [{"Calibration Data": 50, "Stop time (ns)": with_offset(50), "Response (pN/V)": 10}],
+        ),
+    )
+
+    def make_item(calibration_factor):
+        return psc.CalibrationResults(
+            model=None,
+            ps_model=None,
+            ps_data=None,
+            params={},
+            results={
+                "Rf": psc.CalibrationParameter("Rf", calibration_factor, "val"),
+                "kappa": psc.CalibrationParameter("kappa", 5, "val"),
+            },
+            fitted_params=[],
+        )
+
+    slc_recalibrated = slc.recalibrate_force(make_item(5))
+    np.testing.assert_allclose(slc.data, np.arange(100))
+    np.testing.assert_allclose(slc_recalibrated.data, np.arange(100) * 0.5)
+    assert slc.calibration[0].force_sensitivity == 10
+    assert slc_recalibrated.calibration[0].force_sensitivity == 5
+
+    with mock_datetime("lumicks.pylake.calibration.datetime.datetime"):
+        assert str(slc_recalibrated.calibration) == textwrap.dedent(
+            """\
+              #  Applied at    Kind       Stiffness (pN/nm)    Force sens. (pN/V)  Disp. sens. (µm/V)    Hydro    Surface    Data?
+            ---  ------------  -------  -------------------  --------------------  --------------------  -------  ---------  -------
+              0  %x %X         Unknown                    5                     5  N/A                   False    False      False"""
+        )
+    slc_recalibrated.calibration._repr_html_()
+
+    slc_recalibrate_twice = slc.recalibrate_force(make_item(10))
+    np.testing.assert_allclose(slc_recalibrated.data, np.arange(100) * 0.5)
+    np.testing.assert_allclose(slc_recalibrate_twice.data, np.arange(100))
+    assert slc_recalibrated.calibration[0].force_sensitivity == 5
+    assert slc_recalibrate_twice.calibration[0].force_sensitivity == 10

lumicks/pylake/tests/test_file/test_file.py

@@ -75,25 +75,10 @@ def test_redirect_list(h5_file):
         assert f["Point Scan"]["PointScan1"].start == np.int64(20e9)
 
 
-def test_calibration_str(h5_file, monkeypatch):
-    # Representation of time is timezone and locale dependent, hence we monkeypatch it
-    class FakeDateTime:
-        @classmethod
-        def fromtimestamp(cls, timestamp, *args, **kwargs):
-            """Inject a timezone"""
-            return FakeDateTime()
-
-        def strftime(self, str_format):
-            return str_format
-
+def test_calibration_str(h5_file, mock_datetime):
     f = pylake.File.from_h5py(h5_file)
     if f.format_version == 2:
-        with monkeypatch.context() as m:
-            m.setattr(
-                "lumicks.pylake.calibration.datetime.datetime",
-                FakeDateTime,
-            )
-
+        with mock_datetime("lumicks.pylake.calibration.datetime.datetime"):
             assert str(f.force1x.calibration) == dedent(
                 (
                     """\