Merge pull request #186 from ISISComputingGroup/negative_trapezoids

Negative trapezoids

Author: jackbdoughty

doc/fitting/standard_fits.md

@@ -160,4 +160,26 @@ g(x) = \max(f(x), \text{background})
 y = \min(g(x), \text{background} + \text{height})
 ```
 
-![TrapezoidModel](../_static/images_fits/trapezoid.png)
+![TrapezoidModel](../_static/images_fits/trapezoid.png)
+
+## Negative Trapezoid
+
+API Reference: [`NegativeTrapezoid`](ibex_bluesky_core.fitting.NegativeTrapezoid)
+
+This model is the same shape as the trapezoid described above, but with a negative height.s
+
+- `cen` - The centre (x) of the model
+- `gradient` - How steep the edges of the trapezoid are
+- `height` - The maximum height of the model below `background`
+- `background` - The maximum value (y) of the model
+- `y_offset` - Acts as a width factor for the trapezoid. If you extrapolate the sides of the trapezoid until they meet, this value represents the y coord of this point minus height and background.
+
+```{math}
+f(x) = \text{y_offset} - \text{height} + \text{background} + \text{gradient} * |x - \text{cen}|
+```
+```{math}
+g(x) = \max(f(x), \text{background} - \text{height})
+```
+```{math}
+y = \min(g(x), \text{background})
+```

src/ibex_bluesky_core/fitting.py

@@ -20,6 +20,7 @@
     "Gaussian",
     "Linear",
     "Lorentzian",
+    "NegativeTrapezoid",
     "Polynomial",
     "SlitScan",
     "TopHat",
@@ -475,6 +476,57 @@ 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)
+    y_range = np.max(y) - np.min(y)
+    if y_range == 0.0:
+        width = (np.max(x) - np.min(x)) / 2
+    else:
+        width = total_area / y_range
+    return com, width
+
+
 class TopHat(Fit):
     """Top Hat Fitting."""
 
@@ -502,26 +554,31 @@ def guess(
         def guess(
             x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
         ) -> dict[str, lmfit.Parameter]:
-            top = np.where(y > np.mean(y))[0]
-            # Guess that any value above the mean is the top part
-
-            if len(top) > 0:
-                width = x[np.max(top)] - x[np.min(top)]
-            else:
-                width = (max(x) - min(x)) / 2
+            cen, width = _guess_cen_and_width(x, y)
 
             init_guess = {
-                "cen": lmfit.Parameter("cen", np.mean(x)),
-                "width": lmfit.Parameter("width", width),
-                "height": lmfit.Parameter("height", (max(y) - min(y))),
-                "background": lmfit.Parameter("background", min(y)),
+                "cen": lmfit.Parameter("cen", cen),
+                "width": lmfit.Parameter("width", width, min=0),
+                "height": lmfit.Parameter(
+                    "height",
+                    np.max(y) - np.min(y),
+                ),
+                "background": lmfit.Parameter("background", np.min(y)),
             }
 
             return init_guess
 
         return guess
 
 
+def _guess_trapezoid_gradient(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> float:
+    gradients = np.zeros_like(x[1:], dtype=np.float64)
+    x_diffs = x[:-1] - x[1:]
+    y_diffs = y[:-1] - y[1:]
+    np.divide(y_diffs, x_diffs, out=gradients, where=x_diffs != 0)
+    return np.max(np.abs(gradients))
+
+
 class Trapezoid(Fit):
     """Trapezoid Fitting."""
 
@@ -557,37 +614,74 @@ def guess(
         def guess(
             x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
         ) -> dict[str, lmfit.Parameter]:
-            top = np.where(y > np.mean(y))[0]
-            # Guess that any value above the y mean is the top part
+            cen, width = _guess_cen_and_width(x, y)
+            gradient_guess = _guess_trapezoid_gradient(x, y)
 
-            cen = np.mean(x)
+            height = np.max(y) - np.min(y)
             background = np.min(y)
-            height = np.max(y) - background
+            y_offset = gradient_guess * width / 2.0
 
-            if top.size > 0:
-                i = np.min(top)
-                x1 = x[i]  # x1 is the left of the top part
-            else:
-                width_top = (np.max(x) - np.min(x)) / 2
-                x1 = cen - width_top / 2
+            init_guess = {
+                "cen": lmfit.Parameter("cen", cen, min=np.min(x), max=np.max(x)),
+                "gradient": lmfit.Parameter("gradient", gradient_guess, min=0),
+                "height": lmfit.Parameter("height", height, min=0),
+                "background": lmfit.Parameter("background", background),
+                "y_offset": lmfit.Parameter("y_offset", y_offset),
+            }
 
-            x0 = 0.5 * (np.min(x) + x1)  # Guess that x0 is half way between min(x) and x1
+            return init_guess
 
-            if height == 0.0:
-                gradient = 0.0
-            else:
-                gradient = height / (x1 - x0)
+        return guess
+
+
+class NegativeTrapezoid(Fit):
+    """Negative Trapezoid Fitting."""
+
+    equation = """
+    y = clip(y_offset - height + background + gradient * abs(x - cen),
+     background - height, background)"""
+
+    @classmethod
+    def model(cls, *args: int) -> lmfit.Model:
+        """Negative Trapezoid Model."""
 
-            y_intercept0 = np.max(y) - gradient * x1  # To find the slope function
-            y_tip = gradient * cen + y_intercept0
-            y_offset = y_tip - height - background
+        def model(
+            x: npt.NDArray[np.float64],
+            cen: float,
+            gradient: float,
+            height: float,
+            background: float,
+            y_offset: float,  # Acts as a width multiplier
+        ) -> npt.NDArray[np.float64]:
+            y = y_offset - height + background + gradient * np.abs(x - cen)
+            y = np.maximum(y, background - height)
+            y = np.minimum(y, background)
+            return y
+
+        return lmfit.Model(model, name=f"{cls.__name__}  [{cls.equation}]")
+
+    @classmethod
+    def guess(
+        cls, *args: int
+    ) -> Callable[[npt.NDArray[np.float64], npt.NDArray[np.float64]], dict[str, lmfit.Parameter]]:
+        """Negative Trapezoid Guessing."""
+
+        def guess(
+            x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
+        ) -> dict[str, lmfit.Parameter]:
+            cen, width = _guess_cen_and_width(x, -y)
+            gradient_guess = _guess_trapezoid_gradient(x, y)
+
+            height = np.max(y) - np.min(y)
+            background = np.max(y)
+            y_offset = -gradient_guess * width / 2.0
 
             init_guess = {
-                "cen": lmfit.Parameter("cen", cen),
-                "gradient": lmfit.Parameter("gradient", gradient, min=0),
+                "cen": lmfit.Parameter("cen", cen, min=np.min(x), max=np.max(x)),
+                "gradient": lmfit.Parameter("gradient", gradient_guess, min=0),
                 "height": lmfit.Parameter("height", height, min=0),
                 "background": lmfit.Parameter("background", background),
-                "y_offset": lmfit.Parameter("y_offset", y_offset, min=0),
+                "y_offset": lmfit.Parameter("y_offset", y_offset),
             }
 
             return init_guess

tests/callbacks/fitting/test_fitting_methods.py

@@ -17,6 +17,7 @@
     Gaussian,
     Linear,
     Lorentzian,
+    NegativeTrapezoid,
     Polynomial,
     SlitScan,
     TopHat,
@@ -599,7 +600,7 @@ def test_width_guess(self):
 
             outp = TopHat.guess()(x, y)
 
-            assert outp["width"] == pytest.approx(1.0, rel=1e-2)
+            assert outp["width"] == pytest.approx(2.0, rel=1e-2)
 
         def test_guess_given_flat_data(self):
             x = np.arange(-5.0, 5.0, 1.0, dtype=np.float64)
@@ -688,15 +689,15 @@ def test_gradient_guess(self):
 
             outp = Trapezoid.guess()(x, y)
 
-            assert outp["gradient"] == pytest.approx(10.0, rel=1e-2)
+            assert outp["gradient"] == pytest.approx(8.0, rel=1e-2)
 
         def test_y_offset_guess(self):
-            x = np.arange(-5.0, 5.0, 1.0, dtype=np.float64)
+            x = np.linspace(-5.0, 5.0, num=11, dtype=np.float64)
             y = np.array([1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0, 1.0], dtype=np.float64)
 
             outp = Trapezoid.guess()(x, y)
 
-            x1 = np.arange(-6.0, 6.0, 1.0, dtype=np.float64)
+            x1 = np.linspace(-6.0, 6.0, num=13, dtype=np.float64)
             y1 = np.array(
                 [1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0, 1.0], dtype=np.float64
             )
@@ -713,3 +714,112 @@ def test_guess_given_flat_data(self):
             outp = Trapezoid.guess()(x, y)
             # check that with flat data gradient guess is 0
             assert outp["gradient"] == pytest.approx(0.0, rel=1e-2)
+
+
+class TestNegativeTrapezoid:
+    class TestNegativeTrapezoidModel:
+        def test_negative_trapezoid_model(self):
+            x = np.arange(-5.0, 6.0, 1.0, dtype=np.float64)
+            cen = 0
+            y_offset = -1
+            height = 1
+            background = 1
+            gradient = 1
+
+            outp = NegativeTrapezoid.model().func(
+                x,
+                cen=cen,
+                y_offset=y_offset,
+                height=height,
+                background=background,
+                gradient=gradient,
+            )
+
+            assert background == pytest.approx(np.max(outp), rel=1e-2)
+            assert background - height == pytest.approx(np.min(outp), rel=1e-2)
+
+            outp1 = NegativeTrapezoid.model().func(
+                x,
+                cen=cen + 3,
+                y_offset=y_offset,
+                height=height,
+                background=background,
+                gradient=gradient - 0.5,
+            )
+
+            # check centre moves when data is shifted
+            assert np.mean(x[np.where(outp < background)]) < np.mean(
+                x[np.where(outp1 < background)]
+            )
+
+            # check gradient: a greater gradient means smaller average y values as wider
+            assert np.mean(outp) > np.mean(outp1)
+
+            outp2 = NegativeTrapezoid.model().func(
+                x,
+                cen=cen,
+                y_offset=y_offset - 5,
+                height=height,
+                background=background,
+                gradient=gradient,
+            )
+
+            # check y_offset: a smaller y_offset means smaller average y values as wider
+            assert np.mean(outp) > np.mean(outp2)
+
+    class TestNegativeTrapezoidGuess:
+        def test_background_guess(self):
+            x = np.array([-1.0, 0.0, 1.0, 2.0, 3.0], dtype=np.float64)
+            y = np.array([-1.0, -2.0, -2.0, -2.0, -1.0], dtype=np.float64)
+
+            outp = NegativeTrapezoid.guess()(x, y)
+
+            assert outp["background"] == pytest.approx(-1.0, rel=1e-2)
+
+        def test_cen_height_guess(self):
+            x = np.array([-1.0, 0.0, 1.0, 2.0, 3.0], dtype=np.float64)
+            y = np.array([-1.0, -1.0, -2.0, -1.0, -1.0], dtype=np.float64)
+
+            outp = NegativeTrapezoid.guess()(x, y)
+
+            assert outp["cen"] == pytest.approx(1.0, rel=1e-2)
+            assert outp["height"] == pytest.approx(1.0, rel=1e-2)
+
+        def test_gradient_guess(self):
+            x = np.array([-4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0], dtype=np.float64)
+            y = np.array([1.0, 2.0, 4.0, 8.0, 16.0, 8.0, 4.0, 2.0, 1.0], dtype=np.float64)
+
+            # Should choose x = -1.0 as x1 and x = -2.5 as x0
+            # height = 16 - 1 = 15
+            # gradient = 15 / (-1 - -2.5) = 10
+
+            outp = NegativeTrapezoid.guess()(x, y)
+
+            assert outp["gradient"] == pytest.approx(8.0, rel=1e-2)
+
+        def test_y_offset_guess(self):
+            x = np.linspace(-5.0, 5.0, num=11, dtype=np.float64)
+            y = np.array(
+                [-1.0, -1.0, -1.0, -2.0, -3.0, -3.0, -3.0, -2.0, -1.0, -1.0, -1.0], dtype=np.float64
+            )
+
+            outp = NegativeTrapezoid.guess()(x, y)
+
+            x1 = np.linspace(-6.0, 6.0, num=13, dtype=np.float64)
+            y1 = np.array(
+                [-1.0, -1.0, -1.0, -2.0, -3.0, -3.0, -3.0, -3.0, -3.0, -2.0, -1.0, -1.0, -1.0],
+                dtype=np.float64,
+            )
+
+            outp1 = NegativeTrapezoid.guess()(x1, y1)
+
+            # Assert that with a greater top width, y_offset decreases
+            assert outp["y_offset"] > outp1["y_offset"]
+
+        def test_guess_given_flat_data(self):
+            x = np.arange(-5.0, 5.0, 1.0, dtype=np.float64)
+            y = np.zeros_like(x, dtype=np.float64) + 1
+
+            outp = NegativeTrapezoid.guess()(x, y)
+            # check that with flat data gradient guess is 0
+            assert outp["gradient"] == pytest.approx(0.0, rel=1e-2)