From d35aa1132cd6300ca97e61151bd6ea6a9fe7d4a7 Mon Sep 17 00:00:00 2001 From: Joep Vanlier Date: Mon, 20 Jan 2025 18:32:16 +0100 Subject: [PATCH] calibration: add convenience api for items that have the calibration inside, you can make it even easier for users --- changelog.md | 1 + .../force_calibration/calibration_items.rst | 67 +++++++++++++++++++ .../figures/in_item_less_blocking.png | 3 + .../figures/in_item_no_hyco.png | 3 + .../figures/in_item_no_hyco_trace.png | 3 + .../figures/in_item_plot.png | 3 + .../figures/in_item_raw_data.png | 3 + lumicks/pylake/calibration.py | 2 +- .../force_calibration/calibration_item.py | 43 +++++++++++- .../tests/test_calibration_item.py | 43 ++++++++++++ .../pylake/tests/test_file/test_file_items.py | 26 +++++++ 11 files changed, 193 insertions(+), 4 deletions(-) create mode 100644 docs/tutorial/force_calibration/figures/in_item_less_blocking.png create mode 100644 docs/tutorial/force_calibration/figures/in_item_no_hyco.png create mode 100644 docs/tutorial/force_calibration/figures/in_item_no_hyco_trace.png create mode 100644 docs/tutorial/force_calibration/figures/in_item_plot.png create mode 100644 docs/tutorial/force_calibration/figures/in_item_raw_data.png diff --git a/changelog.md b/changelog.md index 979b6fb10..530c4eb51 100644 --- a/changelog.md +++ b/changelog.md @@ -10,6 +10,7 @@ * Force `num_points_per_block` to be integer when down-sampling a [`PowerSpectrum`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.force_calibration.power_spectrum.PowerSpectrum.html). * Added [voltage](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html#lumicks.pylake.calibration.ForceCalibrationItem.voltage), [sum_voltage](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html#lumicks.pylake.calibration.ForceCalibrationItem.sum_voltage) and [driving](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html#lumicks.pylake.calibration.ForceCalibrationItem.driving) properties which return raw calibration data. Note that these are only available when the option to export the raw data has explicitly been selected in Bluelake. When unavailable, these properties will return empty slices. +* Added [ForceCalibrationItem.recalibrate_with()](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html#lumicks.pylake.calibration.ForceCalibrationItem.recalibrate_with) and [ForceCalibrationItem.plot()](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.calibration.ForceCalibrationItem.html#lumicks.pylake.calibration.ForceCalibrationItem.plot) to recalibrate and/or plot calibrations performed in Bluelake. Note that the raw data needs to have been exported inside the item (Bluelake 2.7.x feature) for this feature to work. ## v1.6.1 | 2025-01-29 diff --git a/docs/tutorial/force_calibration/calibration_items.rst b/docs/tutorial/force_calibration/calibration_items.rst index 2e5ed54ea..c137d5fa0 100644 --- a/docs/tutorial/force_calibration/calibration_items.rst +++ b/docs/tutorial/force_calibration/calibration_items.rst @@ -69,11 +69,78 @@ These parameters are properties and can be extracted as such:: Redoing a Bluelake calibration ------------------------------ +Starting from Bluelake `2.7.0`, it is possible to export the raw data used for calibration with the +calibration item. You can enable this in the settings panel in Bluelake. +With this feature, recalibrating your data differently becomes a lot easier. +If you don't have this setting enabled, you can still recalibrate your data, but you will have to +ensure that you manually export the raw calibration data and slice the appropriate calibration data. +To find out how to do this, please refer to the section :ref:`Recalibrating using timeline data`. + +.. _recalibrating_simple: + +Recalibrating using calibration items with raw data +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Let's start by loading our calibration item:: + + f = lk.File("test_data/raw_data_in_item.h5") + calibration = f.force1x.calibration[0] + +We can quickly plot this calibration:: + + calibration.plot() + +.. image:: figures/in_item_plot.png + +If you wish to obtain the raw calibration data, you can simply access the properties :meth:`~lumicks.pylake.force_calibration.calibration_item.ForceCalibrationItem.calibration.voltage` and for active calibration :meth:`~lumicks.pylake.force_calibration.calibration_item.ForceCalibrationItem.calibration.driving`. +These return slices you can plot and interact with:: + + calibration.voltage.plot() + +.. image:: figures/in_item_raw_data.png + +We can easily re-perform this calibration by invoking :meth:`~lumicks.pylake.force_calibration.calibration_item.ForceCalibrationItem.recalibrate_with()`. +Let's see what this spectrum would have looked like with less blocking:: + + recalibrated = calibration.recalibrate_with(num_points_per_block=200) + recalibrated.plot() + +.. image:: figures/in_item_less_blocking.png + +Note that any of the calibration parameters can easily be changed this way. +To investigate the effect of the hydrodynamically correct model for example, we can try turning it off:: + + recalibrated_no_hyco = calibration.recalibrate_with(hydrodynamically_correct=False) + recalibrated_no_hyco.plot() + +.. image:: figures/in_item_no_hyco.png + +This clearly fits the data poorly. +To see what it would do to our timeline data, we can simply recalibrate that by applying the calibration item to :meth:`~lumicks.pylake.channel.Slice.recalibrate_force`:: + + recalibrated_force = f.force1x.recalibrate_force(recalibrated_no_hyco) + + f.force1x.plot() + recalibrated_force.plot() + +.. image:: figures/in_item_no_hyco_trace.png + +That's all there is to it. + +Recalibrating using timeline data +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. _recalibrating_manually: + .. important:: In order to redo a Bluelake calibration, the force data that was used for the calibration has to be included in the `.h5` file. *Note that this force data is not exported nor marked by default*; it has to be explicitly added to the exported file. +We start by loading the calibration item:: + + f = lk.File("test_data/passive_calibration.h5") + We can directly slice the channel by the calibration item we want to reproduce to extract the relevant data:: force1x_slice = f.force1x[f.force1x.calibration[1]] diff --git a/docs/tutorial/force_calibration/figures/in_item_less_blocking.png b/docs/tutorial/force_calibration/figures/in_item_less_blocking.png new file mode 100644 index 000000000..585a02a3d --- /dev/null +++ b/docs/tutorial/force_calibration/figures/in_item_less_blocking.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3a50b0bc2c28c89cee7c8edc4e1008eaceef7a0eb11d8922c303890a31283712 +size 35090 diff --git a/docs/tutorial/force_calibration/figures/in_item_no_hyco.png b/docs/tutorial/force_calibration/figures/in_item_no_hyco.png new file mode 100644 index 000000000..3fb5aed43 --- /dev/null +++ b/docs/tutorial/force_calibration/figures/in_item_no_hyco.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:448976ec8bc432b3549bc41b9edaf001c61cf2d9eee0c153c50bbbe56eed507f +size 27785 diff --git a/docs/tutorial/force_calibration/figures/in_item_no_hyco_trace.png b/docs/tutorial/force_calibration/figures/in_item_no_hyco_trace.png new file mode 100644 index 000000000..ea1766ce7 --- /dev/null +++ b/docs/tutorial/force_calibration/figures/in_item_no_hyco_trace.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c1e0c090e8ebc5fb8e1c22f1fe011dd69a556a425386ca2e2d63d9fa249f3f16 +size 21989 diff --git a/docs/tutorial/force_calibration/figures/in_item_plot.png b/docs/tutorial/force_calibration/figures/in_item_plot.png new file mode 100644 index 000000000..ce069f7f2 --- /dev/null +++ b/docs/tutorial/force_calibration/figures/in_item_plot.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c4aebd892979f1eb5399a53b62292db3e4ad170e3faaf1062cd0d976eb499b77 +size 26165 diff --git a/docs/tutorial/force_calibration/figures/in_item_raw_data.png b/docs/tutorial/force_calibration/figures/in_item_raw_data.png new file mode 100644 index 000000000..3d93b03ae --- /dev/null +++ b/docs/tutorial/force_calibration/figures/in_item_raw_data.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d2ddeceafd06752f97dd0b2caac605907d3c88dd6501b7913f7f544480dae66f +size 31193 diff --git a/lumicks/pylake/calibration.py b/lumicks/pylake/calibration.py index d086924ec..3c058218b 100644 --- a/lumicks/pylake/calibration.py +++ b/lumicks/pylake/calibration.py @@ -103,7 +103,7 @@ def from_field(hdf5, force_channel) -> "ForceCalibrationList": def make_slice(dset, field, y_label, title) -> Slice: """Fetch raw data from the dataset""" - if field in dset: + if field in dset and dset[field].size > 0: return Slice( Continuous.from_dataset(dset[field]), labels={"x": "Time (s)", "y": y_label, "title": title}, diff --git a/lumicks/pylake/force_calibration/calibration_item.py b/lumicks/pylake/force_calibration/calibration_item.py index 9e3dd4842..f79b9cee6 100644 --- a/lumicks/pylake/force_calibration/calibration_item.py +++ b/lumicks/pylake/force_calibration/calibration_item.py @@ -28,17 +28,17 @@ def __init__( @property def voltage(self): """Uncalibrated voltage reading on the detector""" - return self._voltage if self._voltage else empty_slice + return self._voltage if self._voltage is not None else empty_slice @property def sum_voltage(self): """Uncalibrated sum voltage on the detector""" - return self._sum_voltage if self._sum_voltage else empty_slice + return self._sum_voltage if self._sum_voltage is not None else empty_slice @property def driving(self): """Driving signal used for active calibration""" - return self._driving if self._driving else empty_slice + return self._driving if self._driving is not None else empty_slice @staticmethod def _verify_full(method): @@ -160,6 +160,43 @@ def __repr__(self): ) return f"{self.__class__.__name__}({properties})" + def _check_has_data(self): + if not len(self.voltage): + raise ValueError( + "This calibration item does not contain the raw data. If you still have the " + "timeline force data, you can de-calibrate that and perform the re-calibration" + "manually. See the pylake tutorial on force calibration for more information." + ) + + def plot(self): + self._check_has_data() + self.recalibrate_with().plot() + + def plot_spectrum_residual(self): + """Plot the residuals of the fitted spectrum. + + This diagnostic plot can be used to determine how well the spectrum fits the data. While + it cannot be used to diagnose over-fitting (being unable to reliably estimate parameters + due to insufficient information in the data), it can be used to diagnose under-fitting (the + model not fitting the data adequately). + + In an ideal situation, the residual plot should show a noise band around 1 without any + systematic deviations. + """ + self._check_has_data() + self.recalibrate_with().plot_spectrum_residual() + + def recalibrate_with(self, **params): + """Returns a calibration structure with some parameters overridden. + + For a full list of parameters to override, please see + :func:`~lumicks.pylake.calibrate_force()`""" + self._check_has_data() + active_data = {"driving_data": self.driving.data} if self.active_calibration else {} + return calibrate_force( + self.voltage.data, **(self.calibration_params() | active_data | params) + ) + @_verify_full def _model_params(self): """Returns parameters with which to create an active or passive calibration model""" diff --git a/lumicks/pylake/force_calibration/tests/test_calibration_item.py b/lumicks/pylake/force_calibration/tests/test_calibration_item.py index a7a254287..826602ceb 100644 --- a/lumicks/pylake/force_calibration/tests/test_calibration_item.py +++ b/lumicks/pylake/force_calibration/tests/test_calibration_item.py @@ -297,6 +297,49 @@ def test_non_full(compare_to_reference_dict): ) +def test_plot_item(active_ref_data): + item = ForceCalibrationItem(ref_active) + with pytest.raises(ValueError, match="This calibration item does not contain the raw data"): + item.plot() + + with pytest.raises(ValueError, match="This calibration item does not contain the raw data"): + item.plot_spectrum_residual() + + with pytest.raises(ValueError, match="This calibration item does not contain the raw data"): + item.recalibrate_with(hydrodynamically_correct=False) + + voltage, driving = active_ref_data + + item = ForceCalibrationItem(ref_active, voltage=voltage) + with pytest.raises( + ValueError, match="Active calibration requires the driving_data to be defined" + ): + item.plot() + + item = ForceCalibrationItem(ref_active, voltage=voltage, driving=driving) + item.plot() + item.plot_spectrum_residual() + + +def test_recalibrate_item(active_ref_data): + voltage, driving = active_ref_data + item = ForceCalibrationItem(ref_active, voltage=voltage, driving=driving) + np.testing.assert_allclose(item.stiffness, 0.50338, rtol=1e-4) + + same_calibration = item.recalibrate_with() + np.testing.assert_allclose(same_calibration.stiffness, 0.521374, rtol=1e-4) + assert same_calibration.fitted_diode + + recalibrated = item.recalibrate_with( + hydrodynamically_correct=False, + fixed_diode=17498.229, + fixed_alpha=0.135, + distance_to_surface=1.5, + ) + np.testing.assert_allclose(recalibrated.stiffness, 0.498051, rtol=1e-4) + assert not recalibrated.fitted_diode + + def test_force_calibration_handling(): def timestamp(time): return 1714391268938540100 + int(1e9) * time diff --git a/lumicks/pylake/tests/test_file/test_file_items.py b/lumicks/pylake/tests/test_file/test_file_items.py index f406c0414..93a46798c 100644 --- a/lumicks/pylake/tests/test_file/test_file_items.py +++ b/lumicks/pylake/tests/test_file/test_file_items.py @@ -89,6 +89,32 @@ def test_calibration(h5_file): assert f.force1y.calibration[idx].driving.labels["title"] == "Driving data for axis y" +def test_calibration_lazy(monkeypatch, h5_file): + import h5py + + f = pylake.File.from_h5py(h5_file) + if f.format_version == 2: + data_accessed = False + h5_array = h5py.Dataset.__array__ + + def patched_array(self, *args, **kwargs): + nonlocal data_accessed + data_accessed = True + return h5_array(self, *args, **kwargs) + + with monkeypatch.context() as m: + m.setattr("h5py.Dataset.__array__", patched_array) + + assert not f.force1x.calibration[2].has_data + assert f.force1y.calibration[2].voltage.sample_rate == 78125 + assert not data_accessed # Verify that we didn't read the slice yet + + # Read the data + _ = f.force1y.calibration[2].voltage.data + + assert data_accessed + + def test_marker(h5_file): f = pylake.File.from_h5py(h5_file)