Fix threading

packages/essreduce/src/ess/reduce/unwrap/__init__.py

@@ -28,6 +28,7 @@
 from .types import (
     DetectorLtotal,
     ErrorLimitedLookupTable,
+    FrameUnwrapBackend,
     LookupTable,
     LookupTableFilename,
     LookupTableRelativeErrorThreshold,
@@ -48,6 +49,7 @@
     "DiskChoppers",
     "DistanceResolution",
     "ErrorLimitedLookupTable",
+    "FrameUnwrapBackend",
     "GenericUnwrapWorkflow",
     "LookupTable",
     "LookupTableFilename",

packages/essreduce/src/ess/reduce/unwrap/to_wavelength.py

@@ -7,6 +7,7 @@
 event_time_offset coordinates to data with a time-of-flight coordinate.
 """
 
+import warnings
 from collections.abc import Callable
 from dataclasses import asdict
 
@@ -16,11 +17,6 @@
 import scippnexus as snx
 from scippneutron._utils import elem_unit
 
-try:
-    from .interpolator_numba import Interpolator as InterpolatorImpl
-except ImportError:
-    from .interpolator_scipy import Interpolator as InterpolatorImpl
-
 from ..nexus.types import (
     Component,
     EmptyDetector,
@@ -37,6 +33,7 @@
 from .types import (
     DetectorLtotal,
     ErrorLimitedLookupTable,
+    FrameUnwrapBackend,
     KeepEventTimeOffset,
     LookupTable,
     LookupTableRelativeErrorThreshold,
@@ -47,13 +44,47 @@
 )
 
 
+def _get_interpolator_class(backend: FrameUnwrapBackend) -> type:
+    if backend == FrameUnwrapBackend.scipy:
+        from .interpolator_scipy import Interpolator
+
+        return Interpolator
+
+    try:
+        from numba import get_num_threads, threading_layer
+
+        # Initialize the threading layer now so we can inspect the selected backend
+        # and fall back before interpolation.
+        get_num_threads()
+        if threading_layer() in {'omp', 'tbb'}:
+            from .interpolator_numba import Interpolator
+
+            return Interpolator
+    except (ImportError, ValueError):
+        pass
+
+    warnings.warn(
+        "The 'numba' frame-unwrapping backend was requested, but Numba is "
+        "unavailable or did not select a thread-safe threading layer. Falling "
+        "back to the 'scipy' backend. This fallback is deprecated and will be "
+        "an error in a future release. Select the SciPy backend explicitly with "
+        "backend='scipy'.",
+        FutureWarning,
+        stacklevel=3,
+    )
+    from .interpolator_scipy import Interpolator
+
+    return Interpolator
+
+
 class WavelengthInterpolator:
     def __init__(
         self,
         lookup: sc.DataArray,
         distance_unit: str,
         time_unit: str,
         wavelength_unit: str = 'angstrom',
+        backend: FrameUnwrapBackend = FrameUnwrapBackend.numba,
     ):
         """
         Interpolator object that converts event_time_offset and distances to
@@ -96,7 +127,7 @@ def __init__(
 
         distances = lookup.coords["distance"].to(unit=distance_unit, copy=False)
 
-        self._interpolator = InterpolatorImpl(
+        self._interpolator = _get_interpolator_class(backend)(
             time_edges=time_coord,
             distance_edges=distances.values,
             values=(
@@ -137,7 +168,10 @@ def __call__(
 
 
 def _compute_wavelength_histogram(
-    da: sc.DataArray, lookup: LookupTable, ltotal: sc.Variable
+    da: sc.DataArray,
+    lookup: LookupTable,
+    ltotal: sc.Variable,
+    backend: FrameUnwrapBackend,
 ) -> sc.DataArray:
     # In NeXus, 'time_of_flight' is the canonical name in NXmonitor, but in some files,
     # it may be called 'tof' or 'frame_time'.
@@ -164,7 +198,10 @@ def _compute_wavelength_histogram(
 
     # Create linear interpolator
     interp = WavelengthInterpolator(
-        lookup.array, distance_unit=ltotal.unit, time_unit=eto_unit
+        lookup.array,
+        distance_unit=ltotal.unit,
+        time_unit=eto_unit,
+        backend=backend,
     )
 
     # Compute wavelengths of the bin edges using the interpolator
@@ -250,6 +287,7 @@ def _prepare_wavelength_interpolation_inputs(
     lookup: LookupTable,
     ltotal: sc.Variable,
     pulse_stride_offset: int | None,
+    backend: FrameUnwrapBackend,
 ) -> dict:
     """
     Prepare the inputs required for the wavelength interpolation.
@@ -269,13 +307,18 @@ def _prepare_wavelength_interpolation_inputs(
         When pulse-skipping, the offset of the first pulse in the stride. This is
         typically zero but can be a small integer < pulse_stride.
         If None, a guess is made.
+    backend:
+        Backend used to interpolate the wavelength lookup table.
     """
     etos = da.bins.coords["event_time_offset"].to(dtype=float, copy=False)
     eto_unit = elem_unit(etos)
 
     # Create linear interpolator
     interp = WavelengthInterpolator(
-        lookup.array, distance_unit=ltotal.unit, time_unit=eto_unit
+        lookup.array,
+        distance_unit=ltotal.unit,
+        time_unit=eto_unit,
+        backend=backend,
     )
 
     # Operate on events (broadcast distances to all events)
@@ -343,13 +386,15 @@ def _compute_wavelength_events(
     lookup: LookupTable,
     ltotal: sc.Variable,
     pulse_stride_offset: int | None,
+    backend: FrameUnwrapBackend,
     keep_event_time_offset: bool,
 ) -> sc.DataArray:
     inputs = _prepare_wavelength_interpolation_inputs(
         da=da,
         lookup=lookup,
         ltotal=ltotal,
         pulse_stride_offset=pulse_stride_offset,
+        backend=backend,
     )
 
     # Compute wavelength for all neutrons using the interpolator
@@ -476,17 +521,21 @@ def _compute_wavelength_data(
     lookup: ErrorLimitedLookupTable[RunType, Component],
     ltotal: sc.Variable,
     pulse_stride_offset: int,
+    backend: FrameUnwrapBackend,
     keep_event_time_offset: bool,
 ) -> sc.DataArray:
     if da.bins is None:
-        data = _compute_wavelength_histogram(da=da, lookup=lookup, ltotal=ltotal)
+        data = _compute_wavelength_histogram(
+            da=da, lookup=lookup, ltotal=ltotal, backend=backend
+        )
         out = rebin_strictly_increasing(data, dim='wavelength')
     else:
         out = _compute_wavelength_events(
             da=da,
             lookup=lookup,
             ltotal=ltotal,
             pulse_stride_offset=pulse_stride_offset,
+            backend=backend,
             keep_event_time_offset=keep_event_time_offset,
         )
     return out.assign_coords(Ltotal=ltotal)
@@ -498,6 +547,7 @@ def detector_wavelength_data(
     ltotal: DetectorLtotal[RunType],
     pulse_stride_offset: PulseStrideOffset,
     keep_event_time_offset: KeepEventTimeOffset,
+    backend: FrameUnwrapBackend = FrameUnwrapBackend.numba,
 ) -> WavelengthDetector[RunType]:
     """
     Convert the time-of-arrival (event_time_offset) data to wavelength data using a
@@ -520,13 +570,16 @@ def detector_wavelength_data(
     keep_event_time_offset:
         Whether to keep the event_time_offset coordinate after converting to
         wavelength.
+    backend:
+        Backend used to interpolate the wavelength lookup table.
     """
     return WavelengthDetector[RunType](
         _compute_wavelength_data(
             da=detector_data,
             lookup=lookup,
             ltotal=ltotal,
             pulse_stride_offset=pulse_stride_offset,
+            backend=backend,
             keep_event_time_offset=keep_event_time_offset,
         )
     )
@@ -538,6 +591,7 @@ def monitor_wavelength_data(
     ltotal: MonitorLtotal[RunType, MonitorType],
     pulse_stride_offset: PulseStrideOffset,
     keep_event_time_offset: KeepEventTimeOffset,
+    backend: FrameUnwrapBackend = FrameUnwrapBackend.numba,
 ) -> WavelengthMonitor[RunType, MonitorType]:
     """
     Convert the time-of-arrival (event_time_offset) data to wavelength data using a
@@ -560,13 +614,16 @@ def monitor_wavelength_data(
     keep_event_time_offset:
         Whether to keep the event_time_offset coordinate after converting to
         wavelength.
+    backend:
+        Backend used to interpolate the wavelength lookup table.
     """
     return WavelengthMonitor[RunType, MonitorType](
         _compute_wavelength_data(
             da=monitor_data,
             lookup=lookup,
             ltotal=ltotal,
             pulse_stride_offset=pulse_stride_offset,
+            backend=backend,
             keep_event_time_offset=keep_event_time_offset,
         )
     )

packages/essreduce/src/ess/reduce/unwrap/types.py

@@ -18,6 +18,13 @@ class WavelengthLutMode(StrEnum):
     file = 'file'
 
 
+class FrameUnwrapBackend(StrEnum):
+    """Backend used to interpolate the frame-unwrapping lookup table."""
+
+    numba = 'numba'
+    scipy = 'scipy'
+
+
 class LookupTableFilename(sl.Scope[RunType, Component, str], str):
     """Filename of the wavelength lookup table."""
 

packages/essreduce/src/ess/reduce/unwrap/workflow.py

@@ -6,7 +6,7 @@
 
 from ..nexus import GenericNeXusWorkflow
 from ..nexus.types import AnyRun
-from . import WavelengthLutMode, lut, to_wavelength
+from . import FrameUnwrapBackend, WavelengthLutMode, lut, to_wavelength
 
 
 def GenericUnwrapWorkflow(
@@ -47,7 +47,6 @@ def GenericUnwrapWorkflow(
         Mode for creating the wavelength lookup table. Possible values are
         'analytical', 'simulation', and 'file'. See
         https://scipp.github.io/ess/reduce/user-guide/unwrap/lut-building-methods.html
-
     Returns
     -------
     :
@@ -62,6 +61,7 @@ def GenericUnwrapWorkflow(
         wf.insert(provider)
     for key, value in lut.default_parameters(wavelength_from=wavelength_from).items():
         wf[key] = value
+    wf[FrameUnwrapBackend] = FrameUnwrapBackend.numba
 
     return wf
 

packages/essreduce/tests/unwrap/interpolator_test.py

@@ -1,8 +1,13 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
 
+from unittest.mock import patch
+
+import numba
 import numpy as np
+import pytest
 
+from ess.reduce.unwrap import FrameUnwrapBackend, to_wavelength
 from ess.reduce.unwrap.interpolator_numba import (
     Interpolator as InterpolatorNumba,
 )
@@ -100,3 +105,55 @@ def test_numba_and_scipy_interpolators_yield_same_results_with_values_on_edges()
     numba_result = numba_interp(times, distances)
     scipy_result = scipy_interp(times, distances)
     assert np.allclose(numba_result, scipy_result, equal_nan=True)
+
+
+def test_scipy_backend_selects_scipy(monkeypatch):
+    def fail_if_called():
+        raise AssertionError('Numba should not be checked for the SciPy backend.')
+
+    monkeypatch.setattr(numba, 'get_num_threads', fail_if_called)
+
+    impl = to_wavelength._get_interpolator_class(FrameUnwrapBackend.scipy)
+
+    assert impl is InterpolatorScipy
+
+
+def test_numba_backend_selects_numba():
+    impl = to_wavelength._get_interpolator_class(FrameUnwrapBackend.numba)
+
+    assert impl is InterpolatorNumba
+
+
+def test_numba_backend_falls_back_to_scipy_without_threadsafe_backend(monkeypatch):
+    def unavailable_backend():
+        raise ValueError('No threading layer could be loaded.')
+
+    monkeypatch.setattr(numba, 'get_num_threads', unavailable_backend)
+
+    with pytest.warns(
+        FutureWarning,
+        match="fallback is deprecated and will be an error in a future release",
+    ):
+        impl = to_wavelength._get_interpolator_class(FrameUnwrapBackend.numba)
+
+    assert impl is InterpolatorScipy
+
+
+def test_numba_backend_falls_back_to_scipy_without_numba():
+    with (
+        patch.dict('sys.modules', {'numba': None}),
+        pytest.warns(FutureWarning, match="Numba is unavailable"),
+    ):
+        impl = to_wavelength._get_interpolator_class(FrameUnwrapBackend.numba)
+
+    assert impl is InterpolatorScipy
+
+
+def test_numba_backend_falls_back_to_scipy_for_workqueue(monkeypatch):
+    monkeypatch.setattr(numba, 'get_num_threads', lambda: 1)
+    monkeypatch.setattr(numba, 'threading_layer', lambda: 'workqueue')
+
+    with pytest.warns(FutureWarning, match="thread-safe threading layer"):
+        impl = to_wavelength._get_interpolator_class(FrameUnwrapBackend.numba)
+
+    assert impl is InterpolatorScipy

packages/essreduce/tests/unwrap/workflow_test.py

@@ -22,6 +22,7 @@
     SampleRun,
 )
 from ess.reduce.unwrap import (
+    FrameUnwrapBackend,
     GenericUnwrapWorkflow,
     fakes,
     simulate_chopper_cascade_using_tof,
@@ -33,8 +34,17 @@
 Monitor0 = NewType("Monitor0", int)
 
 
+def test_GenericUnwrapWorkflow_defaults_to_numba_backend():
+    wf = GenericUnwrapWorkflow(run_types=[SampleRun], monitor_types=[])
+
+    assert wf.compute(FrameUnwrapBackend) == FrameUnwrapBackend.numba
+
+
 def _make_workflow(
-    wavelength_from, *, keep_event_time_offset=False
+    wavelength_from,
+    *,
+    backend: FrameUnwrapBackend = FrameUnwrapBackend.numba,
+    keep_event_time_offset=False,
 ) -> sciline.Pipeline:
     sizes = {'detector_number': 10}
     detector_geometry = sc.DataArray(
@@ -88,6 +98,7 @@ def _make_workflow(
         monitor_types=[Monitor0],
         wavelength_from=wavelength_from,
     )
+    wf[FrameUnwrapBackend] = backend
     wf[NeXusDetectorName] = "detector"
     wf[NeXusName[Monitor0]] = "monitor"
     wf[unwrap.LookupTableRelativeErrorThreshold] = {
@@ -122,10 +133,11 @@ def simulation_results_psc_choppers():
 
 @pytest.mark.parametrize("wavelength_from", ["simulation", "analytical"])
 @pytest.mark.parametrize("detector_or_monitor", ["detector", "monitor"])
+@pytest.mark.parametrize("backend", ["numba", "scipy"])
 def test_GenericUnwrapWorkflow_computes_wavelength(
-    wavelength_from, detector_or_monitor, simulation_results_psc_choppers
+    wavelength_from, detector_or_monitor, backend, simulation_results_psc_choppers
 ):
-    wf = _make_workflow(wavelength_from=wavelength_from)
+    wf = _make_workflow(wavelength_from=wavelength_from, backend=backend)
 
     if wavelength_from == "simulation":
         wf[unwrap.SimulationResults[SampleRun]] = simulation_results_psc_choppers