Polarisation / spin-echo tuning scans

doc/devices/dae.md

@@ -295,6 +295,29 @@ 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
+
+Polarization refers to the property of transverse waves which specifies the geometrical orientation of the 
+oscillations. 
+
+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.
+
+{py:obj}`ibex_bluesky_core.devices.simpledae.polarisingdae.polarization`
+
+See [`PolarisingDae`](#PolarisingDae) and [`PolarisingReducer`](#PolarisingReducer) for how this is integrated into DAE behaviour. 
+
 ## Waiters
 
 A [`waiter`](ibex_bluesky_core.devices.simpledae.Waiter) defines an arbitrary strategy for how long to count at each point.
@@ -347,6 +370,62 @@ Waits for a user-specified time duration, irrespective of DAE state.
 
 Does not publish any additional signals.
 
+## Polarising DAE
+
+The polarising DAE provides specialised functionality for taking data whilst taking into account the polarity of the beam.
+
+### PolarisingDae
+
+[`PolarisingDae`](ibex_bluesky_core.devices.simpledae.polarisingdae.PolarisingDae) is a more complex version of `SimpleDae`, designed specifically for taking polarisation measurements. It requires a flipper device and uses it to flip from one neutron state to the other between runs.
+
+Key features:
+- Controls a flipper device to switch between neutron states
+- Handles two separate reduction strategies (up and down states)
+- Calculates polarisation from the two states after the two runs
+
+### polarising_dae
+
+[`polarising_dae`](ibex_bluesky_core.devices.simpledae.polarisingdae.polarising_dae) is a helper function that creates a configured `PolarisingDae` instance with wavelength binning based normalisation and polarisation calculation capabilities.
+
+The following is how you may want to use `polarising_dae`:
+```python
+import scipp
+
+flipper = block_rw(float, "alice")
+wavelength_interval = scipp.array(dims=["tof"], values=[0, 9999999999.0], unit=scipp.units.angstrom, dtype="float64") # Creates a wavelength interval of the whole sprectrum
+total_flight_path_length = sc.scalar(value=10, unit=sc.units.m)
+
+dae = polarising_dae(det_pixels=[1], frames=500, flipper=flipper, flipper_states=(0.0, 1.0), intervals=[wavelength_interval], total_flight_path_length=total_flight_path_length, monitor=2)
+```
+
+:::{note}
+  Notice how you must define what the `flipper_states` are to the polarising dae. This is so that it knows what to assign to the `flipper` device to move it to the "up state" and "down state"
+  .
+:::
+
+### Polarising Reducers
+
+#### WavelengthBoundedNormalizer
+
+[`WavelengthBoundedNormalizer`](ibex_bluesky_core.devices.simpledae.polarisingdae.WavelengthBoundedNormalizer) sums wavelength-bounded spectra and normalises by monitor intensity.
+
+Published signals:
+- `wavelength_bands`: DeviceVector containing wavelength band measurements
+  - `det_counts`: detector counts in the wavelength band
+  - `mon_counts`: monitor counts in the wavelength band
+  - `intensity`: normalised intensity in the wavelength band
+  - Associated uncertainty measurements for each value
+
+#### PolarisingReducer
+
+[`PolarisingReducer`](ibex_bluesky_core.devices.simpledae.polarisingdae.PolarisingReducer) calculates polarisation from 'spin-up' and 'spin-down' states of a polarising DAE. Uses the [`polarization`](#polarizationasymmetry) algorithm.
+
+Published signals:
+- `wavelength_bands`: DeviceVector containing polarisation measurements
+  - `polarisation`: The calculated polarisation value for that wavelength band
+  - `polarisation_ratio`: Ratio between up and down states for that wavelength band
+  - Associated uncertainty measurements for each value
+
 ---
 
 ## `Dae` (base class, advanced)
@@ -465,9 +544,6 @@ A [`DaeSpectra`](ibex_bluesky_core.devices.dae.DaeSpectra) object provides 3 arr
 The [`Dae`](ibex_bluesky_core.devices.dae) base class does not provide any spectra by default. User-level classes should specify 
 the set of spectra which they are interested in.
 
-
-
-
 Spectra can be summed between two bounds based on time of flight bounds, or wavelength bounds, for both detector and monitor normalizers.
 
 Both Scalar Normalizers (PeriodGoodFramesNormalizer, GoodFramesNormalizer) and MonitorNormalizers
@@ -478,3 +554,25 @@ or wavelength bounds.
 or sums using time of flight bounds, or wavelength bounds.
 
 For both options, the default, if none is specified, is to use pre-existing bounds.
+
+### Wavelength Band Classes
+
+#### WavelengthBand
+
+[`WavelengthBand`](ibex_bluesky_core.devices.dae.WavelengthBand) represents a few measurements within a specific wavelength band. Has a setter method to assign values to the published signals.
+
+Additional Signals:
+- `det_counts`: Detector counts
+- `mon_counts`: Monitor counts
+- `intensity`: Normalized intensity
+- Associated uncertainty measurements for each value
+
+#### PolarisedWavelengthBand
+
+[`PolarisedWavelengthBand`](ibex_bluesky_core.devices.dae.PolarisedWavelengthBand) represents the polarisation information calculated using measurements taken from two `WavelengthBand` objects, one published from an "up state" `WavelengthBoundedNormalizer`
+and the other from a "down state" `WavelengthBoundedNormalizer`. Has a setter method to assign values to the published signals.
+
+Additional signals:
+- `polarisation`: The calculated polarisation value
+- `polarisation_ratio`: Ratio between up and down states
+- Associated uncertainty measurements for polarisation values

src/ibex_bluesky_core/callbacks/__init__.py

@@ -20,7 +20,7 @@
 from ibex_bluesky_core.callbacks._file_logger import (
     HumanReadableFileCallback,
 )
-from ibex_bluesky_core.callbacks._fitting import LiveFit, LiveFitLogger
+from ibex_bluesky_core.callbacks._fitting import ChainedLiveFit, LiveFit, LiveFitLogger
 from ibex_bluesky_core.callbacks._plotting import LivePColorMesh, LivePlot, PlotPNGSaver, show_plot
 from ibex_bluesky_core.callbacks._utils import get_default_output_path
 from ibex_bluesky_core.fitting import FitMethod
@@ -32,6 +32,7 @@
 
 
 __all__ = [
+    "ChainedLiveFit",
     "DocLoggingCallback",
     "HumanReadableFileCallback",
     "ISISCallbacks",

src/ibex_bluesky_core/callbacks/_fitting.py

@@ -4,13 +4,18 @@
 import logging
 import os
 import warnings
+from itertools import zip_longest
 from pathlib import Path
 
+import lmfit
 import numpy as np
-from bluesky.callbacks import CallbackBase
+from bluesky.callbacks import CallbackBase, LiveFitPlot
 from bluesky.callbacks import LiveFit as _DefaultLiveFit
 from bluesky.callbacks.core import make_class_safe
-from event_model import Event, RunStart, RunStop
+from event_model import Event, EventDescriptor, RunStart, RunStop
+from lmfit import Parameter
+from matplotlib.axes import Axes
+from numpy import typing as npt
 
 from ibex_bluesky_core.callbacks._utils import (
     DATA,
@@ -24,7 +29,7 @@
 
 logger = logging.getLogger(__name__)
 
-__all__ = ["LiveFit", "LiveFitLogger"]
+__all__ = ["ChainedLiveFit", "LiveFit", "LiveFitLogger"]
 
 
 @make_class_safe(logger=logger)  # pyright: ignore (pyright doesn't understand this decorator)
@@ -84,20 +89,17 @@ def update_weight(self, weight: float | None = 0.0) -> None:
         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.
-
-        Args:
-            None
-
-        Returns:
-            None
-
-        """
+    def can_fit(self) -> bool:
+        """Check if enough data points have been collected to fit."""
         n = len(self.model.param_names)
-        if len(self.ydata) < n:
+        return len(self.ydata) >= n
+
+    def update_fit(self) -> None:
+        """Use the guess function with the most recent x and y values after every update."""
+        if not self.can_fit():
             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
+                {len(self.model.param_names)} data points""",
                 stacklevel=1,
             )
         else:
@@ -174,7 +176,7 @@ def start(self, doc: RunStart) -> None:
         self.filename = self.output_dir / f"{rb_num}" / file
 
     def event(self, doc: Event) -> Event:
-        """Start collecting, y, x and yerr data.
+        """Start collecting, y, x, and yerr data.
 
         Args:
             doc: (Event): An event document.
@@ -251,3 +253,138 @@ def write_fields_table_uncertainty(self) -> None:
 
         rows = zip(self.x_data, self.y_data, self.yerr_data, self.y_fit_data, strict=True)
         self.csvwriter.writerows(rows)
+
+
+class ChainedLiveFit(CallbackBase):
+    """Processes multiple LiveFits, each fit's results inform the next, with optional plotting.
+
+    This callback handles a sequence of LiveFit instances where the parameters from each
+    completed fit serve as the initial guess for the subsequent fit. Optional plotting
+    is built in using LivePlotLib.
+    """
+
+    def __init__(
+        self,
+        method: FitMethod,
+        y: list[str],
+        x: str,
+        *,
+        yerr: list[str] | None = None,
+        ax: list[Axes] | None = None,
+    ) -> None:
+        """Initialise ChainedLiveFit with multiple LiveFits.
+
+        Args:
+            method: FitMethod instance for fitting
+            y: List of y-axis variable names
+            x: x-axis variable name
+            yerr: Optional list of error values corresponding to y variables
+            ax: A list of axes to plot fits on to. Creates LiveFitPlot instances.
+
+        """
+        super().__init__()
+
+        self._livefits = [
+            LiveFit(method=method, y=y_name, x=x, yerr=yerr_name)
+            for y_name, yerr_name in zip_longest(y, yerr or [])
+        ]
+
+        self._livefitplots = [
+            LiveFitPlot(livefit=livefit, ax=axis)
+            for livefit, axis in zip(self._livefits, ax or [], strict=False)
+        ]
+
+    def _process_doc(
+        self, doc: RunStart | Event | RunStop | EventDescriptor, method_name: str
+    ) -> None:
+        """Process a document for either LivePlots or LiveFits.
+
+        Args:
+            doc: document to process
+            method_name: Name of the method to call ('start', 'descriptor', 'event', or 'stop')
+
+        """
+        callbacks = self._livefitplots or self._livefits
+        for callback in callbacks:
+            getattr(callback, method_name)(doc)
+
+    def start(self, doc: RunStart) -> None:
+        """Process start document for all callbacks.
+
+        Args:
+            doc: RunStart document
+
+        """
+        self._process_doc(doc, "start")
+
+    def descriptor(self, doc: EventDescriptor) -> None:
+        """Process descriptor document for all callbacks.
+
+        Args:
+            doc: EventDescriptor document.
+
+        """
+        self._process_doc(doc, "descriptor")
+
+    def event(self, doc: Event) -> Event:
+        """Process event document for all callbacks.
+
+        Args:
+            doc: Event document
+
+        """
+        init_guess = {}
+
+        for livefit in self._livefits:
+            rem_guess = livefit.method.guess
+            try:
+                if init_guess:
+                    # Use previous fit results as initial guess for next fit
+                    def guess_func(
+                        a: npt.NDArray[np.float64], b: npt.NDArray[np.float64]
+                    ) -> dict[str, lmfit.Parameter]:
+                        return {
+                            name: Parameter(name, value.value)
+                            for name, value in init_guess.items()  # noqa: B023
+                        }
+
+                    # Using value.value means that paramater uncertainty
+                    # is not carried over between fits
+                    livefit.method.guess = guess_func
+
+                if self._livefitplots:
+                    self._livefitplots[self._livefits.index(livefit)].event(doc)
+                else:
+                    livefit.event(doc)
+
+            finally:
+                livefit.method.guess = rem_guess
+
+                if livefit.can_fit():
+                    if livefit.result is None:
+                        raise RuntimeError("LiveFit.result was None. Could not update fit.")
+
+                    init_guess = livefit.result.params
+
+        from ibex_bluesky_core.callbacks import show_plot  # noqa: PLC0415
+
+        show_plot()
+
+        return doc
+
+    def stop(self, doc: RunStop) -> None:
+        """Process stop document and update fitting parameters.
+
+        Args:
+            doc: RunStop document
+
+        """
+        self._process_doc(doc, "stop")
+
+    def get_livefits(self) -> list[LiveFit]:
+        """Return a list of the livefits."""
+        return self._livefits
+
+    def get_livefitplots(self) -> list[LiveFitPlot]:
+        """Return a list of the livefitplots."""
+        return self._livefitplots

src/ibex_bluesky_core/devices/dae/__init__.py

@@ -33,7 +33,11 @@
     SinglePeriodSettings,
 )
 from ibex_bluesky_core.devices.dae._settings import DaeSettings, DaeSettingsData, DaeTimingSource
-from ibex_bluesky_core.devices.dae._spectra import DaeSpectra
+from ibex_bluesky_core.devices.dae._spectra import (
+    DaeSpectra,
+    PolarisedWavelengthBand,
+    WavelengthBand,
+)
 from ibex_bluesky_core.devices.dae._tcb_settings import (
     DaeTCBSettings,
     DaeTCBSettingsData,
@@ -63,13 +67,15 @@
     "DaeTimingSource",
     "PeriodSource",
     "PeriodType",
+    "PolarisedWavelengthBand",
     "RunstateEnum",
     "SinglePeriodSettings",
     "TCBCalculationMethod",
     "TCBTimeUnit",
     "TimeRegime",
     "TimeRegimeMode",
     "TimeRegimeRow",
+    "WavelengthBand",
 ]
 
 T = TypeVar("T", bound=SignalDatatype)