Polarisation / spin-echo tuning scans

doc/devices/dae.md

@@ -13,7 +13,7 @@ This means that [`SimpleDae`](ibex_bluesky_core.devices.simpledae.SimpleDae) is
 example running using either one DAE run per scan point, or one DAE period per scan point.
 
 For complex use-cases, particularly those where the DAE may need to start and stop multiple 
-acquisitions per scan point (e.g. polarization measurements), [`SimpleDae`](ibex_bluesky_core.devices.simpledae.SimpleDae) is unlikely to be 
+acquisitions per scan point (e.g. Polarisation measurements), [`SimpleDae`](ibex_bluesky_core.devices.simpledae.SimpleDae) is unlikely to be 
 suitable; instead the [`Dae`](ibex_bluesky_core.devices.dae.Dae) class should be subclassed directly to allow for finer control.
 
 ## Example configurations
@@ -295,10 +295,10 @@ 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
+### Polarisation/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.
+{py:obj}`ibex_bluesky_core.utils.calculate_polarisation`, for calculating the quantity (a-b)/(a+b). This quantity is used, for example, in neutron polarisation 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.
 
@@ -314,7 +314,7 @@ 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 
+The polarisation function provided will calculate the polarisation between two values, A and B, which 
 have different definitions based on the instrument context.
 
 Instrument-Specific Interpretations
@@ -328,6 +328,10 @@ Similar to LARMOR, A and B represent intensities before and after flipper switch
 Muon Instruments
 A and B refer to Measurements from different detector banks.
 
+{py:obj}`ibex_bluesky_core.utils.calculate_polarisation`
+
+See [`PolarisationReducer`](#PolarisationReducer) 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.
@@ -380,6 +384,63 @@ 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.
+
+### DualRunDae
+
+[`DualRunDae`](ibex_bluesky_core.devices.polarisingdae.DualRunDae) is a more complex version of [`SimpleDae`](ibex_bluesky_core.devices.simpledae.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 three separate reduction strategies 
+  - Up & Down Reducers, which run after each run
+  - Main reducer, which runs after everything else
+
+### polarising_dae
+
+[`polarising_dae`](ibex_bluesky_core.devices.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
+
+#### MultiWavelengthBandNormalizer
+
+[`MultiWavelengthBandNormalizer`](ibex_bluesky_core.devices.polarisingdae.MultiWavelengthBandNormalizer) 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
+
+[`PolarisationReducer`](ibex_bluesky_core.devices.polarisingdae.PolarisationReducer) calculates polarisation from 'spin-up' and 'spin-down' states of a polarising DAE. Uses the [`Polarisation`](#polarisationasymmetry) 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)
@@ -498,9 +559,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

doc/fitting/fitting.md

@@ -1,4 +1,4 @@
-# Fitting Callback
+# Fitting Callbacks
 
 Similar to [`LivePlot`](../callbacks/plotting.md), [`ibex_bluesky_core`](ibex_bluesky_core) provides a thin wrapper around Bluesky's [`LiveFit`](ibex_bluesky_core.callbacks.LiveFit) class, enhancing it with additional functionality to better support real-time data fitting. This wrapper not only offers a wide selection of models to fit your data on, but also introduces guess generation for fit parameters. As new data points are acquired, the wrapper refines these guesses dynamically, improving the accuracy of the fit with each additional piece of data, allowing for more efficient and adaptive real-time fitting workflows.
 
@@ -7,7 +7,8 @@ In order to use the wrapper, import[`LiveFit`](ibex_bluesky_core.callbacks.LiveF
 ```py
 from ibex_bluesky_core.callbacks.fitting import LiveFit
 ```
-**Note:** that you do not *need* [`LivePlot`](ibex_bluesky_core.callbacks.LivePlot)  for [`LiveFit`](ibex_bluesky_core.callbacks.LiveFit) to work but it may be useful to know visaully how well the model fits to the raw data.
+.. note::
+  that you do not *need* [`LivePlot`](ibex_bluesky_core.callbacks.LivePlot)  for [`LiveFit`](ibex_bluesky_core.callbacks.LiveFit) to work but it may be useful to know visaully how well the model fits to the raw data.
 
 ## Configuration
 
@@ -31,7 +32,8 @@ fit_callback = LiveFit(Gaussian.fit(), y="y_signal", x="x_signal", yerr="yerr_si
 fit_plot_callback = LiveFitPlot(fit_callback, ax=ax, color="r")
 ```
 
-**Note:** that the [`LiveFit`](ibex_bluesky_core.callbacks.LiveFit) callback doesn't directly do the plotting, it will return function parameters of the model its trying to fit to; a [`LiveFit`](ibex_bluesky_core.callbacks.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.
+.. note::
+  that the [`LiveFit`](ibex_bluesky_core.callbacks.LiveFit) callback doesn't directly do the plotting, it will return function parameters of the model its trying to fit to; a [`LiveFit`](ibex_bluesky_core.callbacks.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. 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`](ibex_bluesky_core.callbacks.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.
 
@@ -89,7 +91,8 @@ lf = LiveFit([FIT].fit(), y="y_signal", x="x_signal", update_every=0.5)
 
 The `[FIT].fit()` function will pass the [`FitMethod`](ibex_bluesky_core.fitting.FitMethod) object straight to the [`LiveFit`](ibex_bluesky_core.callbacks.LiveFit) class.
 
-**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:
+.. 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_signal", x="x_signal", update_every=0.5)
@@ -146,7 +149,8 @@ lf = LiveFit(fit_method, y="y_signal", x="x_signal", update_every=0.5)
 # Then subscribe to LiveFitPlot(lf, ...)
 ```
 
-**Note:** that the parameters returned from the guess function must allocate to the arguments to the model function, ignoring the independant variable e.g `x` in this case. Array-like structures are not allowed. See the [lmfit documentation](https://lmfit.github.io/lmfit-py/parameters.html) for more information.
+.. note::
+  that the parameters returned from the guess function must allocate to the arguments to the model function, ignoring the independant variable e.g `x` in this case. Array-like structures are not allowed. See the [lmfit documentation](https://lmfit.github.io/lmfit-py/parameters.html) for more information.
 
 #### Option 2: Continued
 
@@ -201,10 +205,60 @@ lf = LiveFit(fit_method, y="y_signal", x="x_signal", update_every=0.5)
 
 Or you can create a completely user-defined fitting method.
 
-**Note:** that for fits that require arguments, you will need to pass values to their respecitive `.model` and `.guess` functions. E.g for `Polynomial` fitting:
+.. note::
+  that for fits that require arguments, you will need to pass values to their respecitive `.model` and `.guess` functions. E.g for `Polynomial` fitting:
 
 ```py
 fit_method = FitMethod(Polynomial.model(3), different_guess) # If using a custom guess function
 lf = LiveFit(fit_method, ...)
 ```
 See the [standard fits](#models) list above for standard fits which require parameters. It gets more complicated if you want to define your own custom model or guess which you want to pass parameters to. You will have to define a function that takes these parameters and returns the model / guess function with the subsituted values.
+
+## Chained Fitting
+
+[`ChainedLiveFit`](ibex_bluesky_core.callbacks.ChainedLiveFit) is a specialised callback that manages multiple LiveFit instances in a chain, where each fit's results inform the next fit's initial parameters. This is particularly useful when dealing with complex data sets where subsequent fits depend on the parameters obtained from previous fits.
+
+This is useful for when you need to be careful with your curve fitting due to the presence of noisy data. It allows you to fit your widest (full) wavelength band first and then using its fit parameters as the initial guess of the parameters for the next fit
+
+# Usage
+To show how we expect this to be used we will use the PolarisingDae and wavelength bands to highlight the need for the carry over of fitting parameters. Below shows two wavelength bands, first bigger than the second, we will fit to the data in the first and carry it over to the data in the second to improve its, otherwise worse, fit.
+
+```python
+# Needed for PolarisingDae
+flipper = block_rw(float, "flipper")
+total_flight_path_length = sc.scalar(value=10, unit=sc.units.m)
+
+x_axis = block_rw(float, "x_axis", write_config=BlockWriteConfig(settle_time_s=0.5))
+wavelength_band_0 = sc.array(dims=["tof"], values=[0, 9999999999.0], unit=sc.units.angstrom, dtype="float64")
+wavelength_band_1 = sc.array(dims=["tof"], values=[0.0, 0.07], unit=sc.units.angstrom, dtype="float64")
+
+dae = polarising_dae(det_pixels=[1], frames=50, flipper=flipper, flipper_states=(0.0, 1.0),
+                     intervals=[wavelength_band_0, wavelength_band_1],
+                     total_flight_path_length=total_flight_path_length, monitor=2)
+
+
+def plan() -> Generator[Msg, None, None]:
+  fig, (ax1, ax2) = yield from call_qt_aware(plt.subplots, 2)
+  chained_fit = ChainedLiveFit(method=Linear.fit(), y=[dae.reducer.wavelength_bands[0].calculate_polarisation.name,
+                                                       dae.reducer.wavelength_bands[1].calculate_polarisation.name],
+                               x=bob.name, ax=[ax1, ax2])
+
+  # Subscribe chained_fit to RE and run do a scan for example
+  # chained_fit.get_livefits()[-1].result will give you the fitting results for the last wavelength band
+```
+
+- You are expected to pass in the list of signal names for each independent variable to `y` in order of how you want the subsequent fitting to go.
+- You may also pass in a list of matplotlib axes, which will mean that LiveFitPlots are created per LiveFit, and it will plot the each respective fit to an axis. LiveFitPlots are not created if you do not pass `ax`.
+- Similar to the `y` parameter, you may pass signal names which correspond to uncertainty values for each independent variable.
+
+```{hint}
+The method for fitting is the same across all independent variables.
+```
+
+```{note}
+Parameter uncertainties are not carried over between fits 
+```
+
+```{important}
+If a fit fails to converge, subsequent fits will use their default guess functions
+```

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",