Merge pull request #245 from scipp/toa-tof-rebin

Resample TOA-TOF-converted histogram data using `rebin`

Author: SimonHeybrock

src/ess/reduce/time_of_flight/__init__.py

@@ -13,7 +13,6 @@
     resample_monitor_time_of_flight_data,
 )
 from .simulation import simulate_beamline
-from .to_events import to_events
 from .types import (
     DetectorLtotal,
     DetectorTofData,
@@ -58,5 +57,4 @@
     "resample_detector_time_of_flight_data",
     "resample_monitor_time_of_flight_data",
     "simulate_beamline",
-    "to_events",
 ]

src/ess/reduce/time_of_flight/eto_to_tof.py

@@ -27,7 +27,7 @@
     MonitorType,
     RunType,
 )
-from .to_events import to_events
+from .resample import rebin_strictly_increasing
 from .types import (
     DetectorLtotal,
     DetectorTofData,
@@ -664,51 +664,13 @@ def monitor_time_of_flight_data(
     )
 
 
-def _resample_tof_data(da: sc.DataArray) -> sc.DataArray:
-    """
-    Histogrammed data that has been converted to `tof` will typically have
-    unsorted bin edges (due to either wrapping of `time_of_flight` or wavelength
-    overlap between subframes).
-    This function re-histograms the data to ensure that the bin edges are sorted.
-    It makes use of the ``to_events`` helper which generates a number of events in each
-    bin with a uniform distribution. The new events are then histogrammed using a set of
-    sorted bin edges.
-
-    WARNING:
-    This function is highly experimental, has limitations and should be used with
-    caution. It is a workaround to the issue that rebinning data with unsorted bin
-    edges is not supported in scipp.
-    As such, this function is not part of the default set of providers, and needs to be
-    inserted manually into the workflow.
-
-    Parameters
-    ----------
-    da:
-        Histogrammed data with the time-of-flight coordinate.
-    """
-    dim = next(iter(set(da.dims) & {"time_of_flight", "tof"}))
-    data = da.rename_dims({dim: "tof"}).drop_coords(
-        [name for name in da.coords if name != "tof"]
-    )
-    events = to_events(data, "event")
-
-    # Define a new bin width, close to the original bin width.
-    # TODO: this could be a workflow parameter
-    coord = da.coords["tof"]
-    bin_width = (coord[dim, 1:] - coord[dim, :-1]).nanmedian()
-    rehist = events.hist(tof=bin_width)
-    return rehist.assign_coords(
-        {key: var for key, var in da.coords.items() if dim not in var.dims}
-    )
-
-
 def resample_detector_time_of_flight_data(
     da: DetectorTofData[RunType],
 ) -> ResampledDetectorTofData[RunType]:
     """
     Resample the detector time-of-flight data to ensure that the bin edges are sorted.
     """
-    return ResampledDetectorTofData(_resample_tof_data(da))
+    return ResampledDetectorTofData[RunType](rebin_strictly_increasing(da, dim='tof'))
 
 
 def resample_monitor_time_of_flight_data(
@@ -717,7 +679,9 @@ def resample_monitor_time_of_flight_data(
     """
     Resample the monitor time-of-flight data to ensure that the bin edges are sorted.
     """
-    return ResampledMonitorTofData(_resample_tof_data(da))
+    return ResampledMonitorTofData[RunType, MonitorType](
+        rebin_strictly_increasing(da, dim='tof')
+    )
 
 
 def default_parameters() -> dict:

src/ess/reduce/time_of_flight/resample.py

@@ -0,0 +1,97 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
+
+
+import numpy as np
+import scipp as sc
+
+
+def find_strictly_increasing_sections(var: sc.Variable) -> list[slice]:
+    """
+    Find strictly increasing sections in a coordinate dimension (minimum length 2).
+
+    Parameters
+    ----------
+    var:
+        The variable to analyze, which should be one-dimensional.
+
+    Returns
+    -------
+    sections:
+        Slice objects that can be used extract strictly increasing sections.
+    """
+    values = var.values
+    finite = np.isfinite(values)
+    increasing = (np.sign(np.diff(values)) > 0) & finite[:-1] & finite[1:]
+    # 1 marks the start of an increasing section, -1 marks the end
+    transitions = np.diff(np.concatenate(([False], increasing, [False])).astype(int))
+    section_starts = np.where(transitions == 1)[0]
+    section_ends = np.where(transitions == -1)[0] + np.array(1)
+    return [
+        slice(start, end)
+        for start, end in zip(section_starts, section_ends, strict=True)
+        if end - start >= 2  # Ensure section has at least 2 points
+    ]
+
+
+def get_min_max(
+    var: sc.Variable, *, dim: str, slices: list[slice]
+) -> tuple[sc.Variable, sc.Variable]:
+    if not slices:
+        raise ValueError("No strictly increasing sections found.")
+    combined = sc.concat([var[dim, slice] for slice in slices], dim)
+    return combined.min(), combined.max()
+
+
+def make_regular_grid(
+    var: sc.Variable, *, dim: str, slices: list[slice]
+) -> sc.Variable:
+    """
+    Create a regular grid variable based on the min and max of the slices.
+
+    The grid is constructed such that it includes the minimum and maximum values
+    of the strictly increasing sections, with a step size equal to the difference
+    between the first two values of the section with the minimum start value (which is
+    not necessarily the first section).
+    """
+    min_val, max_val = get_min_max(var, dim=dim, slices=slices)
+    first: sc.Variable | None = None
+    for s in slices:
+        first = var[dim, s]
+        if sc.identical(first[0], min_val):
+            break
+    if first is None:
+        # This should not happen if slices are correctly identified and passed from
+        # find_strictly_increasing_sections.
+        raise ValueError("Section is not strictly increasing.")
+    step = first[1] - first[0]
+    return sc.arange(
+        dim=dim,
+        start=min_val.value,
+        stop=max_val.value + step.value,  # Ensure the last bin edge is included
+        step=step.value,
+        unit=step.unit,
+        dtype=step.dtype,
+    )
+
+
+def rebin_strictly_increasing(da: sc.DataArray, dim: str) -> sc.DataArray:
+    """
+    Find strictly monotonic sections in a coordinate dimension and rebin the data array
+    into a regular grid based on these sections.
+    """
+    # Ensure the dimension is named like the coordinate.
+    da = da.rename_dims({da.coords[dim].dim: dim})
+    slices = find_strictly_increasing_sections(da.coords[dim])
+    if len(slices) == 1:
+        return da[dim, slices[0]]
+    if not slices:
+        raise ValueError("No strictly increasing sections found.")
+    if da.coords[dim].dtype not in (sc.DType.float64, sc.DType.float32):
+        # rebin does not like integer coords.
+        da = da.assign_coords({dim: da.coords[dim].to(dtype='float64')})
+    # Slices refer to the indices in the coord, which are bin edges. For slicing data
+    # we need to stop at the last index minus one.
+    sections = [da[dim, section.start : section.stop - 1] for section in slices]
+    edges = make_regular_grid(da.coords[dim], dim=dim, slices=slices)
+    return sc.reduce([sc.rebin(section, {dim: edges}) for section in sections]).sum()