Merge pull request #83 from ISISComputingGroup/uncertainties

Set uncertainties more correctly

Author: rerpha

src/ibex_bluesky_core/callbacks/reflectometry/_det_map.py

@@ -33,6 +33,7 @@
 from matplotlib.axes import Axes
 
 from ibex_bluesky_core.callbacks.plotting import show_plot
+from ibex_bluesky_core.devices.simpledae.reducers import VARIANCE_ADDITION
 
 logger = logging.getLogger(__name__)
 
@@ -60,23 +61,23 @@ def event(self, doc: Event, **kwargs: dict[str, Any]) -> Event:
         det_data = doc["data"][self._det_name]
         mon_data = doc["data"][self._mon_name]
 
-        det = sc.Variable(
-            dims=["spectrum"], values=det_data, variances=det_data + 0.5, dtype="float64"
-        )
-        mon = sc.Variable(
-            dims=["spectrum"], values=mon_data, variances=mon_data + 0.5, dtype="float64"
-        )
+        det = sc.Variable(dims=["spectrum"], values=det_data, variances=det_data, dtype="float64")
+        mon = sc.Variable(dims=["spectrum"], values=mon_data, variances=mon_data, dtype="float64")
 
         det_sum = det.sum()
         mon_sum = mon.sum()
 
+        # See doc\architectural_decisions\005-variance-addition.md
+        # for justification of this addition to variances.
+        det_sum.variance += VARIANCE_ADDITION
+
         if mon_sum.value == 0.0:
             raise ValueError(
                 "No monitor counts. Check beamline setup & beam status. "
                 "I/I_0 normalization not possible."
             )
-        else:
-            normalized = det_sum / mon_sum
+
+        normalized = det_sum / mon_sum
 
         doc["data"][self._out_name] = normalized.value
         doc["data"][self._out_name + "_err"] = math.sqrt(normalized.variance)

src/ibex_bluesky_core/devices/dae/dae_spectra.py

@@ -10,7 +10,6 @@
 from ophyd_async.core import Array1D, SignalR, StandardReadable
 from ophyd_async.epics.core import epics_signal_r
 
-VARIANCE_ADDITION = 0.5
 logger = logging.getLogger(__name__)
 
 
@@ -116,14 +115,11 @@ async def read_spectrum_dataarray(self) -> sc.DataArray:
         if unit is None:
             raise ValueError("Could not determine engineering units of tof edges.")
 
-        # See doc\architectural_decisions\005-variance-addition.md
-        # for justfication of the VARIANCE_ADDITION to variances
-
         return sc.DataArray(
             data=sc.Variable(
                 dims=["tof"],
                 values=counts,
-                variances=counts + VARIANCE_ADDITION,
+                variances=counts,
                 unit=sc.units.counts,
                 dtype="float64",
             ),

src/ibex_bluesky_core/devices/simpledae/reducers.py

@@ -30,6 +30,7 @@
 
 
 INTENSITY_PRECISION = 6
+VARIANCE_ADDITION = 0.5
 
 
 async def sum_spectra(spectra: Collection[DaeSpectra]) -> sc.Variable | sc.DataArray:
@@ -176,20 +177,21 @@ async def reduce_data(self, dae: "SimpleDae") -> None:
             self.sum_detector(self.detectors.values()), self.denominator(dae).get_value()
         )
 
-        self._det_counts_setter(float(summed_counts.value))
+        if denominator == 0.0:
+            raise ValueError("Cannot normalize; denominator is zero. Check beamline configuration.")
+
+        # See doc\architectural_decisions\005-variance-addition.md
+        # for justification of this addition to variances.
+        summed_counts.variance += VARIANCE_ADDITION
 
-        if denominator == 0.0:  # To avoid zero division
-            self._intensity_setter(0.0)
-            intensity_var = 0.0
-        else:
-            intensity = summed_counts / denominator
-            self._intensity_setter(intensity.value)
-            intensity_var = intensity.variance if intensity.variance is not None else 0.0
+        intensity = summed_counts / denominator
 
-        detector_counts_var = 0.0 if summed_counts.variance is None else summed_counts.variance
+        self._det_counts_setter(float(summed_counts.value))
+        self._det_counts_stddev_setter(math.sqrt(summed_counts.variance))
+
+        self._intensity_setter(intensity.value)
+        self._intensity_stddev_setter(math.sqrt(intensity.variance))
 
-        self._det_counts_stddev_setter(math.sqrt(detector_counts_var))
-        self._intensity_stddev_setter(math.sqrt(intensity_var))
         logger.info("reduction complete")
 
     def additional_readable_signals(self, dae: "SimpleDae") -> list[Device]:
@@ -284,25 +286,25 @@ async def reduce_data(self, dae: "SimpleDae") -> None:
             self.sum_monitor(self.monitors.values()),
         )
 
-        if monitor_counts.value == 0.0:  # To avoid zero division
-            self._intensity_setter(0.0)
-            intensity_var = 0.0
+        if monitor_counts.value == 0.0:
+            raise ValueError(
+                "Cannot normalize; got zero monitor counts. Check beamline configuration."
+            )
 
-        else:
-            intensity = detector_counts / monitor_counts
-            self._intensity_setter(float(intensity.value))
-            intensity_var = intensity.variance if intensity.variance is not None else 0.0
+        # See doc\architectural_decisions\005-variance-addition.md
+        # for justification of this addition to variances.
+        detector_counts.variance += VARIANCE_ADDITION
 
-        self._intensity_stddev_setter(math.sqrt(intensity_var))
+        intensity = detector_counts / monitor_counts
 
+        self._intensity_setter(float(intensity.value))
         self._det_counts_setter(float(detector_counts.value))
         self._mon_counts_setter(float(monitor_counts.value))
 
-        detector_counts_var = 0.0 if detector_counts.variance is None else detector_counts.variance
-        monitor_counts_var = 0.0 if monitor_counts.variance is None else monitor_counts.variance
+        self._intensity_stddev_setter(math.sqrt(intensity.variance))
+        self._det_counts_stddev_setter(math.sqrt(detector_counts.variance))
+        self._mon_counts_stddev_setter(math.sqrt(monitor_counts.variance))
 
-        self._det_counts_stddev_setter(math.sqrt(detector_counts_var))
-        self._mon_counts_stddev_setter(math.sqrt(monitor_counts_var))
         logger.info("reduction complete")
 
     def additional_readable_signals(self, dae: "SimpleDae") -> list[Device]:

tests/devices/simpledae/test_reducers.py

@@ -1,5 +1,6 @@
 # pyright: reportMissingParameterType=false
 import math
+import re
 from unittest.mock import AsyncMock
 
 import numpy as np
@@ -9,6 +10,7 @@
 
 from ibex_bluesky_core.devices.simpledae import SimpleDae
 from ibex_bluesky_core.devices.simpledae.reducers import (
+    VARIANCE_ADDITION,
     GoodFramesNormalizer,
     MonitorNormalizer,
     PeriodGoodFramesNormalizer,
@@ -256,6 +258,7 @@ def spectra_bins_easy_to_test() -> sc.DataArray:
         data=sc.Variable(
             dims=["tof"],
             values=[1000.0, 2000.0, 3000.0, 2000.0, 1000.0],
+            variances=[1000.0, 2000.0, 3000.0, 2000.0, 1000.0],
             unit=sc.units.counts,
             dtype="float64",
         ),
@@ -273,6 +276,7 @@ def spectra_bins_tof_convert_to_wavelength_easy() -> sc.DataArray:
         data=sc.Variable(
             dims=["tof"],
             values=[1000.0, 2000.0, 3000.0, 2000.0, 1000.0],
+            variances=[1000.0, 2000.0, 3000.0, 2000.0, 1000.0],
             unit=sc.units.counts,
             dtype="float64",
         ),
@@ -293,6 +297,7 @@ def spectra_bins_tof_convert_to_wavelength_easy_copy() -> sc.DataArray:
         data=sc.Variable(
             dims=["tof"],
             values=[1000.0, 2000.0, 3000.0, 2000.0, 1000.0],
+            variances=[1000.0, 2000.0, 3000.0, 2000.0, 1000.0],
             unit=sc.units.counts,
             dtype="float64",
         ),
@@ -361,13 +366,23 @@ async def test_period_good_frames_normalizer(
 
     period_good_frames_reducer.detectors[1].read_spectrum_dataarray = AsyncMock(
         return_value=sc.DataArray(
-            data=sc.Variable(dims=["tof"], values=[1000.0, 2000.0, 3000.0], unit=sc.units.counts),
+            data=sc.Variable(
+                dims=["tof"],
+                values=[1000.0, 2000.0, 3000.0],
+                variances=[1000.0, 2000.0, 3000.0],
+                unit=sc.units.counts,
+            ),
             coords={"tof": sc.array(dims=["tof"], values=[0, 1, 2, 3])},
         )
     )
     period_good_frames_reducer.detectors[2].read_spectrum_dataarray = AsyncMock(
         return_value=sc.DataArray(
-            data=sc.Variable(dims=["tof"], values=[4000.0, 5000.0, 6000.0], unit=sc.units.counts),
+            data=sc.Variable(
+                dims=["tof"],
+                values=[4000.0, 5000.0, 6000.0],
+                variances=[1000.0, 2000.0, 3000.0],
+                unit=sc.units.counts,
+            ),
             coords={"tof": sc.array(dims=["tof"], values=[0, 1, 2, 3])},
         )
     )
@@ -395,6 +410,7 @@ async def test_period_good_frames_normalizer_uncertainties(
                 values=[1000.0, 2000.0, 3000.0],
                 variances=[1000.0, 2000.0, 3000.0],
                 unit=sc.units.counts,
+                dtype="float64",
             ),
             coords={"tof": sc.array(dims=["tof"], values=[0, 1, 2, 3])},
         )
@@ -406,6 +422,7 @@ async def test_period_good_frames_normalizer_uncertainties(
                 values=[4000.0, 5000.0, 6000.0],
                 variances=[4000.0, 5000.0, 6000.0],
                 unit=sc.units.counts,
+                dtype="float64",
             ),
             coords={
                 "tof": sc.array(
@@ -420,8 +437,10 @@ async def test_period_good_frames_normalizer_uncertainties(
     det_counts_stddev = await period_good_frames_reducer.det_counts_stddev.get_value()
     intensity_stddev = await period_good_frames_reducer.intensity_stddev.get_value()
 
-    assert det_counts_stddev == math.sqrt(21000)
-    assert intensity_stddev == pytest.approx(math.sqrt((21000 + (123**2 / 21000)) / 123**2), 1e-4)
+    assert det_counts_stddev == math.sqrt(21000 + VARIANCE_ADDITION)
+    assert intensity_stddev == pytest.approx(
+        math.sqrt((21000 + VARIANCE_ADDITION) / (123**2)), 1e-4
+    )
 
 
 async def test_period_good_frames_normalizer_zero_counts(
@@ -458,13 +477,11 @@ async def test_period_good_frames_normalizer_zero_counts(
         )
     )
 
-    await period_good_frames_reducer.reduce_data(simpledae)
-
-    det_counts_stddev = await period_good_frames_reducer.det_counts_stddev.get_value()
-    intensity_stddev = await period_good_frames_reducer.intensity_stddev.get_value()
-
-    assert det_counts_stddev == 0
-    assert intensity_stddev == 0
+    with pytest.raises(
+        ValueError,
+        match=re.escape("Cannot normalize; denominator is zero. Check beamline configuration."),
+    ):
+        await period_good_frames_reducer.reduce_data(simpledae)
 
 
 async def test_scalar_normalizer_tof_bounded_zero_to_one_half(
@@ -542,13 +559,23 @@ async def test_monitor_normalizer(
 ):  # 1, 1
     monitor_normalizer.detectors[1].read_spectrum_dataarray = AsyncMock(
         return_value=sc.DataArray(
-            data=sc.Variable(dims=["tof"], values=[1000.0, 2000.0, 3000.0], unit=sc.units.counts),
+            data=sc.Variable(
+                dims=["tof"],
+                values=[1000.0, 2000.0, 3000.0],
+                variances=[1000.0, 2000.0, 3000.0],
+                unit=sc.units.counts,
+            ),
             coords={"tof": sc.array(dims=["tof"], values=[0, 1, 2, 3])},
         )
     )
     monitor_normalizer.monitors[2].read_spectrum_dataarray = AsyncMock(
         return_value=sc.DataArray(
-            data=sc.Variable(dims=["tof"], values=[4000.0, 5000.0, 6000.0], unit=sc.units.counts),
+            data=sc.Variable(
+                dims=["tof"],
+                values=[4000.0, 5000.0, 6000.0],
+                variances=[4000.0, 5000.0, 6000.0],
+                unit=sc.units.counts,
+            ),
             coords={"tof": sc.array(dims=["tof"], values=[0, 1, 2, 3])},
         )
     )
@@ -590,15 +617,11 @@ async def test_monitor_normalizer_zero_counts(
         )
     )
 
-    await monitor_normalizer.reduce_data(simpledae)
-
-    det_counts_stddev = await monitor_normalizer.det_counts_stddev.get_value()
-    mon_counts_stddev = await monitor_normalizer.mon_counts_stddev.get_value()
-    intensity_stddev = await monitor_normalizer.intensity_stddev.get_value()
-
-    assert det_counts_stddev == 0
-    assert mon_counts_stddev == 0
-    assert intensity_stddev == 0
+    with pytest.raises(
+        ValueError,
+        match=re.escape("Cannot normalize; got zero monitor counts. Check beamline configuration."),
+    ):
+        await monitor_normalizer.reduce_data(simpledae)
 
 
 async def test_monitor_normalizer_uncertainties(
@@ -633,9 +656,11 @@ async def test_monitor_normalizer_uncertainties(
     mon_counts_stddev = await monitor_normalizer.mon_counts_stddev.get_value()
     intensity_stddev = await monitor_normalizer.intensity_stddev.get_value()
 
-    assert det_counts_stddev == math.sqrt(6000)
+    assert det_counts_stddev == math.sqrt(6000 + VARIANCE_ADDITION)
     assert mon_counts_stddev == math.sqrt(15000)
-    assert intensity_stddev == pytest.approx(math.sqrt((6000 + (6000**2 / 15000)) / 15000**2), 1e-4)
+    assert intensity_stddev == pytest.approx(
+        (6000 / 15000) * math.sqrt((6000.5 / 6000**2) + (15000 / 15000**2)), 1e-8
+    )
 
 
 def test_monitor_normalizer_publishes_raw_and_normalized_counts(
@@ -658,7 +683,7 @@ def test_monitor_normalizer_publishes_raw_and_normalized_count_uncertainties(
     assert monitor_normalizer.mon_counts_stddev in readables
 
 
-async def test_monitor_normalizer_det_sum_normal_mon_sum_tof_bound_(  # 1, 2
+async def test_monitor_normalizer_det_sum_normal_mon_sum_tof_bound(  # 1, 2
     simpledae: SimpleDae,
     monitor_normalizer_zero_to_one_half_det_norm_mon_tof: MonitorNormalizer,
     spectra_bins_easy_to_test: sc.DataArray,

tests/devices/test_dae.py

@@ -12,20 +12,23 @@
 from ophyd_async.testing import get_mock_put, set_mock_value
 
 from ibex_bluesky_core.devices import compress_and_hex, dehex_and_decompress
+from ibex_bluesky_core.devices.dae import convert_xml_to_names_and_values, set_value_in_dae_xml
 from ibex_bluesky_core.devices.dae.dae import Dae, RunstateEnum
+from ibex_bluesky_core.devices.dae.dae_controls import BeginRunExBits
 from ibex_bluesky_core.devices.dae.dae_period_settings import (
     DaePeriodSettings,
     DaePeriodSettingsData,
     PeriodSource,
     PeriodType,
     SinglePeriodSettings,
+    _convert_period_settings_to_xml,
 )
 from ibex_bluesky_core.devices.dae.dae_settings import (
     DaeSettings,
     DaeSettingsData,
     TimingSource,
 )
-from ibex_bluesky_core.devices.dae.dae_spectra import VARIANCE_ADDITION, DaeSpectra
+from ibex_bluesky_core.devices.dae.dae_spectra import DaeSpectra
 from ibex_bluesky_core.devices.dae.dae_tcb_settings import (
     CalculationMethod,
     DaeTCBSettings,
@@ -34,11 +37,8 @@
     TimeRegimeMode,
     TimeRegimeRow,
     TimeUnit,
+    _convert_tcb_settings_to_xml,
 )
-from src.ibex_bluesky_core.devices.dae import convert_xml_to_names_and_values, set_value_in_dae_xml
-from src.ibex_bluesky_core.devices.dae.dae_controls import BeginRunExBits
-from src.ibex_bluesky_core.devices.dae.dae_period_settings import _convert_period_settings_to_xml
-from src.ibex_bluesky_core.devices.dae.dae_tcb_settings import _convert_tcb_settings_to_xml
 from tests.conftest import MOCK_PREFIX
 from tests.devices.dae_testing_data import (
     dae_settings_template,
@@ -980,9 +980,9 @@ async def test_read_spectrum_dataarray(spectrum: DaeSpectra):
                 dims=["tof"],
                 values=[1000, 2000, 3000],
                 variances=[
-                    1000 + VARIANCE_ADDITION,
-                    2000 + VARIANCE_ADDITION,
-                    3000 + VARIANCE_ADDITION,
+                    1000,
+                    2000,
+                    3000,
                 ],
                 unit=sc.units.counts,
                 dtype="float32",