Merge pull request #218 from ISISComputingGroup/centre_of_mass_callback

Centre of mass callback

Author: rerpha

doc/architectural_decisions/008-centre-of-mass.md

@@ -0,0 +1,21 @@
+# Using our own Centre of Mass Callback
+
+## Status
+
+Current
+
+## Context
+
+A decision needs to be made about whether to make changes to upstream Bluesky so that their `CoM` callback works for us, or we make our own.
+
+## Decision
+
+We will be making our own `CoM` callback.
+
+## Justification & Consequences
+
+We attempted to make changes to upstream Bluesky which were rejected, as it adds limits to the functionality of the callback. We also found other limitations with using their callback, such as not being able to have disordered and non-continuous data sent to it without it skewing the calculated value- we need it to work with disordered and non-continuous data as we need to be able to run continuous scans.
+
+This will mean that...
+- Our version of the callback will not be supported by Bluesky and may need changes as Bluesky updates.
+- We can have a version of the callback that is made bespokely for our use cases.

doc/fitting/fitting.md

@@ -214,6 +214,40 @@ 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.
 
+# Centre of Mass
+
+[`CentreOfMass`](ibex_bluesky_core.callbacks.CentreOfMass) is a callback that provides functionality for calculating our definition of Centre of Mass. We calculate centre of mass from the 2D region bounded by min(y), min(x), max(x), and straight-line segments joining (x, y) data points with their nearest neighbours along the x axis.
+
+[`CentreOfMass`](ibex_bluesky_core.callbacks.CentreOfMass) has a property, `result` which stores the centre of mass value once the callback has finished.
+
+In order to use the callback, import `CentreOfMass` from `ibex_bluesky_core.callbacks`.
+```py
+from ibex_bluesky_core.callbacks import CentreOfMass
+```
+
+## Our CoM Algorithm
+
+Given non-continuous arrays of collected data `x` and `y`, ({py:obj}`ibex_bluesky_core.callbacks.CentreOfMass`) returns the `x` value of the centre of mass.
+
+Our use cases require that our algorithm abides to the following rules:
+- Any background on data does not skew the centre of mass
+- The order in which data is received does not skew the centre of mass
+- Should support non-constant point spacing without skewing the centre of mass
+
+*Note that this is designed for only **positive** peaks.*
+
+### Step-by-step
+
+1) Sort `x` and `y` arrays in respect of `x` ascending. This is so that data can be received in any order.
+2) From each `y` element, subtract `min(y)`. This means that any constant background over data is ignored. (Does not work for negative peaks)
+3) Calculate weight/widths for each point; based on it's `x` distances from neighbouring points. This ensures non-constant point spacing is accounted for in our calculation.
+4) For each decomposed shape that makes up the total area under the curve, `CoM` is calculated as the following:
+```{math}
+com_x = \frac{\sum_{}^{}x * y * \text{weight}}{\sum_{}^{}y * \text{weight}}
+```
+
+[`CentreOfMass`](ibex_bluesky_core.callbacks.CentreOfMass) can be used from our callbacks collection. See [ISISCallbacks](ibex_bluesky_core.callbacks.ISISCallbacks).
+
 ## 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.

src/ibex_bluesky_core/callbacks/__init__.py

@@ -16,11 +16,18 @@
 from event_model import RunStart
 from matplotlib.axes import Axes
 
+from ibex_bluesky_core.callbacks._centre_of_mass import (
+    CentreOfMass,
+)
 from ibex_bluesky_core.callbacks._document_logger import DocLoggingCallback
 from ibex_bluesky_core.callbacks._file_logger import (
     HumanReadableFileCallback,
 )
-from ibex_bluesky_core.callbacks._fitting import ChainedLiveFit, 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 +39,7 @@
 
 
 __all__ = [
+    "CentreOfMass",
     "ChainedLiveFit",
     "DocLoggingCallback",
     "HumanReadableFileCallback",
@@ -49,7 +57,7 @@
 class ISISCallbacks:
     """ISIS standard callbacks for use within plans."""
 
-    def __init__(  # noqa: PLR0912
+    def __init__(  # noqa: PLR0912, PLR0915
         self,
         *,
         x: str,
@@ -66,6 +74,7 @@ def __init__(  # noqa: PLR0912
         fit: FitMethod | None = None,
         show_fit_on_plot: bool = True,
         add_peak_stats: bool = True,
+        add_centre_of_mass: bool = True,
         add_live_fit_logger: bool = True,
         live_fit_logger_output_dir: str | PathLike[str] | None = None,
         live_fit_logger_postfix: str = "",
@@ -129,6 +138,7 @@ def _inner():
             fit: The fit method to use when fitting.
             show_fit_on_plot: whether to show fit on plot.
             add_peak_stats: whether to add a peak stats callback.
+            add_centre_of_mass: whether to add a centre of mass callback.
             add_live_fit_logger: whether to add a live fit logger.
             live_fit_logger_output_dir: the output directory for live fit logger.
             live_fit_logger_postfix: the postfix to add to live fit logger.
@@ -142,6 +152,7 @@ def _inner():
         fig = None
         self._subs = []
         self._peak_stats = None
+        self._com = None
         self._live_fit = None
         if measured_fields is None:
             measured_fields = []
@@ -177,6 +188,10 @@ def _inner():
             self._peak_stats = PeakStats(x=x, y=y)
             self._subs.append(self._peak_stats)
 
+        if add_centre_of_mass:
+            self._com = CentreOfMass(x=x, y=y)
+            self._subs.append(self._com)
+
         if (add_plot_cb or show_fit_on_plot) and not ax:
             logger.debug("No axis provided, creating a new one")
             fig, ax = None, None
@@ -260,11 +275,18 @@ def live_fit(self) -> LiveFit:
 
     @property
     def peak_stats(self) -> PeakStats:
-        """The peak stats object containing statistics ie. centre of mass."""
+        """The peak stats object containing statistics ie. bluesky's centre of mass."""
         if self._peak_stats is None:
             raise ValueError("peak stats was not added as a callback.")
         return self._peak_stats
 
+    @property
+    def com(self) -> CentreOfMass:
+        """The centre of mass object containing ibex_bluesky_core's centre of mass."""
+        if self._com is None:
+            raise ValueError("centre of mass was not added as a callback.")
+        return self._com
+
     @property
     def subs(self) -> list[CallbackBase]:
         """The list of subscribed callbacks."""

src/ibex_bluesky_core/callbacks/_centre_of_mass.py

@@ -0,0 +1,59 @@
+import logging
+
+import numpy as np
+from bluesky.callbacks import CollectThenCompute
+
+from ibex_bluesky_core.utils import center_of_mass_of_area_under_curve
+
+logger = logging.getLogger(__name__)
+
+__all__ = ["CentreOfMass"]
+
+
+class CentreOfMass(CollectThenCompute):
+    """Compute centre of mass after a run finishes.
+
+    Calculates the CoM of the 2D region bounded by min(y), min(x), max(x),
+    and straight-line segments joining (x, y) data points with their nearest
+    neighbours along the x axis.
+    """
+
+    def __init__(self, x: str, y: str) -> None:
+        """Initialise the callback.
+
+        Args:
+            x: Name of independent variable in event data
+            y: Name of dependent variable in event data
+
+        """
+        super().__init__()
+        self.x: str = x
+        self.y: str = y
+        self._result: float | None = None
+
+    @property
+    def result(self) -> float | None:
+        """The centre-of-mass calculated by this callback."""
+        return self._result
+
+    def compute(self) -> None:
+        """Calculate statistics at the end of the run."""
+        x_values = []
+        y_values = []
+
+        for event in self._events:
+            if self.x not in event["data"]:
+                raise ValueError(f"{self.x} is not in event document.")
+
+            if self.y not in event["data"]:
+                raise ValueError(f"{self.y} is not in event document.")
+
+            x_values.append(event["data"][self.x])
+            y_values.append(event["data"][self.y])
+
+        if not x_values:
+            return
+
+        x_data = np.array(x_values, dtype=np.float64)
+        y_data = np.array(y_values, dtype=np.float64)
+        (self._result, _) = center_of_mass_of_area_under_curve(x_data, y_data)

src/ibex_bluesky_core/fitting.py

@@ -26,6 +26,8 @@
     "Trapezoid",
 ]
 
+from ibex_bluesky_core.utils import center_of_mass_of_area_under_curve
+
 
 class FitMethod:
     """Tell LiveFit how to fit to a scan. Has a Model function and a Guess function.
@@ -545,49 +547,11 @@ def guess(
         return guess
 
 
-def _center_of_mass_and_mass(
-    x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
-) -> tuple[float, float]:
-    """Compute the centre of mass of the area under a curve defined by a series of (x, y) points.
-
-    The "area under the curve" is a shape bounded by:
-    - min(y), along the bottom edge
-    - min(x), on the left-hand edge
-    - max(x), on the right-hand edge
-    - straight lines joining (x, y) data points to their nearest neighbours
-        along the x-axis, along the top edge
-    This is implemented by geometric decomposition of the shape into a series of trapezoids,
-    which are further decomposed into rectangular and triangular regions.
-    """
-    sort_indices = np.argsort(x, kind="stable")
-    x = np.take_along_axis(x, sort_indices, axis=None)
-    y = np.take_along_axis(y - np.min(y), sort_indices, axis=None)
-    widths = np.diff(x)
-
-    # Area under the curve for two adjacent points is a right trapezoid.
-    # Split that trapezoid into a rectangular region, plus a right triangle.
-    # Find area and effective X CoM for each.
-    rect_areas = widths * np.minimum(y[:-1], y[1:])
-    rect_x_com = (x[:-1] + x[1:]) / 2.0
-    triangle_areas = widths * np.abs(y[:-1] - y[1:]) / 2.0
-    triangle_x_com = np.where(
-        y[:-1] > y[1:], x[:-1] + (widths / 3.0), x[:-1] + (2.0 * widths / 3.0)
-    )
-
-    total_area = np.sum(rect_areas + triangle_areas)
-    if total_area == 0.0:
-        # If all data was flat, return central x
-        return (x[0] + x[-1]) / 2.0, 0
-
-    com = np.sum(rect_areas * rect_x_com + triangle_areas * triangle_x_com) / total_area
-    return com, total_area
-
-
 def _guess_cen_and_width(
     x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
 ) -> tuple[float, float]:
     """Guess the center and width of a positive peak."""
-    com, total_area = _center_of_mass_and_mass(x, y)
+    com, total_area = center_of_mass_of_area_under_curve(x, y)
     y_range = np.max(y) - np.min(y)
     if y_range == 0.0:
         width = (np.max(x) - np.min(x)) / 2

src/ibex_bluesky_core/utils.py

@@ -6,12 +6,15 @@
 from typing import Any, Protocol
 
 import matplotlib
+import numpy as np
+import numpy.typing as npt
 import scipp as sc
 from bluesky.protocols import NamedMovable, Readable
 
 __all__ = [
     "NamedReadableAndMovable",
     "calculate_polarisation",
+    "center_of_mass_of_area_under_curve",
     "centred_pixel",
     "get_pv_prefix",
     "is_matplotlib_backend_qt",
@@ -52,6 +55,67 @@ class NamedReadableAndMovable(Readable[Any], NamedMovable[Any], Protocol):
     """Abstract class for type checking that an object is readable, named and movable."""
 
 
+def center_of_mass_of_area_under_curve(
+    x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
+) -> tuple[float, float]:
+    """Compute the centre of mass of the area under a curve defined by a series of (x, y) points.
+
+    The "area under the curve" is a shape bounded by:
+    - min(y), along the bottom edge
+    - min(x), on the left-hand edge
+    - max(x), on the right-hand edge
+    - straight lines joining (x, y) data points to their nearest neighbours
+    along the x-axis, along the top edge
+
+    This is implemented by geometric decomposition of the shape into a series of trapezoids,
+    which are further decomposed into rectangular and triangular regions.
+
+    Returns a tuple of the centre of mass and the total area under the curve.
+
+    """
+    # Sorting here avoids special-cases with disordered points, which may occur
+    # from a there-and-back scan, or from an adaptive scan.
+    sort_indices = np.argsort(x, kind="stable")
+    x = np.take_along_axis(x, sort_indices, axis=None)
+    y = np.take_along_axis(y - np.min(y), sort_indices, axis=None)
+
+    # If the data points are "fence-posts", this calculates the x width of
+    # each "fence panel".
+    widths = np.diff(x)
+
+    # Area under the curve for two adjacent points is a right trapezoid.
+    # Split that trapezoid into a rectangular region, plus a right triangle.
+    # Find area and effective X CoM for each.
+
+    # We want the area of the rectangular part of the right trapezoid.
+    # This is width * [height of either left or right point, whichever is lowest]
+    rect_areas = widths * np.minimum(y[:-1], y[1:])
+    # CoM of a rectangle in x is simply the average x.
+    rect_x_com = (x[:-1] + x[1:]) / 2.0
+
+    # Now the area of the triangular part of the right trapezoid - this is
+    # width * height / 2, where height is the absolute difference between the
+    # two y values.
+    triangle_areas = widths * np.abs(y[:-1] - y[1:]) / 2.0
+    # CoM of a right triangle is 1/3 along the base, from the right angle
+    # y[:-1] > y[1:] is true if y_[n] > y_[n+1], (i.e. if the right angle is on the
+    # left-hand side of the triangle).
+    # If that's true, the CoM lies 1/3 of the way along the x axis
+    # Otherwise, the CoM lies 2/3 of the way along the x axis (1/3 from the right angle)
+    triangle_x_com = np.where(
+        y[:-1] > y[1:], x[:-1] + (widths / 3.0), x[:-1] + (2.0 * widths / 3.0)
+    )
+
+    total_area = np.sum(rect_areas + triangle_areas)
+    if total_area == 0.0:
+        # If all data was flat, return central x
+        return (x[0] + x[-1]) / 2.0, total_area
+
+    return np.sum(
+        rect_areas * rect_x_com + triangle_areas * triangle_x_com
+    ) / total_area, total_area
+
+
 def calculate_polarisation(
     a: sc.Variable | sc.DataArray, b: sc.Variable | sc.DataArray
 ) -> sc.Variable | sc.DataArray: