From beea00088c2e1d79b0e3c426dccdb6422e361d08 Mon Sep 17 00:00:00 2001 From: Jack Harper Date: Wed, 6 Nov 2024 13:18:28 +0000 Subject: [PATCH 01/16] backing up what we have so far with bluesky prs --- manual_system_tests/dae_scan.py | 60 +++++++++++-------- .../callbacks/fitting/__init__.py | 4 +- .../callbacks/fitting/fitting_utils.py | 2 + src/ibex_bluesky_core/callbacks/plotting.py | 5 ++ .../devices/dae/dae_spectra.py | 5 +- 5 files changed, 49 insertions(+), 27 deletions(-) diff --git a/manual_system_tests/dae_scan.py b/manual_system_tests/dae_scan.py index f4d12cc..7b9f20c 100644 --- a/manual_system_tests/dae_scan.py +++ b/manual_system_tests/dae_scan.py @@ -8,15 +8,17 @@ import bluesky.plans as bp import matplotlib import matplotlib.pyplot as plt -from bluesky.callbacks import LiveTable +from bluesky.callbacks import LiveTable, LiveFitPlot from bluesky.preprocessors import subs_decorator from bluesky.utils import Msg from ophyd_async.plan_stubs import ensure_connected from ibex_bluesky_core.callbacks.file_logger import HumanReadableFileCallback +from ibex_bluesky_core.callbacks.fitting import LiveFit +from ibex_bluesky_core.callbacks.fitting.fitting_utils import Linear, Gaussian from ibex_bluesky_core.callbacks.plotting import LivePlot from ibex_bluesky_core.devices import get_pv_prefix -from ibex_bluesky_core.devices.block import block_rw_rbv +from ibex_bluesky_core.devices.block import block_rw_rbv, BlockWriteConfig from ibex_bluesky_core.devices.simpledae import SimpleDae from ibex_bluesky_core.devices.simpledae.controllers import ( RunPerPointController, @@ -47,10 +49,12 @@ def dae_scan_plan() -> Generator[Msg, None, None]: - The DAE waited for at least 500 good frames at each point """ prefix = get_pv_prefix() - block = block_rw_rbv(float, "mot") + block = block_rw_rbv(float, "mot", write_config=BlockWriteConfig(settle_time_s=0.5)) + blocka = block_rw_rbv(float, "alice") + blockb = block_rw_rbv(float, "bob", write_config=BlockWriteConfig(settle_time_s=0.5)) controller = RunPerPointController(save_run=True) - waiter = GoodFramesWaiter(500) + waiter = GoodFramesWaiter(1) reducer = GoodFramesNormalizer( prefix=prefix, detector_spectra=[i for i in range(1, 100)], @@ -67,38 +71,44 @@ def dae_scan_plan() -> Generator[Msg, None, None]: controller.run_number.set_name("run number") reducer.intensity.set_name("normalized counts") - yield from ensure_connected(block, dae, force_reconnect=True) + yield from ensure_connected(block, blocka, blockb, force_reconnect=True) + print(reducer.intensity.name) + + + + lf = LiveFit(Gaussian.fit(), y=blockb.name, x=block.name, yerr=blocka.name) + fig, ax = plt.subplots() + @subs_decorator( [ - HumanReadableFileCallback( - Path("C:\\") / "instrument" / "var" / "logs" / "bluesky" / "output_files", - [ - block.name, - controller.run_number.name, - reducer.intensity.name, - reducer.det_counts.name, - dae.good_frames.name, - ], - ), - LivePlot(y=reducer.intensity.name, x=block.name, marker="x", linestyle="none"), + # HumanReadableFileCallback( + # Path("C:\\") / "instrument" / "var" / "logs" / "bluesky" / "output_files", + # [ + # block.name, + # controller.run_number.name, + # reducer.intensity.name, + # reducer.det_counts.name, + # dae.good_frames.name, + # ], + # ), + LiveFitPlot(lf, ax=ax, color="r", num_points=1000), + LivePlot(y=blockb.name, x=block.name, yerr=blocka.name, marker="x", linestyle="none", ax=ax), LiveTable( [ block.name, - controller.run_number.name, - reducer.intensity.name, - reducer.intensity_stddev.name, - reducer.det_counts.name, - reducer.det_counts_stddev.name, - dae.good_frames.name, + blocka.name, + blockb.name ] ), + ] ) def _inner() -> Generator[Msg, None, None]: - num_points = 3 - yield from bps.mv(dae.number_of_periods, num_points) - yield from bp.scan([dae], block, 0, 10, num=num_points) + num_points = 10 + # yield from bps.mv(dae.number_of_periods, num_points) + yield from bps.read(blocka) + yield from bp.scan([blockb, blocka], block, 0, 10, num=num_points) yield from _inner() diff --git a/src/ibex_bluesky_core/callbacks/fitting/__init__.py b/src/ibex_bluesky_core/callbacks/fitting/__init__.py index 6ef3bd9..e3ad034 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/__init__.py +++ b/src/ibex_bluesky_core/callbacks/fitting/__init__.py @@ -56,6 +56,7 @@ def __init__( x: str, *, update_every: int = 1, + yerr = None ) -> None: """Call Bluesky LiveFit with assumption that there is only one independant variable. @@ -73,6 +74,7 @@ def __init__( y=y, independent_vars={"x": x}, update_every=update_every, + yerr=yerr ) def update_fit(self) -> None: @@ -91,4 +93,4 @@ def update_fit(self) -> None: # Calls the guess function on the set of data already collected in the run ) - super().update_fit() + super().update_fit() \ No newline at end of file diff --git a/src/ibex_bluesky_core/callbacks/fitting/fitting_utils.py b/src/ibex_bluesky_core/callbacks/fitting/fitting_utils.py index 0fc0957..88cda45 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/fitting_utils.py +++ b/src/ibex_bluesky_core/callbacks/fitting/fitting_utils.py @@ -197,6 +197,8 @@ def guess( numerator = sum(x * y) - sum(x) * sum(y) denominator = sum(x**2) - sum(x) ** 2 + print(f"N: {numerator} D: {denominator}") + m = numerator / denominator c = (sum(y) - m * sum(x)) / len(x) diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index 4299b27..a4f2119 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -15,6 +15,10 @@ class LivePlot(_DefaultLivePlot): """Live plot, customized for IBEX.""" + def __init__(self, y, x=None, yerr=None, *args, **kwargs): + super().__init__(y=y, x=x, yerr=yerr, *args, **kwargs) + + def _show_plot(self) -> None: # Play nicely with the "normal" backends too - only force show if we're # actually using our custom backend. @@ -28,5 +32,6 @@ def start(self, doc: RunStart) -> None: def event(self, doc: Event) -> None: """Process an event document (delegate to superclass, then show the plot).""" + super().event(doc) self._show_plot() diff --git a/src/ibex_bluesky_core/devices/dae/dae_spectra.py b/src/ibex_bluesky_core/devices/dae/dae_spectra.py index 0c90469..e0320e2 100644 --- a/src/ibex_bluesky_core/devices/dae/dae_spectra.py +++ b/src/ibex_bluesky_core/devices/dae/dae_spectra.py @@ -107,7 +107,10 @@ async def read_spectrum_dataarray(self) -> sc.DataArray: if unit is None: raise ValueError("Could not determine engineering units of tof edges.") + # TODO add reference to ADR + VARIANCE_ADDITION = 0.5 + return sc.DataArray( - data=sc.Variable(dims=["tof"], values=counts, variances=counts, unit=sc.units.counts), + data=sc.Variable(dims=["tof"], values=counts, variances=counts+VARIANCE_ADDITION, unit=sc.units.counts), coords={"tof": sc.array(dims=["tof"], values=tof_edges, unit=sc.Unit(unit))}, ) From decb23429da2abd6760c856907c7461b245ab6d7 Mon Sep 17 00:00:00 2001 From: Jack Harper Date: Wed, 6 Nov 2024 14:18:57 +0000 Subject: [PATCH 02/16] add yerr in liveplot wrapper --- .../callbacks/fitting/__init__.py | 2 +- src/ibex_bluesky_core/callbacks/plotting.py | 23 +++++++++++++++++-- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/fitting/__init__.py b/src/ibex_bluesky_core/callbacks/fitting/__init__.py index c15c873..44bf4b2 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/__init__.py +++ b/src/ibex_bluesky_core/callbacks/fitting/__init__.py @@ -64,7 +64,7 @@ def __init__( self.method = method super().__init__( - model=method.model, y=y, independent_vars={"x": x}, update_every=update_every, yerr=yerr + model=method.model, y=y, independent_vars={"x": x}, update_every=update_every ) def update_fit(self) -> None: diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index 4013c5d..5d23cf7 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -5,7 +5,7 @@ import matplotlib import matplotlib.pyplot as plt from bluesky.callbacks import LivePlot as _DefaultLivePlot -from bluesky.callbacks.core import make_class_safe +from bluesky.callbacks.core import make_class_safe, get_obj_fields from event_model.documents import Event, RunStart logger = logging.getLogger(__name__) @@ -16,7 +16,12 @@ class LivePlot(_DefaultLivePlot): """Live plot, customized for IBEX.""" def __init__(self, y, x=None, yerr=None, *args, **kwargs): - super().__init__(y=y, x=x, yerr=yerr, *args, **kwargs) + super().__init__(y=y, x=x, *args, **kwargs) + if yerr is not None: + self.yerr, *others = get_obj_fields([yerr]) + else: + self.yerr = None + self.yerr_data = [] def _show_plot(self) -> None: # Play nicely with the "normal" backends too - only force show if we're @@ -24,6 +29,20 @@ def _show_plot(self) -> None: if "genie_python" in matplotlib.get_backend(): plt.show() + def event(self, doc): + new_yerr = None if self.yerr is None else doc["data"][self.yerr] + self.update_yerr(new_yerr) + super().event(doc) + + def update_plot(self): + if self.yerr is not None: + self.ax.errorbar(x=self.x_data, y=self.y_data, yerr=self.yerr_data, fmt="none") + super().update_plot() + + def update_yerr(self, y_err): + # super.update_caches(x, y) + self.yerr_data.append(y_err) + def start(self, doc: RunStart) -> None: """Process an start document (delegate to superclass, then show the plot).""" super().start(doc) From 984a1f723df95f1a52c00cbd8763976278086cc7 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Wed, 6 Nov 2024 14:43:19 +0000 Subject: [PATCH 03/16] updates fitting for uncertainties --- .../callbacks/fitting/__init__.py | 44 +++++++++++++++++-- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/fitting/__init__.py b/src/ibex_bluesky_core/callbacks/fitting/__init__.py index 44bf4b2..2e98162 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/__init__.py +++ b/src/ibex_bluesky_core/callbacks/fitting/__init__.py @@ -2,12 +2,14 @@ import logging from typing import Callable +import warnings import lmfit import numpy as np import numpy.typing as npt from bluesky.callbacks import LiveFit as _DefaultLiveFit from bluesky.callbacks.core import make_class_safe +from event_model.documents.event import Event class FitMethod: @@ -50,7 +52,7 @@ class LiveFit(_DefaultLiveFit): """Live fit, customized for IBEX.""" def __init__( - self, method: FitMethod, y: str, x: str, *, update_every: int = 1, yerr=None + self, method: FitMethod, y: str, x: str, *, update_every: int = 1, yerr: str | None = None ) -> None: """Call Bluesky LiveFit with assumption that there is only one independant variable. @@ -58,15 +60,38 @@ def __init__( method (FitMethod): The FitMethod (Model & Guess) to use when fitting. y (str): The name of the dependant variable. x (str): The name of the independant variable. - update_every (int): How often to update the fit. (seconds) + update_every (int, optional): How often to update the fit. (seconds) + yerr (str or None, optional): Name of field in the Event document that provides standard deviation for each Y value """ self.method = method + self.yerr = yerr + self.weight_data = [] super().__init__( model=method.model, y=y, independent_vars={"x": x}, update_every=update_every ) + def event(self, doc: Event): + + weight = None + if self.yerr is not None: + try: + weight = 1 / doc["data"][self.yerr] + except ZeroDivisionError: + warnings.warn(f"standard deviation for y is 0, therefore applying weight of 0 on fit", + stacklevel=1, + ) + weight = 0.0 + + self.update_weight(weight) + super().event(doc) + + def update_weight(self, weight: float | None = 0.0): + if self.yerr is not None: + self.weight_data.append(weight) + + def update_fit(self) -> None: """Use the provided guess function with the most recent x and y values after every update. @@ -83,4 +108,17 @@ def update_fit(self) -> None: # Calls the guess function on the set of data already collected in the run ) - super().update_fit() + N = len(self.model.param_names) + if len(self.ydata) < N: + warnings.warn( + f"LiveFitPlot cannot update fit until there are at least {N} data points", + stacklevel=1, + ) + else: + kwargs = {} + kwargs.update(self.independent_vars_data) + kwargs.update(self.init_guess) + self.result = self.model.fit( + self.ydata, weights=None if self.yerr is None else self.weight_data, **kwargs + ) + self.__stale = False From 4bc0764e384cf5b32a84390f968863535bbea7a2 Mon Sep 17 00:00:00 2001 From: Jack Harper Date: Wed, 6 Nov 2024 16:14:12 +0000 Subject: [PATCH 04/16] add test for fitting when weights are given --- src/ibex_bluesky_core/callbacks/plotting.py | 10 ++--- .../devices/dae/dae_spectra.py | 3 +- .../fitting/test_fitting_callback.py | 40 +++++++++++++++++++ tests/devices/test_dae.py | 4 +- 4 files changed, 47 insertions(+), 10 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index 5d23cf7..1d13a68 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -29,10 +29,12 @@ def _show_plot(self) -> None: if "genie_python" in matplotlib.get_backend(): plt.show() - def event(self, doc): + def event(self, doc: Event): + """Process an event document (delegate to superclass, then show the plot).""" new_yerr = None if self.yerr is None else doc["data"][self.yerr] self.update_yerr(new_yerr) super().event(doc) + self._show_plot() def update_plot(self): if self.yerr is not None: @@ -47,9 +49,3 @@ def start(self, doc: RunStart) -> None: """Process an start document (delegate to superclass, then show the plot).""" super().start(doc) self._show_plot() - - def event(self, doc: Event) -> None: - """Process an event document (delegate to superclass, then show the plot).""" - - super().event(doc) - self._show_plot() diff --git a/src/ibex_bluesky_core/devices/dae/dae_spectra.py b/src/ibex_bluesky_core/devices/dae/dae_spectra.py index 0fc31e9..07a2535 100644 --- a/src/ibex_bluesky_core/devices/dae/dae_spectra.py +++ b/src/ibex_bluesky_core/devices/dae/dae_spectra.py @@ -9,6 +9,7 @@ from ophyd_async.core import SignalR, StandardReadable from ophyd_async.epics.signal import epics_signal_r +VARIANCE_ADDITION = 0.5 class DaeSpectra(StandardReadable): """Subdevice for a single DAE spectra.""" @@ -108,7 +109,7 @@ async def read_spectrum_dataarray(self) -> sc.DataArray: raise ValueError("Could not determine engineering units of tof edges.") # TODO add reference to ADR - VARIANCE_ADDITION = 0.5 + return sc.DataArray( data=sc.Variable( diff --git a/tests/callbacks/fitting/test_fitting_callback.py b/tests/callbacks/fitting/test_fitting_callback.py index a59da49..c0f5402 100644 --- a/tests/callbacks/fitting/test_fitting_callback.py +++ b/tests/callbacks/fitting/test_fitting_callback.py @@ -1,4 +1,5 @@ from unittest import mock +from unittest.mock import patch, MagicMock import lmfit import numpy as np @@ -92,3 +93,42 @@ def model(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: ) model_mock.assert_called() + +# @patch('ibex_bluesky_core.callbacks.fitting.lmfit') +def test_model_called_with_weights_if_yerr_is_given(): + # model_mock = mock.MagicMock() + # + def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, lmfit.Parameter]: + return {} + # + # def model(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + # model_mock() + # return x + + # fit = mock.MagicMock() + # type(fit).independent_vars = {"x"} + # Does not pass the function to lmfit before FitMethod + model = lmfit.Model(lambda x:x) + model.fit = MagicMock() + method = FitMethod(model=model, guess=guess) + lf = LiveFit(method, y="y", x="x", yerr="yerr") + + + x = 1 + y = 2 + yerr = 3 + + # Fake just enough of an event document. + lf.event( + { + "data": { # type: ignore + "y": y, + "x": x, + "yerr": yerr + } + } + ) + + lf.update_fit() + + model.fit.assert_called_with([y], weights = [1/yerr], x=[1]) diff --git a/tests/devices/test_dae.py b/tests/devices/test_dae.py index 4a4bc7d..f18773a 100644 --- a/tests/devices/test_dae.py +++ b/tests/devices/test_dae.py @@ -25,7 +25,7 @@ DaeSettingsData, TimingSource, ) -from ibex_bluesky_core.devices.dae.dae_spectra import DaeSpectra +from ibex_bluesky_core.devices.dae.dae_spectra import DaeSpectra, VARIANCE_ADDITION from ibex_bluesky_core.devices.dae.dae_tcb_settings import ( CalculationMethod, DaeTCBSettings, @@ -979,7 +979,7 @@ async def test_read_spectrum_dataarray(spectrum: DaeSpectra): data=sc.Variable( dims=["tof"], values=[1000, 2000, 3000], - variances=[1000, 2000, 3000], + variances=[1000 + VARIANCE_ADDITION, 2000 + VARIANCE_ADDITION, 3000+ VARIANCE_ADDITION], unit=sc.units.counts, dtype="float32", ), From 768fd0f4b6839a5f9ca6889fc7406744bf387175 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Wed, 6 Nov 2024 17:20:09 +0000 Subject: [PATCH 05/16] unit tests for full coverage --- .../callbacks/fitting/__init__.py | 11 ++-- .../fitting/test_fitting_callback.py | 59 ++++++++++++++----- tests/callbacks/test_plotting_callback.py | 56 ++++++++++++++++++ 3 files changed, 107 insertions(+), 19 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/fitting/__init__.py b/src/ibex_bluesky_core/callbacks/fitting/__init__.py index 2e98162..7c58e2a 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/__init__.py +++ b/src/ibex_bluesky_core/callbacks/fitting/__init__.py @@ -102,11 +102,6 @@ def update_fit(self) -> None: None """ - self.init_guess = self.method.guess( - np.array(next(iter(self.independent_vars_data.values()))), - np.array(self.ydata), - # Calls the guess function on the set of data already collected in the run - ) N = len(self.model.param_names) if len(self.ydata) < N: @@ -115,6 +110,12 @@ def update_fit(self) -> None: stacklevel=1, ) else: + self.init_guess = self.method.guess( + np.array(next(iter(self.independent_vars_data.values()))), + np.array(self.ydata), + # Calls the guess function on the set of data already collected in the run + ) + kwargs = {} kwargs.update(self.independent_vars_data) kwargs.update(self.init_guess) diff --git a/tests/callbacks/fitting/test_fitting_callback.py b/tests/callbacks/fitting/test_fitting_callback.py index c0f5402..4d297a1 100644 --- a/tests/callbacks/fitting/test_fitting_callback.py +++ b/tests/callbacks/fitting/test_fitting_callback.py @@ -1,11 +1,13 @@ from unittest import mock from unittest.mock import patch, MagicMock +import warnings import lmfit import numpy as np import numpy.typing as npt from ibex_bluesky_core.callbacks.fitting import FitMethod, LiveFit +from ibex_bluesky_core.callbacks.fitting.fitting_utils import Linear def test_guess_called(): @@ -94,31 +96,20 @@ def model(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: model_mock.assert_called() -# @patch('ibex_bluesky_core.callbacks.fitting.lmfit') def test_model_called_with_weights_if_yerr_is_given(): - # model_mock = mock.MagicMock() - # + def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, lmfit.Parameter]: return {} - # - # def model(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: - # model_mock() - # return x - # fit = mock.MagicMock() - # type(fit).independent_vars = {"x"} - # Does not pass the function to lmfit before FitMethod model = lmfit.Model(lambda x:x) model.fit = MagicMock() method = FitMethod(model=model, guess=guess) lf = LiveFit(method, y="y", x="x", yerr="yerr") - x = 1 y = 2 yerr = 3 - # Fake just enough of an event document. lf.event( { "data": { # type: ignore @@ -129,6 +120,46 @@ def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, l } ) - lf.update_fit() - model.fit.assert_called_with([y], weights = [1/yerr], x=[1]) + + +def test_warning_given_if_yerr_is_0(): + + def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, lmfit.Parameter]: + return {} + + model = lmfit.Model(lambda x:x) + model.fit = MagicMock() + method = FitMethod(model=model, guess=guess) + lf = LiveFit(method, y="y", x="x", yerr="yerr") + + x = 1 + y = 2 + yerr = 0 + + with warnings.catch_warnings(record=True) as w: + lf.event( + { + "data": { # type: ignore + "y": y, + "x": x, + "yerr": yerr + } + } + ) + + model.fit.assert_called_with([y], weights=[0.0], x=[1]) + + assert len(w) == 1 + assert "standard deviation for y is 0, therefore applying weight of 0 on fit" in str(w[-1].message) + + +def test_warning_if_no_y_data(): + + with warnings.catch_warnings(record=True) as w: + + lf = LiveFit(Linear.fit(), y="y", x="x", yerr="yerr") + lf.update_fit() + + assert len(w) == 1 + assert "LiveFitPlot cannot update fit until there are at least" in str(w[-1].message) \ No newline at end of file diff --git a/tests/callbacks/test_plotting_callback.py b/tests/callbacks/test_plotting_callback.py index e5b486a..ce0fcbb 100644 --- a/tests/callbacks/test_plotting_callback.py +++ b/tests/callbacks/test_plotting_callback.py @@ -1,6 +1,8 @@ from typing import Any from unittest.mock import MagicMock, patch +from matplotlib import pyplot as plt + from ibex_bluesky_core.callbacks.plotting import LivePlot @@ -42,3 +44,57 @@ def test_show_plot_only_shows_if_backend_is_genie(): mock_get_backend.return_value = "simulated_genie_python_backend" lp._show_plot() mock_plt_show.assert_called_once() + + +def test_errorbars_created_if_yerr_is_given(): + + _, ax = plt.subplots() + ax.errorbar = MagicMock() + + lp = LivePlot(y="y", x="x", yerr="yerr", ax=ax) + + x = 1 + y = 2 + yerr = 3 + + # Fake just enough of an event document. + lp.start( # type: ignore + { + "time": 0, + "plan_name": "scan", + "uid": 0, + "scan_id": 0, + } + ) + + # Fake just enough of an event document. + lp.event( + { + "data": { # type: ignore + "y": y, + "x": x, + "yerr": yerr + } + } + ) + + lp.ax.errorbar.assert_called_with(y=[y], x=[x], yerr=[yerr], fmt="none") # type: ignore + +def test_errorbars_not_created_if_no_yerr(): + + _, ax = plt.subplots() + ax.errorbar = MagicMock() + + lp = LivePlot(y="y", x="x", ax=ax) + + lp.start( # type: ignore + { + "time": 0, + "plan_name": "scan", + "uid": 0, + "scan_id": 0, + } + ) + + lp.update_plot() + assert not lp.ax.errorbar.called \ No newline at end of file From 2254809bfc9b93fe559e4e3976e0f89f411b5687 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Thu, 7 Nov 2024 09:49:53 +0000 Subject: [PATCH 06/16] some ruff fixes --- .../callbacks/fitting/__init__.py | 23 ++++++++------- src/ibex_bluesky_core/callbacks/plotting.py | 2 +- .../fitting/test_fitting_callback.py | 29 +++++++++---------- tests/callbacks/test_plotting_callback.py | 13 ++++----- 4 files changed, 33 insertions(+), 34 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/fitting/__init__.py b/src/ibex_bluesky_core/callbacks/fitting/__init__.py index 7c58e2a..4300d03 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/__init__.py +++ b/src/ibex_bluesky_core/callbacks/fitting/__init__.py @@ -1,8 +1,8 @@ """For IBEX Bluesky scan fitting.""" import logging -from typing import Callable import warnings +from typing import Callable import lmfit import numpy as np @@ -61,7 +61,8 @@ def __init__( y (str): The name of the dependant variable. x (str): The name of the independant variable. update_every (int, optional): How often to update the fit. (seconds) - yerr (str or None, optional): Name of field in the Event document that provides standard deviation for each Y value + yerr (str or None, optional): Name of field in the Event document + that provides standard deviation for each Y value """ self.method = method @@ -72,14 +73,15 @@ def __init__( model=method.model, y=y, independent_vars={"x": x}, update_every=update_every ) - def event(self, doc: Event): - + def event(self, doc: Event) -> None: + """When an event is received, update caches.""" weight = None if self.yerr is not None: try: weight = 1 / doc["data"][self.yerr] except ZeroDivisionError: - warnings.warn(f"standard deviation for y is 0, therefore applying weight of 0 on fit", + warnings.warn( + "standard deviation for y is 0, therefore applying weight of 0 on fit", stacklevel=1, ) weight = 0.0 @@ -87,10 +89,10 @@ def event(self, doc: Event): self.update_weight(weight) super().event(doc) - def update_weight(self, weight: float | None = 0.0): + def update_weight(self, weight: float | None = 0.0) -> None: + """Update uncertainties cache.""" if self.yerr is not None: self.weight_data.append(weight) - def update_fit(self) -> None: """Use the provided guess function with the most recent x and y values after every update. @@ -102,11 +104,10 @@ def update_fit(self) -> None: None """ - - N = len(self.model.param_names) - if len(self.ydata) < N: + n = len(self.model.param_names) + if len(self.ydata) < n: warnings.warn( - f"LiveFitPlot cannot update fit until there are at least {N} data points", + f"LiveFitPlot cannot update fit until there are at least {n} data points", stacklevel=1, ) else: diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index 1d13a68..c998cbc 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -5,7 +5,7 @@ import matplotlib import matplotlib.pyplot as plt from bluesky.callbacks import LivePlot as _DefaultLivePlot -from bluesky.callbacks.core import make_class_safe, get_obj_fields +from bluesky.callbacks.core import get_obj_fields, make_class_safe from event_model.documents import Event, RunStart logger = logging.getLogger(__name__) diff --git a/tests/callbacks/fitting/test_fitting_callback.py b/tests/callbacks/fitting/test_fitting_callback.py index 4d297a1..05d2d31 100644 --- a/tests/callbacks/fitting/test_fitting_callback.py +++ b/tests/callbacks/fitting/test_fitting_callback.py @@ -1,6 +1,6 @@ -from unittest import mock -from unittest.mock import patch, MagicMock import warnings +from unittest import mock +from unittest.mock import MagicMock import lmfit import numpy as np @@ -96,12 +96,12 @@ def model(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: model_mock.assert_called() -def test_model_called_with_weights_if_yerr_is_given(): +def test_model_called_with_weights_if_yerr_is_given(): def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, lmfit.Parameter]: return {} - model = lmfit.Model(lambda x:x) + model = lmfit.Model(lambda x: x) model.fit = MagicMock() method = FitMethod(model=model, guess=guess) lf = LiveFit(method, y="y", x="x", yerr="yerr") @@ -115,20 +115,19 @@ def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, l "data": { # type: ignore "y": y, "x": x, - "yerr": yerr + "yerr": yerr, } } ) - model.fit.assert_called_with([y], weights = [1/yerr], x=[1]) + model.fit.assert_called_with([y], weights=[1 / yerr], x=[1]) def test_warning_given_if_yerr_is_0(): - def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, lmfit.Parameter]: return {} - model = lmfit.Model(lambda x:x) + model = lmfit.Model(lambda x: x) model.fit = MagicMock() method = FitMethod(model=model, guess=guess) lf = LiveFit(method, y="y", x="x", yerr="yerr") @@ -143,23 +142,23 @@ def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, l "data": { # type: ignore "y": y, "x": x, - "yerr": yerr + "yerr": yerr, } } ) model.fit.assert_called_with([y], weights=[0.0], x=[1]) - + assert len(w) == 1 - assert "standard deviation for y is 0, therefore applying weight of 0 on fit" in str(w[-1].message) + assert "standard deviation for y is 0, therefore applying weight of 0 on fit" in str( + w[-1].message + ) def test_warning_if_no_y_data(): - with warnings.catch_warnings(record=True) as w: - lf = LiveFit(Linear.fit(), y="y", x="x", yerr="yerr") lf.update_fit() - + assert len(w) == 1 - assert "LiveFitPlot cannot update fit until there are at least" in str(w[-1].message) \ No newline at end of file + assert "LiveFitPlot cannot update fit until there are at least" in str(w[-1].message) diff --git a/tests/callbacks/test_plotting_callback.py b/tests/callbacks/test_plotting_callback.py index ce0fcbb..916563b 100644 --- a/tests/callbacks/test_plotting_callback.py +++ b/tests/callbacks/test_plotting_callback.py @@ -47,7 +47,6 @@ def test_show_plot_only_shows_if_backend_is_genie(): def test_errorbars_created_if_yerr_is_given(): - _, ax = plt.subplots() ax.errorbar = MagicMock() @@ -58,7 +57,7 @@ def test_errorbars_created_if_yerr_is_given(): yerr = 3 # Fake just enough of an event document. - lp.start( # type: ignore + lp.start( # type: ignore { "time": 0, "plan_name": "scan", @@ -73,21 +72,21 @@ def test_errorbars_created_if_yerr_is_given(): "data": { # type: ignore "y": y, "x": x, - "yerr": yerr + "yerr": yerr, } } ) - lp.ax.errorbar.assert_called_with(y=[y], x=[x], yerr=[yerr], fmt="none") # type: ignore + lp.ax.errorbar.assert_called_with(y=[y], x=[x], yerr=[yerr], fmt="none") # type: ignore -def test_errorbars_not_created_if_no_yerr(): +def test_errorbars_not_created_if_no_yerr(): _, ax = plt.subplots() ax.errorbar = MagicMock() lp = LivePlot(y="y", x="x", ax=ax) - lp.start( # type: ignore + lp.start( # type: ignore { "time": 0, "plan_name": "scan", @@ -97,4 +96,4 @@ def test_errorbars_not_created_if_no_yerr(): ) lp.update_plot() - assert not lp.ax.errorbar.called \ No newline at end of file + assert not lp.ax.errorbar.called From 426aa24e852447e840061d04ed712b377539b897 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Thu, 7 Nov 2024 14:05:35 +0000 Subject: [PATCH 07/16] ADR draft 1 --- .../005-variance-addition.md | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 doc/architectural_decisions/005-variance-addition.md diff --git a/doc/architectural_decisions/005-variance-addition.md b/doc/architectural_decisions/005-variance-addition.md new file mode 100644 index 0000000..74d7b07 --- /dev/null +++ b/doc/architectural_decisions/005-variance-addition.md @@ -0,0 +1,17 @@ +# Why do we add 0.5 to variance? + +## Status + +Current + +## Context + +If we pass counts data to `scipp` as variance this is correct- but if one of the counts is 0 then its variance is 0, leading to a division by 0 error when calculating its weight for fitting `weight = 1 / doc["data"][self.yerr]`. + +## Decision + +Our solution was to add 0.5 (`VARIANCE_ADDITION`) to each count to calculate variance. The actual data should be unchanged, the +0.5 is only for uncertainty calculation. + +## Justification + +The above approach is both "smooth" and converges towards sqrt(N) in the limit with high counts, and should also mean that we never get an uncertainty of zero in the fitting side. \ No newline at end of file From 603d530e8aca23d7776492e46adc1a0bac8732ca Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Thu, 7 Nov 2024 14:15:28 +0000 Subject: [PATCH 08/16] update variance addition adr --- doc/architectural_decisions/005-variance-addition.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/architectural_decisions/005-variance-addition.md b/doc/architectural_decisions/005-variance-addition.md index 74d7b07..c0b5842 100644 --- a/doc/architectural_decisions/005-variance-addition.md +++ b/doc/architectural_decisions/005-variance-addition.md @@ -6,7 +6,7 @@ Current ## Context -If we pass counts data to `scipp` as variance this is correct- but if one of the counts is 0 then its variance is 0, leading to a division by 0 error when calculating its weight for fitting `weight = 1 / doc["data"][self.yerr]`. +If we pass counts data to `scipp` as variance this is correct- but if one of the counts is 0 then its variance is 0, leading to a division by 0 error when calculating its weight for fitting `weight = 1 / doc["data"][self.yerr]`. See `src\ibex_bluesky_core\devices\dae\dae_spectra.py line 118`. ## Decision From 0904ca0c3ecd68b96973bc61a7894a4e8c5a184d Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Thu, 7 Nov 2024 16:31:03 +0000 Subject: [PATCH 09/16] ruff --- src/ibex_bluesky_core/callbacks/plotting.py | 34 ++++++++++++++++----- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index c998cbc..3b7912c 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -1,6 +1,7 @@ """IBEX plotting callbacks.""" import logging +from typing import Any import matplotlib import matplotlib.pyplot as plt @@ -15,7 +16,24 @@ class LivePlot(_DefaultLivePlot): """Live plot, customized for IBEX.""" - def __init__(self, y, x=None, yerr=None, *args, **kwargs): + def __init__( + self, + y: str, + x: str | None = None, + yerr: str | None = None, + *args: Any, # noqa: ANN401 + **kwargs: Any, # noqa: ANN401 + ) -> None: + """Initialise LivePlot. + + Args: + y (str): The name of the dependant variable. + x (str or None, optional): The name of the independant variable. + yerr (str or None, optional): Name of uncertainties signal. + *args: As per mpl_plotting.py + **kwargs: As per mpl_plotting.py + + """ super().__init__(y=y, x=x, *args, **kwargs) if yerr is not None: self.yerr, *others = get_obj_fields([yerr]) @@ -24,26 +42,28 @@ def __init__(self, y, x=None, yerr=None, *args, **kwargs): self.yerr_data = [] def _show_plot(self) -> None: + """Call plt.show().""" # Play nicely with the "normal" backends too - only force show if we're # actually using our custom backend. if "genie_python" in matplotlib.get_backend(): plt.show() - def event(self, doc: Event): + def event(self, doc: Event) -> None: """Process an event document (delegate to superclass, then show the plot).""" new_yerr = None if self.yerr is None else doc["data"][self.yerr] self.update_yerr(new_yerr) super().event(doc) self._show_plot() - def update_plot(self): + def update_plot(self) -> None: + """Create error bars if needed, then update plot.""" if self.yerr is not None: - self.ax.errorbar(x=self.x_data, y=self.y_data, yerr=self.yerr_data, fmt="none") + self.ax.errorbar(x=self.x_data, y=self.y_data, yerr=self.yerr_data, fmt="none") # type: ignore super().update_plot() - def update_yerr(self, y_err): - # super.update_caches(x, y) - self.yerr_data.append(y_err) + def update_yerr(self, yerr: float | None) -> None: + """Update uncertainties data.""" + self.yerr_data.append(yerr) def start(self, doc: RunStart) -> None: """Process an start document (delegate to superclass, then show the plot).""" From 043d1a485a98c63c612d9436e26f18042ad488d3 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Thu, 7 Nov 2024 16:37:56 +0000 Subject: [PATCH 10/16] Ruff --- src/ibex_bluesky_core/devices/dae/dae_spectra.py | 2 +- tests/devices/test_dae.py | 8 ++++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/ibex_bluesky_core/devices/dae/dae_spectra.py b/src/ibex_bluesky_core/devices/dae/dae_spectra.py index 07a2535..535d75c 100644 --- a/src/ibex_bluesky_core/devices/dae/dae_spectra.py +++ b/src/ibex_bluesky_core/devices/dae/dae_spectra.py @@ -11,6 +11,7 @@ VARIANCE_ADDITION = 0.5 + class DaeSpectra(StandardReadable): """Subdevice for a single DAE spectra.""" @@ -110,7 +111,6 @@ async def read_spectrum_dataarray(self) -> sc.DataArray: # TODO add reference to ADR - return sc.DataArray( data=sc.Variable( dims=["tof"], diff --git a/tests/devices/test_dae.py b/tests/devices/test_dae.py index f18773a..5b43daf 100644 --- a/tests/devices/test_dae.py +++ b/tests/devices/test_dae.py @@ -25,7 +25,7 @@ DaeSettingsData, TimingSource, ) -from ibex_bluesky_core.devices.dae.dae_spectra import DaeSpectra, VARIANCE_ADDITION +from ibex_bluesky_core.devices.dae.dae_spectra import VARIANCE_ADDITION, DaeSpectra from ibex_bluesky_core.devices.dae.dae_tcb_settings import ( CalculationMethod, DaeTCBSettings, @@ -979,7 +979,11 @@ async def test_read_spectrum_dataarray(spectrum: DaeSpectra): data=sc.Variable( dims=["tof"], values=[1000, 2000, 3000], - variances=[1000 + VARIANCE_ADDITION, 2000 + VARIANCE_ADDITION, 3000+ VARIANCE_ADDITION], + variances=[ + 1000 + VARIANCE_ADDITION, + 2000 + VARIANCE_ADDITION, + 3000 + VARIANCE_ADDITION, + ], unit=sc.units.counts, dtype="float32", ), From 086c876884697c8e1c676654730ac5a3089f6342 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Fri, 8 Nov 2024 11:30:01 +0000 Subject: [PATCH 11/16] restore dae_scan --- manual_system_tests/dae_scan.py | 60 ++++++++++++++++----------------- 1 file changed, 29 insertions(+), 31 deletions(-) diff --git a/manual_system_tests/dae_scan.py b/manual_system_tests/dae_scan.py index 6b25729..f4d12cc 100644 --- a/manual_system_tests/dae_scan.py +++ b/manual_system_tests/dae_scan.py @@ -8,17 +8,15 @@ import bluesky.plans as bp import matplotlib import matplotlib.pyplot as plt -from bluesky.callbacks import LiveTable, LiveFitPlot +from bluesky.callbacks import LiveTable from bluesky.preprocessors import subs_decorator from bluesky.utils import Msg from ophyd_async.plan_stubs import ensure_connected from ibex_bluesky_core.callbacks.file_logger import HumanReadableFileCallback -from ibex_bluesky_core.callbacks.fitting import LiveFit -from ibex_bluesky_core.callbacks.fitting.fitting_utils import Linear, Gaussian from ibex_bluesky_core.callbacks.plotting import LivePlot from ibex_bluesky_core.devices import get_pv_prefix -from ibex_bluesky_core.devices.block import block_rw_rbv, BlockWriteConfig +from ibex_bluesky_core.devices.block import block_rw_rbv from ibex_bluesky_core.devices.simpledae import SimpleDae from ibex_bluesky_core.devices.simpledae.controllers import ( RunPerPointController, @@ -49,12 +47,10 @@ def dae_scan_plan() -> Generator[Msg, None, None]: - The DAE waited for at least 500 good frames at each point """ prefix = get_pv_prefix() - block = block_rw_rbv(float, "mot", write_config=BlockWriteConfig(settle_time_s=0.5)) - blocka = block_rw_rbv(float, "alice") - blockb = block_rw_rbv(float, "bob", write_config=BlockWriteConfig(settle_time_s=0.5)) + block = block_rw_rbv(float, "mot") controller = RunPerPointController(save_run=True) - waiter = GoodFramesWaiter(1) + waiter = GoodFramesWaiter(500) reducer = GoodFramesNormalizer( prefix=prefix, detector_spectra=[i for i in range(1, 100)], @@ -71,36 +67,38 @@ def dae_scan_plan() -> Generator[Msg, None, None]: controller.run_number.set_name("run number") reducer.intensity.set_name("normalized counts") - yield from ensure_connected(block, blocka, blockb, force_reconnect=True) - print(reducer.intensity.name) - - lf = LiveFit(Gaussian.fit(), y=blockb.name, x=block.name, yerr=blocka.name) - fig, ax = plt.subplots() + yield from ensure_connected(block, dae, force_reconnect=True) @subs_decorator( [ - # HumanReadableFileCallback( - # Path("C:\\") / "instrument" / "var" / "logs" / "bluesky" / "output_files", - # [ - # block.name, - # controller.run_number.name, - # reducer.intensity.name, - # reducer.det_counts.name, - # dae.good_frames.name, - # ], - # ), - LiveFitPlot(lf, ax=ax, color="r", num_points=1000), - LivePlot( - y=blockb.name, x=block.name, yerr=blocka.name, marker="x", linestyle="none", ax=ax + HumanReadableFileCallback( + Path("C:\\") / "instrument" / "var" / "logs" / "bluesky" / "output_files", + [ + block.name, + controller.run_number.name, + reducer.intensity.name, + reducer.det_counts.name, + dae.good_frames.name, + ], + ), + LivePlot(y=reducer.intensity.name, x=block.name, marker="x", linestyle="none"), + LiveTable( + [ + block.name, + controller.run_number.name, + reducer.intensity.name, + reducer.intensity_stddev.name, + reducer.det_counts.name, + reducer.det_counts_stddev.name, + dae.good_frames.name, + ] ), - LiveTable([block.name, blocka.name, blockb.name]), ] ) def _inner() -> Generator[Msg, None, None]: - num_points = 10 - # yield from bps.mv(dae.number_of_periods, num_points) - yield from bps.read(blocka) - yield from bp.scan([blockb, blocka], block, 0, 10, num=num_points) + num_points = 3 + yield from bps.mv(dae.number_of_periods, num_points) + yield from bp.scan([dae], block, 0, 10, num=num_points) yield from _inner() From 15b310226988dd6285184a3f862f4431f799549a Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Fri, 8 Nov 2024 15:44:05 +0000 Subject: [PATCH 12/16] Ruff --- src/ibex_bluesky_core/callbacks/plotting.py | 2 +- tests/callbacks/test_plotting_callback.py | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index 3b7912c..655ec11 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -34,7 +34,7 @@ def __init__( **kwargs: As per mpl_plotting.py """ - super().__init__(y=y, x=x, *args, **kwargs) + super().__init__(y=y, x=x, *args, **kwargs) # noqa: B026 if yerr is not None: self.yerr, *others = get_obj_fields([yerr]) else: diff --git a/tests/callbacks/test_plotting_callback.py b/tests/callbacks/test_plotting_callback.py index 916563b..b53e5df 100644 --- a/tests/callbacks/test_plotting_callback.py +++ b/tests/callbacks/test_plotting_callback.py @@ -57,11 +57,11 @@ def test_errorbars_created_if_yerr_is_given(): yerr = 3 # Fake just enough of an event document. - lp.start( # type: ignore + lp.start( { "time": 0, "plan_name": "scan", - "uid": 0, + "uid": 0, # type: ignore "scan_id": 0, } ) @@ -69,12 +69,12 @@ def test_errorbars_created_if_yerr_is_given(): # Fake just enough of an event document. lp.event( { - "data": { # type: ignore + "data": { "y": y, "x": x, "yerr": yerr, } - } + } # type: ignore ) lp.ax.errorbar.assert_called_with(y=[y], x=[x], yerr=[yerr], fmt="none") # type: ignore @@ -86,14 +86,14 @@ def test_errorbars_not_created_if_no_yerr(): lp = LivePlot(y="y", x="x", ax=ax) - lp.start( # type: ignore + lp.start( { "time": 0, "plan_name": "scan", - "uid": 0, + "uid": 0, # type: ignore "scan_id": 0, } ) lp.update_plot() - assert not lp.ax.errorbar.called + assert not lp.ax.errorbar.called # type: ignore From 02ffb3e15bb146385a7cb93532a1b92923e38d84 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Fri, 8 Nov 2024 15:49:00 +0000 Subject: [PATCH 13/16] ruff & docs update --- .../005-variance-addition.md | 49 +++++++++++++++++-- src/ibex_bluesky_core/callbacks/plotting.py | 2 +- 2 files changed, 46 insertions(+), 5 deletions(-) diff --git a/doc/architectural_decisions/005-variance-addition.md b/doc/architectural_decisions/005-variance-addition.md index c0b5842..d9ee84f 100644 --- a/doc/architectural_decisions/005-variance-addition.md +++ b/doc/architectural_decisions/005-variance-addition.md @@ -1,4 +1,4 @@ -# Why do we add 0.5 to variance? +# Variance addition to counts data ## Status @@ -6,12 +6,53 @@ Current ## Context -If we pass counts data to `scipp` as variance this is correct- but if one of the counts is 0 then its variance is 0, leading to a division by 0 error when calculating its weight for fitting `weight = 1 / doc["data"][self.yerr]`. See `src\ibex_bluesky_core\devices\dae\dae_spectra.py line 118`. +For counts data, the uncertainty on counts is typically defined by poisson counting statistics, i.e. the standard deviation on `N` counts is `sqrt(N)`. + +This can be problematic in cases where zero counts have been collected, as the standard deviation will then be zero, which will subsequently lead to "infinite" point weightings in downstream fitting routines for example. + +A number of possible approaches were considered: + +| Option | Description | +| --- | --- | +| A | Reject data with zero counts, i.e. explicitly throw an exception if any data with zero counts is seen as part of a scan. | +| B | Use a standard deviation of `NaN` for points with zero counts. | +| C | Define the standard deviation of `N` counts as `1` if counts are zero, otherwise `sqrt(N)`. This is one of the approaches available in mantid for example. | +| D | Define the standard deviation of `N` counts as `sqrt(N+0.5)` unconditionally - on the basis that "half a count" is smaller than the smallest possible actual measurement which can be taken. | +| E | No special handling, calculate std. dev. as `sqrt(N)`. | + +For clarity, the following table shows the value and associated uncertainty for each option: + +| Counts | Std. Dev. (A) | Std. Dev. (B) | Std. Dev. (C) | Std. Dev. (D) | Std. Dev. (E) | +| ------- | ------ | ------- | ------- | ------- | --- | +| 0 | raise exception | NaN | 1 | 0.707 | 0 | +| 1 | 1 | 1 | 1 | 1.224745 | 1 | +| 2 | 1.414214 | 1.414214 | 1.414214 | 1.581139 | 1.414214 | +| 3 | 1.732051 | 1.732051 | 1.732051 | 1.870829 | 1.732051 | +| 4 | 2 | 2 | 2 | 2.12132 | 2 | +| 5 | 2.236068 | 2.236068 | 2.236068 | 2.345208 | 2.236068 | +| 10 | 3.162278 | 3.162278 | 3.162278 | 3.24037 | 3.162278 | +| 50 | 7.071068 | 7.071068 | 7.071068 | 7.106335 | 7.071068 | +| 100 | 10 | 10 | 10 | 10.02497 | 10 | +| 500 | 22.36068 | 22.36068 | 22.36068 | 22.37186 | 22.36068 | +| 1000 | 31.62278 | 31.62278 | 31.62278 | 31.63068 | 31.62278 | +| 5000 | 70.71068 | 70.71068 | 70.71068 | 70.71421 | 70.71068 | +| 10000 | 100 | 100 | 100 | 100.0025 | 100 | + +## Present + +These approaches were discussed in a regular project update meeting including +- TW & FA (Experiment controls) +- CK (Reflectometry) +- JL (Muons) +- RD (SANS) ## Decision -Our solution was to add 0.5 (`VARIANCE_ADDITION`) to each count to calculate variance. The actual data should be unchanged, the +0.5 is only for uncertainty calculation. +The consensus was to go with Option D. ## Justification -The above approach is both "smooth" and converges towards sqrt(N) in the limit with high counts, and should also mean that we never get an uncertainty of zero in the fitting side. \ No newline at end of file +- Option A will cause real-life scans to crash in low counts regions. +- Option B involves `NaN`s, which have many surprising floating-point characteristics and are highly likely to be a source of future bugs. +- Option D was preferred to option C by scientists present. +- Option E causes surprising results and/or crashes downstream, for example fitting may consider points with zero uncertainty to have "infinite" weight, therefore effectively disregarding all other data. diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index 655ec11..58022a0 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -34,7 +34,7 @@ def __init__( **kwargs: As per mpl_plotting.py """ - super().__init__(y=y, x=x, *args, **kwargs) # noqa: B026 + super().__init__(y=y, x=x, *args, **kwargs) # noqa: B026 if yerr is not None: self.yerr, *others = get_obj_fields([yerr]) else: From a0cafb2b38435738d8dfb0f4fb45d22cf6d091f4 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Tue, 12 Nov 2024 10:49:57 +0000 Subject: [PATCH 14/16] updates docs & adds error bars in manual test --- doc/callbacks/plotting.md | 3 ++- doc/fitting/fitting.md | 17 ++++++++++------- manual_system_tests/dae_scan.py | 23 ++++++++++++++++++++--- 3 files changed, 32 insertions(+), 11 deletions(-) diff --git a/doc/callbacks/plotting.md b/doc/callbacks/plotting.md index 0bca52d..3b53392 100644 --- a/doc/callbacks/plotting.md +++ b/doc/callbacks/plotting.md @@ -38,7 +38,8 @@ ax = plt.gca() # Set the y-scale to logarithmic ax.set_yscale("log") # Use the above axes in a LivePlot callback -plot_callback = LivePlot(y="y_variable", x="x_variable", ax=ax) +plot_callback = LivePlot(y="y_variable", x="x_variable", ax=ax, yerr="yerr_variable") +# yerr is the uncertanties of each y value, producing error bars ``` The `plot_callback` object can then be subscribed to the run engine, using either: diff --git a/doc/fitting/fitting.md b/doc/fitting/fitting.md index f9acbb7..c870f8a 100644 --- a/doc/fitting/fitting.md +++ b/doc/fitting/fitting.md @@ -24,14 +24,17 @@ plt.figure() ax = plt.gca() # ax is shared by fit_callback and plot_callback -plot_callback = LivePlot(y="y_variable", x="x_variable", ax=ax) -fit_callback = LiveFit(Gaussian.fit(), y="y_variable", x="x_variable", update_every=0.5) +plot_callback = LivePlot(y="y_signal", x="x_signal", ax=ax, yerr="yerr_signal_") +fit_callback = LiveFit(Gaussian.fit(), y="y_signal", x="x_signal", yerr="yerr_signal_", update_every=0.5) +# Using the yerr parameter allows you to use error bars. # update_every = in seconds, how often to recompute the fit. If `None`, do not compute until the end. Default is 1. fit_plot_callback = LiveFitPlot(fit_callback, ax=ax, color="r") ``` **Note:** that the `LiveFit` callback doesn't directly do the plotting, it will return function parameters of the model its trying to fit to; a `LiveFit` object must be passed to `LiveFitPlot` which can then be subscribed to the `RunEngine`. See the [Bluesky Documentation](https://blueskyproject.io/bluesky/main/callbacks.html#livefitplot) for information on the various arguments that can be passed to the `LiveFitPlot` class. +Using the `yerr` argument allows you to pass uncertainties via a signal to LiveFit, so that the "weight" of each point influences the fit produced. + The `plot_callback` and `fit_plot_callback` objects can then be subscribed to the `RunEngine`, using the same methods as described in [`LivePlot`](../callbacks/plotting.md). See the following example using `@subs_decorator`: ```py @@ -79,7 +82,7 @@ from bluesky.callbacks import LiveFitPlot from ibex_bluesky_core.callbacks.fitting.fitting_utils import [FIT] # Pass [FIT].fit() to the first parameter of LiveFit -lf = LiveFit([FIT].fit(), y="y_variable", x="x_variable", update_every=0.5) +lf = LiveFit([FIT].fit(), y="y_signal", x="x_signal", update_every=0.5) # Then subscribe to LiveFitPlot(lf, ...) ``` @@ -89,7 +92,7 @@ The `[FIT].fit()` function will pass the `FitMethod` object straight to the `Liv **Note:** that for the fits in the above table that require parameters, you will need to pass value(s) to their `.fit` method. For example Polynomial fitting: ```py -lf = LiveFit(Polynomial.fit(3), y="y_variable", x="x_variable", update_every=0.5) +lf = LiveFit(Polynomial.fit(3), y="y_signal", x="x_signal", update_every=0.5) # For a polynomial of degree 3 ``` @@ -138,7 +141,7 @@ def guess(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> dict[str, l fit_method = FitMethod(model, guess) #Pass the model and guess function to FitMethod -lf = LiveFit(fit_method, y="y_variable", x="x_variable", update_every=0.5) +lf = LiveFit(fit_method, y="y_signal", x="x_signal", update_every=0.5) # Then subscribe to LiveFitPlot(lf, ...) ``` @@ -163,7 +166,7 @@ def different_model(x: float, c1: float, c0: float) -> float: fit_method = FitMethod(different_model, Linear.guess()) # Uses the user defined model and the standard Guessing. function for linear models -lf = LiveFit(fit_method, y="y_variable", x="x_variable", update_every=0.5) +lf = LiveFit(fit_method, y="y_signal", x="x_signal", update_every=0.5) # Then subscribe to LiveFitPlot(lf, ...) ``` @@ -188,7 +191,7 @@ def different_guess(x: float, c1: float, c0: float) -> float: fit_method = FitMethod(Linear.model(), different_guess) # Uses the standard linear model and the user defined Guessing. function -lf = LiveFit(fit_method, y="y_variable", x="x_variable", update_every=0.5) +lf = LiveFit(fit_method, y="y_signal", x="x_signal", update_every=0.5) # Then subscribe to LiveFitPlot(lf, ...) ``` diff --git a/manual_system_tests/dae_scan.py b/manual_system_tests/dae_scan.py index f4d12cc..9a342e9 100644 --- a/manual_system_tests/dae_scan.py +++ b/manual_system_tests/dae_scan.py @@ -8,12 +8,14 @@ import bluesky.plans as bp import matplotlib import matplotlib.pyplot as plt -from bluesky.callbacks import LiveTable +from bluesky.callbacks import LiveFitPlot, LiveTable from bluesky.preprocessors import subs_decorator from bluesky.utils import Msg from ophyd_async.plan_stubs import ensure_connected from ibex_bluesky_core.callbacks.file_logger import HumanReadableFileCallback +from ibex_bluesky_core.callbacks.fitting import LiveFit +from ibex_bluesky_core.callbacks.fitting.fitting_utils import Linear from ibex_bluesky_core.callbacks.plotting import LivePlot from ibex_bluesky_core.devices import get_pv_prefix from ibex_bluesky_core.devices.block import block_rw_rbv @@ -27,6 +29,8 @@ from ibex_bluesky_core.devices.simpledae.waiters import GoodFramesWaiter from ibex_bluesky_core.run_engine import get_run_engine +NUM_POINTS = 3 + def dae_scan_plan() -> Generator[Msg, None, None]: """Manual system test which moves a block and reads the DAE. @@ -67,6 +71,11 @@ def dae_scan_plan() -> Generator[Msg, None, None]: controller.run_number.set_name("run number") reducer.intensity.set_name("normalized counts") + _, ax = plt.subplots() + lf = LiveFit( + Linear.fit(), y=reducer.intensity.name, x=block.name, yerr=reducer.intensity_stddev.name + ) + yield from ensure_connected(block, dae, force_reconnect=True) @subs_decorator( @@ -81,7 +90,15 @@ def dae_scan_plan() -> Generator[Msg, None, None]: dae.good_frames.name, ], ), - LivePlot(y=reducer.intensity.name, x=block.name, marker="x", linestyle="none"), + LiveFitPlot(livefit=lf, ax=ax), + LivePlot( + y=reducer.intensity.name, + x=block.name, + marker="x", + linestyle="none", + ax=ax, + yerr=reducer.intensity_stddev.name, + ), LiveTable( [ block.name, @@ -96,7 +113,7 @@ def dae_scan_plan() -> Generator[Msg, None, None]: ] ) def _inner() -> Generator[Msg, None, None]: - num_points = 3 + num_points = NUM_POINTS yield from bps.mv(dae.number_of_periods, num_points) yield from bp.scan([dae], block, 0, 10, num=num_points) From 8afeeba4c41635c6decaa2f2189cb762925d808f Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Tue, 12 Nov 2024 10:59:17 +0000 Subject: [PATCH 15/16] pyright --- manual_system_tests/dae_scan.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/manual_system_tests/dae_scan.py b/manual_system_tests/dae_scan.py index 9a342e9..8b72d14 100644 --- a/manual_system_tests/dae_scan.py +++ b/manual_system_tests/dae_scan.py @@ -29,7 +29,7 @@ from ibex_bluesky_core.devices.simpledae.waiters import GoodFramesWaiter from ibex_bluesky_core.run_engine import get_run_engine -NUM_POINTS = 3 +NUM_POINTS: int = 3 def dae_scan_plan() -> Generator[Msg, None, None]: @@ -114,7 +114,8 @@ def dae_scan_plan() -> Generator[Msg, None, None]: ) def _inner() -> Generator[Msg, None, None]: num_points = NUM_POINTS - yield from bps.mv(dae.number_of_periods, num_points) + yield from bps.mv(dae.number_of_periods, num_points) # type: ignore + # Pyright does not understand as bluesky isn't typed yet yield from bp.scan([dae], block, 0, 10, num=num_points) yield from _inner() From bba3a6b8343badd9412c62b967d3f0abe9cd4233 Mon Sep 17 00:00:00 2001 From: Jack Doughty Date: Wed, 20 Nov 2024 14:26:11 +0000 Subject: [PATCH 16/16] Requested changes --- doc/callbacks/plotting.md | 2 ++ doc/fitting/fitting.md | 6 +++--- manual_system_tests/dae_scan.py | 5 ++--- src/ibex_bluesky_core/callbacks/fitting/__init__.py | 5 +++-- src/ibex_bluesky_core/callbacks/plotting.py | 9 ++++++--- src/ibex_bluesky_core/devices/dae/dae_spectra.py | 3 ++- tests/callbacks/test_plotting_callback.py | 10 ++++------ 7 files changed, 22 insertions(+), 18 deletions(-) diff --git a/doc/callbacks/plotting.md b/doc/callbacks/plotting.md index 3b53392..5bfffba 100644 --- a/doc/callbacks/plotting.md +++ b/doc/callbacks/plotting.md @@ -42,6 +42,8 @@ plot_callback = LivePlot(y="y_variable", x="x_variable", ax=ax, yerr="yerr_varia # yerr is the uncertanties of each y value, producing error bars ``` +By providing a signal name to the `yerr` argument you can pass uncertainties to LivePlot, by not providing anything for this argument means that no errorbars will be drawn. Errorbars are drawn after each point collected, displaying their standard deviation- uncertainty data is collected from Bluesky event documents and errorbars are updated after every new point added. + The `plot_callback` object can then be subscribed to the run engine, using either: - An explicit callback when calling the run engine: `RE(some_plan(), plot_callback)` - Be subscribed in a plan using `@subs_decorator` from bluesky **(recommended)** diff --git a/doc/fitting/fitting.md b/doc/fitting/fitting.md index c870f8a..d192b64 100644 --- a/doc/fitting/fitting.md +++ b/doc/fitting/fitting.md @@ -24,8 +24,8 @@ plt.figure() ax = plt.gca() # ax is shared by fit_callback and plot_callback -plot_callback = LivePlot(y="y_signal", x="x_signal", ax=ax, yerr="yerr_signal_") -fit_callback = LiveFit(Gaussian.fit(), y="y_signal", x="x_signal", yerr="yerr_signal_", update_every=0.5) +plot_callback = LivePlot(y="y_signal", x="x_signal", ax=ax, yerr="yerr_signal") +fit_callback = LiveFit(Gaussian.fit(), y="y_signal", x="x_signal", yerr="yerr_signal", update_every=0.5) # Using the yerr parameter allows you to use error bars. # update_every = in seconds, how often to recompute the fit. If `None`, do not compute until the end. Default is 1. fit_plot_callback = LiveFitPlot(fit_callback, ax=ax, color="r") @@ -33,7 +33,7 @@ fit_plot_callback = LiveFitPlot(fit_callback, ax=ax, color="r") **Note:** that the `LiveFit` callback doesn't directly do the plotting, it will return function parameters of the model its trying to fit to; a `LiveFit` object must be passed to `LiveFitPlot` which can then be subscribed to the `RunEngine`. See the [Bluesky Documentation](https://blueskyproject.io/bluesky/main/callbacks.html#livefitplot) for information on the various arguments that can be passed to the `LiveFitPlot` class. -Using the `yerr` argument allows you to pass uncertainties via a signal to LiveFit, so that the "weight" of each point influences the fit produced. +Using the `yerr` argument allows you to pass uncertainties via a signal to LiveFit, so that the "weight" of each point influences the fit produced. By not providing a signal name you choose not to use uncertainties/weighting in the fitting calculation. Each weight is computed as `1/(standard deviation at point)` and is taken into account to determine how much a point affects the overall fit of the data. Same as the rest of `LiveFit`, the fit will be updated after every new point collected now taking into account the weights of each point. Uncertainty data is collected from Bluesky event documents after each new point. The `plot_callback` and `fit_plot_callback` objects can then be subscribed to the `RunEngine`, using the same methods as described in [`LivePlot`](../callbacks/plotting.md). See the following example using `@subs_decorator`: diff --git a/manual_system_tests/dae_scan.py b/manual_system_tests/dae_scan.py index 8b72d14..eb9e32a 100644 --- a/manual_system_tests/dae_scan.py +++ b/manual_system_tests/dae_scan.py @@ -113,10 +113,9 @@ def dae_scan_plan() -> Generator[Msg, None, None]: ] ) def _inner() -> Generator[Msg, None, None]: - num_points = NUM_POINTS - yield from bps.mv(dae.number_of_periods, num_points) # type: ignore + yield from bps.mv(dae.number_of_periods, NUM_POINTS) # type: ignore # Pyright does not understand as bluesky isn't typed yet - yield from bp.scan([dae], block, 0, 10, num=num_points) + yield from bp.scan([dae], block, 0, 10, num=NUM_POINTS) yield from _inner() diff --git a/src/ibex_bluesky_core/callbacks/fitting/__init__.py b/src/ibex_bluesky_core/callbacks/fitting/__init__.py index a0ada1f..3d18ea3 100644 --- a/src/ibex_bluesky_core/callbacks/fitting/__init__.py +++ b/src/ibex_bluesky_core/callbacks/fitting/__init__.py @@ -61,7 +61,8 @@ def __init__( x (str): The name of the independant variable. update_every (int, optional): How often to update the fit. (seconds) yerr (str or None, optional): Name of field in the Event document - that provides standard deviation for each Y value + that provides standard deviation for each Y value. None meaning + do not use uncertainties in fit. """ self.method = method @@ -116,7 +117,7 @@ def update_fit(self) -> None: np.array(self.ydata), # Calls the guess function on the set of data already collected in the run ) - + logger.info("new guess for %s: %s", self.method, self.init_guess) kwargs = {} diff --git a/src/ibex_bluesky_core/callbacks/plotting.py b/src/ibex_bluesky_core/callbacks/plotting.py index db1cf1f..5547346 100644 --- a/src/ibex_bluesky_core/callbacks/plotting.py +++ b/src/ibex_bluesky_core/callbacks/plotting.py @@ -30,6 +30,7 @@ def __init__( y (str): The name of the dependant variable. x (str or None, optional): The name of the independant variable. yerr (str or None, optional): Name of uncertainties signal. + Providing None means do not plot uncertainties. *args: As per mpl_plotting.py **kwargs: As per mpl_plotting.py @@ -42,9 +43,11 @@ def __init__( self.yerr_data = [] def _show_plot(self) -> None: - """Call plt.show().""" - # Play nicely with the "normal" backends too - only force show if we're - # actually using our custom backend. + """Call plt.show(). + + Play nicely with the "normal" backends too + - only force show if we're actually using our custom backend. + """ if "genie_python" in matplotlib.get_backend(): logger.debug("Explicitly show()ing plot for IBEX") plt.show() diff --git a/src/ibex_bluesky_core/devices/dae/dae_spectra.py b/src/ibex_bluesky_core/devices/dae/dae_spectra.py index 8e98cfa..2c6ed14 100644 --- a/src/ibex_bluesky_core/devices/dae/dae_spectra.py +++ b/src/ibex_bluesky_core/devices/dae/dae_spectra.py @@ -116,7 +116,8 @@ async def read_spectrum_dataarray(self) -> sc.DataArray: if unit is None: raise ValueError("Could not determine engineering units of tof edges.") - # TODO add reference to ADR + # See doc\architectural_decisions\005-variance-addition.md + # for justfication of the VARIANCE_ADDITION to variances return sc.DataArray( data=sc.Variable( diff --git a/tests/callbacks/test_plotting_callback.py b/tests/callbacks/test_plotting_callback.py index b53e5df..c15658c 100644 --- a/tests/callbacks/test_plotting_callback.py +++ b/tests/callbacks/test_plotting_callback.py @@ -60,8 +60,7 @@ def test_errorbars_created_if_yerr_is_given(): lp.start( { "time": 0, - "plan_name": "scan", - "uid": 0, # type: ignore + "uid": "0", "scan_id": 0, } ) @@ -77,7 +76,7 @@ def test_errorbars_created_if_yerr_is_given(): } # type: ignore ) - lp.ax.errorbar.assert_called_with(y=[y], x=[x], yerr=[yerr], fmt="none") # type: ignore + ax.errorbar.assert_called_with(y=[y], x=[x], yerr=[yerr], fmt="none") def test_errorbars_not_created_if_no_yerr(): @@ -89,11 +88,10 @@ def test_errorbars_not_created_if_no_yerr(): lp.start( { "time": 0, - "plan_name": "scan", - "uid": 0, # type: ignore + "uid": "0", "scan_id": 0, } ) lp.update_plot() - assert not lp.ax.errorbar.called # type: ignore + assert not ax.errorbar.called