Merge pull request #97 from ISISComputingGroup/58_polarization_asymmetry

Add function which calculates Polarization & propagating uncertainties

Author: jackbdoughty

doc/devices/dae.md

@@ -295,6 +295,39 @@ as a list, and `units` (μs/microseconds for time of flight bounding, and angstr
 
 If you don't specify either of these options, they will default to summing over the entire spectrum.
 
+### Polarization/Asymmetry
+
+ibex_bluesky_core provides a helper method,
+{py:obj}`ibex_bluesky_core.devices.simpledae.polarization`, for calculating the quantity (a-b)/(a+b). This quantity is used, for example, in neutron polarization measurements, and in calculating asymmetry for muon measurements.
+
+For this expression, scipp's default uncertainty propagation rules cannot be used as the uncertainties on (a-b) are correlated with those of (a+b) in the division step - but scipp assumes uncorrelated data. This helper method calculates the uncertainties following linear error propagation theory, using the partial derivatives of the above expression.
+
+The partial derivatives are:
+
+$ \frac{\delta}{\delta a}\Big(\frac{a - b}{a + b}) = \frac{2b}{(a+b)^2} $
+
+$ \frac{\delta}{\delta b}\Big(\frac{a - b}{a + b}) = -\frac{2a}{(a+b)^2} $
+
+
+Which then means the variances computed by this helper function are:
+
+$ Variance = (\frac{\delta}{\delta a}^2 * variance_a) + (\frac{\delta}{\delta b}^2 * variance_b)  $ 
+
+
+The polarization funtion provided will calculate the polarization between two values, A and B, which 
+have different definitions based on the instrument context.
+
+Instrument-Specific Interpretations
+SANS Instruments (e.g., LARMOR)
+A: Intensity in DAE period before switching a flipper.
+B: Intensity in DAE period after switching a flipper.
+
+Reflectometry Instruments (e.g., POLREF)
+Similar to LARMOR, A and B represent intensities before and after flipper switching.
+
+Muon Instruments
+A and B refer to Measurements from different detector banks.
+
 ## Waiters
 
 A [`waiter`](ibex_bluesky_core.devices.simpledae.Waiter) defines an arbitrary strategy for how long to count at each point.

pyproject.toml

@@ -69,6 +69,7 @@ dev = [
   "pytest-cov",
   "pytest-env",
   "pyqt6",  # For dev testing with matplotlib's qt backend.
+  "uncertainties" # For dev testing polarization, variance & uncertainties
 ]
 
 [project.urls]

src/ibex_bluesky_core/devices/simpledae/__init__.py

@@ -23,6 +23,8 @@
     PeriodGoodFramesNormalizer,
     PeriodSpecIntegralsReducer,
     ScalarNormalizer,
+    polarization,
+    sum_spectra,
     tof_bounded_spectra,
     wavelength_bounded_spectra,
 )
@@ -67,6 +69,8 @@
     "Waiter",
     "check_dae_strategies",
     "monitor_normalising_dae",
+    "polarization",
+    "sum_spectra",
     "tof_bounded_spectra",
     "wavelength_bounded_spectra",
 ]

src/ibex_bluesky_core/devices/simpledae/_reducers.py

@@ -37,10 +37,14 @@
 async def sum_spectra(spectra: Collection[DaeSpectra]) -> sc.Variable | sc.DataArray:
     """Read and sum a number of spectra from the DAE.
 
-    Returns a scipp scalar, which has .value and .variance properties for accessing the sum
-    and variance respectively of the summed counts.
+    Args:
+        spectra: a Collection type object of DAE spectra
+
+    Returns:
+        scipp :external+scipp:py:obj:`Variable <scipp.Variable>`
+        or :external+scipp:py:obj:`DataArray <scipp.DataArray>` describing
+        the sum of the provided spectra.
 
-    More info on scipp scalars can be found here: https://scipp.github.io/generated/functions/scipp.scalar.html
     """
     logger.info("Summing %d spectra using scipp", len(spectra))
     summed_counts = sc.scalar(value=0, unit=sc.units.counts, dtype="float64")
@@ -56,13 +60,14 @@ def tof_bounded_spectra(
     """Sum a set of neutron spectra between the specified time of flight bounds.
 
     Args:
-        bounds: A scipp array of size 2, no variances, unit of us,
-            where the second element must be larger than the first.
+        bounds: A scipp :external+scipp:py:obj:`array <scipp.array>` of size 2, no variances, unit
+            of us, where the second element must be larger than the first.
 
-    Returns a scipp scalar, which has .value and .variance properties for accessing the sum
-    and variance respectively of the summed counts.
+    :rtype:
+        scipp :external+scipp:py:obj:`scalar <scipp.scalar>`
 
-    More info on scipp arrays and scalars can be found here: https://scipp.github.io/generated/functions/scipp.scalar.html
+    Returns a scipp :external+scipp:py:obj:`scalar <scipp.scalar>`, which has .value and .variance
+    properties for accessing the sum and variance respectively of the summed counts.
 
     """
     bounds_value = 2
@@ -88,17 +93,22 @@ def wavelength_bounded_spectra(
     """Sum a set of neutron spectra between the specified wavelength bounds.
 
     Args:
-        bounds: A scipp array of size 2 of wavelength bounds, in units of angstrom,
-            where the second element must be larger than the first.
-        total_flight_path_length: A scipp scalar of Ltotal (total flight path length), the path
-            length from neutron source to detector or monitor, in units of meters.
+        bounds: A scipp :external+scipp:py:obj:`array <scipp.array>` of size 2 of wavelength bounds,
+            in units of angstrom, where the second element must be larger than the first.
+        total_flight_path_length: A scipp :external+scipp:py:obj:`scalar <scipp.scalar>` of Ltotal
+            (total flight path length), the path length from neutron source to detector or monitor,
+            in units of meters.
+
+    :rtype:
+        scipp :external+scipp:py:obj:`scalar <scipp.scalar>`
 
     Time of flight is converted to wavelength using scipp neutron's library function
-        `wavelength_from_tof`, more info on which can be found here:
-        https://scipp.github.io/scippneutron/generated/modules/scippneutron.conversion.tof.wavelength_from_tof.html
+    `wavelength_from_tof`, more info on which can be found here:
+    :external+scippneutron:py:obj:`wavelength_from_tof
+    <scippneutron.conversion.tof.wavelength_from_tof>`
 
-    Returns a scipp scalar, which has .value and .variance properties for accessing the sum
-    and variance respectively of the summed counts.
+    Returns a scipp :external+scipp:py:obj:`scalar <scipp.scalar>`, which has .value and .variance
+    properties for accessing the sum and variance respectively of the summed counts.
 
     """
     bounds_value = 2
@@ -125,6 +135,60 @@ async def sum_spectra_with_wavelength(
     return sum_spectra_with_wavelength
 
 
+def polarization(
+    a: sc.Variable | sc.DataArray, b: sc.Variable | sc.DataArray
+) -> sc.Variable | sc.DataArray:
+    """Calculate polarization value and propagate uncertainties.
+
+    This function computes the polarization given by the formula (a-b)/(a+b)
+    and propagates the uncertainties associated with a and b.
+
+    Args:
+        a: scipp :external+scipp:py:obj:`Variable <scipp.Variable>`
+            or :external+scipp:py:obj:`DataArray <scipp.DataArray>`
+        b: scipp :external+scipp:py:obj:`Variable <scipp.Variable>`
+            or :external+scipp:py:obj:`DataArray <scipp.DataArray>`
+
+    Returns:
+        polarization, ``(a - b) / (a + b)``, as a scipp
+        :external+scipp:py:obj:`Variable <scipp.Variable>`
+        or :external+scipp:py:obj:`DataArray <scipp.DataArray>`
+
+    On SANS instruments e.g. LARMOR, A and B correspond to intensity in different DAE
+    periods (before/after switching a flipper) and the output is interpreted as a neutron
+    polarization ratio.
+
+    On reflectometry instruments e.g. POLREF, the situation is the same as on LARMOR.
+
+    On muon instruments, A and B correspond to measuring from forward/backward detector
+    banks, and the output is interpreted as a muon asymmetry.
+
+    """
+    if a.unit != b.unit:
+        raise ValueError("The units of a and b are not equivalent.")
+    if a.sizes != b.sizes:
+        raise ValueError("Dimensions/shape of a and b must match.")
+
+    # This line allows for dims, units, and dtype to be handled by scipp
+    polarization_value = (a - b) / (a + b)
+
+    variances_a = a.variances
+    variances_b = b.variances
+    values_a = a.values
+    values_b = b.values
+
+    # Calculate partial derivatives
+    partial_a = 2 * values_b / (values_a + values_b) ** 2
+    partial_b = -2 * values_a / (values_a + values_b) ** 2
+
+    variance_return = (partial_a**2 * variances_a) + (partial_b**2 * variances_b)
+
+    # Propagate uncertainties
+    polarization_value.variances = variance_return
+
+    return polarization_value
+
+
 class ScalarNormalizer(Reducer, StandardReadable, ABC):
     """Sum a set of user-specified spectra, then normalize by a scalar signal."""
 

tests/devices/simpledae/test_reducers.py

@@ -7,6 +7,7 @@
 import pytest
 import scipp as sc
 from ophyd_async.testing import get_mock_put, set_mock_value
+from uncertainties import ufloat, unumpy
 
 from ibex_bluesky_core.devices.simpledae import (
     VARIANCE_ADDITION,
@@ -19,7 +20,7 @@
     tof_bounded_spectra,
     wavelength_bounded_spectra,
 )
-from ibex_bluesky_core.devices.simpledae._reducers import DSpacingMappingReducer
+from ibex_bluesky_core.devices.simpledae._reducers import DSpacingMappingReducer, polarization
 
 
 @pytest.fixture
@@ -965,6 +966,93 @@ def test_wavelength_bounded_spectra_bounds_missing_or_too_many(data: list[float]
         )
 
 
+# Polarization
+@pytest.mark.parametrize(
+    ("a", "b", "variance_a", "variance_b"),
+    [
+        # Case 1: Symmetric case with equal uncertainties
+        (5.0, 3.0, 0.1, 0.1),
+        # Case 2: Asymmetric case with different uncertainties
+        (10.0, 6.0, 0.2, 0.3),
+        # Case 3: Case with larger values and different uncertainty magnitudes
+        (100.0, 60.0, 1.0, 2.0),
+    ],
+)
+def test_polarization_function_calculates_accurately(a, b, variance_a, variance_b):
+    # 'Uncertainties' library ufloat type; a nominal value and an error value
+    a_ufloat = ufloat(a, variance_a)
+    b_ufloat = ufloat(b, variance_b)
+
+    # polarization value, i.e. (a - b) / (a + b)
+    polarization_ufloat = (a_ufloat.n - b_ufloat.n) / (a_ufloat.n + b_ufloat.n)
+
+    # the partial derivatives of a and b, calculated with 'uncertainties' library's ufloat type
+    partial_a = (2 * b_ufloat.n) / ((a_ufloat.n + b_ufloat.n) ** 2)
+    partial_b = (-2 * a_ufloat.n) / ((a_ufloat.n + b_ufloat.n) ** 2)
+
+    # variance calculated with 'uncertainties' library
+    variance = (partial_a**2 * a_ufloat.s) + (partial_b**2 * b_ufloat.s)
+    uncertainty = variance**0.5  # uncertainty is sqrt of variance
+
+    # Two scipp scalars, to test our polarization function
+    var_a = sc.scalar(value=a, variance=variance_a, unit="", dtype="float64")
+    var_b = sc.scalar(value=b, variance=variance_b, unit="", dtype="float64")
+    result_value = polarization(var_a, var_b)
+    result_uncertainy = (result_value.variance) ** 0.5  # uncertainty is sqrt of variance
+
+    assert result_value.value == pytest.approx(polarization_ufloat)
+    assert result_uncertainy == pytest.approx(uncertainty)
+
+
+# test that arrays are supported
+@pytest.mark.parametrize(
+    ("a", "b", "variances_a", "variances_b"),
+    [
+        ([5.0, 10.0, 100.0], [3.0, 6.0, 60.0], [0.1, 0.2, 1.0], [0.1, 0.3, 2.0]),
+    ],
+)
+def test_polarization_2_arrays(a, b, variances_a, variances_b):
+    # 'Uncertainties' library ufloat type; a nominal value and an error value
+
+    a_arr = unumpy.uarray(a, [v**0.5 for v in variances_a])  # convert variances to std dev
+    b_arr = unumpy.uarray(b, [v**0.5 for v in variances_b])
+
+    # polarization value, i.e. (a - b) / (a + b)
+    polarization_ufloat = (a_arr - b_arr) / (a_arr + b_arr)
+
+    var_a = sc.array(dims="x", values=a, variances=variances_a, unit="", dtype="float64")
+    var_b = sc.array(dims="x", values=b, variances=variances_b, unit="", dtype="float64")
+
+    result_value = polarization(var_a, var_b)
+
+    result_uncertainties = (result_value.variances) ** 0.5
+
+    assert result_value.values == pytest.approx(unumpy.nominal_values(polarization_ufloat))
+    assert result_uncertainties == pytest.approx(unumpy.std_devs(polarization_ufloat))
+
+
+# test that units don't match
+def test_polarization_units_mismatch():
+    var_a = sc.scalar(value=1, variance=0.1, unit="m", dtype="float64")
+    var_b = sc.scalar(value=1, variance=0.1, unit="u", dtype="float64")
+
+    with pytest.raises(
+        expected_exception=ValueError, match=r"The units of a and b are not equivalent."
+    ):
+        polarization(var_a, var_b)
+
+
+# test that arrays are of unmatching sizes
+def test_polarization_arrays_of_different_sizes():
+    var_a = sc.array(dims=["x"], values=[1, 2], variances=[0.1, 0.1], unit="m", dtype="float64")
+    var_b = sc.array(dims=["x"], values=[1], variances=[0.1], unit="m", dtype="float64")
+
+    with pytest.raises(
+        expected_exception=ValueError, match=r"Dimensions/shape of a and b must match."
+    ):
+        polarization(var_a, var_b)
+
+
 @pytest.mark.parametrize(
     ("current_period", "mon_integrals", "det_integrals"),
     [