From f95ba51004582e6cf0d8bb98e0bf1c6bd5a106f7 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 05:47:39 +0000 Subject: [PATCH 1/7] Align powder focussing 2theta bins with requested TwoThetaBins Focussing binned into 1024 uniform 2theta bins spanning the detector range, independently of the 2theta binning requested for the result. Grouping into the requested bins then assigned each focussing bin as a whole to the bin containing its center, so the number of bins per group varied periodically. For BEER this produced spikes of the order of the signal itself. The focussing grid is now built as a refinement of TwoThetaBins: every requested edge occurs exactly in the grid, and each requested bin is subdivided until the sub-bins reach the resolution needed to reconstruct wavelength for monitor normalization. Grouping is then an exact sum over whole sub-bins. Outside the requested range the grid continues at the base resolution so that integrating over 2theta still covers the full detector. TwoThetaBins gains a None default, meaning the data will not be grouped by 2theta, so workflows that only produce 1D d-spacing output are unaffected. Both focussing and grouping now convert between the units of the requested edges and of the two_theta coordinate. Requesting bins in degrees against a coordinate in radians previously raised. Fixes #670 --- .../essdiffraction/src/ess/beer/workflow.py | 2 + .../essdiffraction/src/ess/dream/workflows.py | 2 + .../essdiffraction/src/ess/powder/grouping.py | 109 +++++++++++++--- .../essdiffraction/src/ess/powder/types.py | 6 +- .../src/ess/snspowder/powgen/workflow.py | 2 + .../tests/powder/grouping_test.py | 118 ++++++++++++++++++ .../snspowder/powgen/powgen_reduction_test.py | 1 + 7 files changed, 219 insertions(+), 21 deletions(-) create mode 100644 packages/essdiffraction/tests/powder/grouping_test.py diff --git a/packages/essdiffraction/src/ess/beer/workflow.py b/packages/essdiffraction/src/ess/beer/workflow.py index 40eb096b3..9c4ba1282 100644 --- a/packages/essdiffraction/src/ess/beer/workflow.py +++ b/packages/essdiffraction/src/ess/beer/workflow.py @@ -12,6 +12,7 @@ CaveMonitor, EmptyCanRun, SampleRun, + TwoThetaBins, VanadiumRun, ) @@ -32,6 +33,7 @@ default_parameters = { CalibrationData: None, + TwoThetaBins: None, PulseLength: sc.scalar(0.003, unit='s'), DetectorBankSizes: { 'south_detector': {'y': 200, 'x': 500}, diff --git a/packages/essdiffraction/src/ess/dream/workflows.py b/packages/essdiffraction/src/ess/dream/workflows.py index 813ae228d..149a5e8b4 100644 --- a/packages/essdiffraction/src/ess/dream/workflows.py +++ b/packages/essdiffraction/src/ess/dream/workflows.py @@ -24,6 +24,7 @@ ReducerSoftware, SampleRun, TofMask, + TwoThetaBins, TwoThetaMask, VanadiumRun, WavelengthMask, @@ -158,6 +159,7 @@ def default_parameters() -> dict: KeepEvents[VanadiumRun]: KeepEvents[VanadiumRun](True), KeepEvents[EmptyCanRun]: KeepEvents[EmptyCanRun](True), TofMask: None, + TwoThetaBins: None, WavelengthMask: None, TwoThetaMask: None, CIFAuthors: CIFAuthors([]), diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index a9b58cf91..7b395ad18 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -19,27 +19,94 @@ def _reconstruct_wavelength( - dspacing_bins: DspacingBins, two_theta_bins: TwoThetaBins + dspacing_bins: sc.Variable, two_theta_bins: sc.Variable ) -> sc.Variable: dspacing = dspacing_bins two_theta = sc.midpoints(two_theta_bins) return (2 * dspacing * sc.sin(two_theta / 2)).to(unit='angstrom') +_BASE_TWO_THETA_RESOLUTION = 1024 +"""Number of bins used to cover the full two-theta range of a detector. + +The two-theta binning used for focussing must be fine enough that the wavelength +reconstructed from bin centers (see :func:`_reconstruct_wavelength`) is accurate +enough for a monitor normalization. This is independent of the two-theta binning +requested for the final result, see :func:`_focussing_two_theta_bins`. +""" + + +def _extend_edges(start: float, step: float, limit: float) -> np.ndarray: + """Ascending edges reaching from ``start`` (exclusive) past ``limit``.""" + n = max(int(np.ceil((limit - start) / step)), 0) + edges = start + step * np.arange(1, n + 1) + return edges if step > 0 else edges[::-1] + + +def _focussing_two_theta_bins( + two_theta: sc.Variable, requested: sc.Variable | None +) -> sc.Variable: + """Return the two-theta bin edges to use for focussing. + + The edges cover the full range of ``two_theta`` with bins no wider than that + range divided by :py:data:`_BASE_TWO_THETA_RESOLUTION`. + + If ``requested`` is given, its edges are a subset of the returned edges: each + requested bin is subdivided into equally wide sub-bins. Grouping the focussed + data into ``requested`` is then an exact sum over whole sub-bins. Without + this alignment, each focussing bin is assigned as a whole to the requested + bin containing its center, so the number of bins per group varies + periodically and produces large spikes in the result. + """ + lo = two_theta.nanmin() + hi = two_theta.nanmax() + # Make the upper edge inclusive of the largest two-theta value. + hi.value = np.nextafter(hi.value, np.inf) + base_width = ((hi - lo) / _BASE_TWO_THETA_RESOLUTION).value + if requested is None: + return sc.linspace( + 'two_theta', + start=lo, + stop=hi, + num=_BASE_TWO_THETA_RESOLUTION + 1, + unit=two_theta.unit, + ) + edges = requested.to(unit=two_theta.unit, dtype='float64').values + n_sub = np.ceil(np.diff(edges) / base_width).astype(int) + # Each piece starts at an exact requested edge, so those edges are preserved. + refined = [ + np.linspace(low, high, num=n + 1)[:-1] + for low, high, n in zip(edges[:-1], edges[1:], n_sub, strict=True) + ] + return sc.array( + dims=['two_theta'], + values=np.concatenate( + [ + _extend_edges(edges[0], -base_width, lo.value), + *refined, + edges[-1:], + _extend_edges(edges[-1], base_width, hi.value), + ] + ), + unit=two_theta.unit, + ) + + def focus_data_dspacing_and_two_theta( data: CorrectedDetector[RunType], dspacing_bins: DspacingBins, + two_theta_bins: TwoThetaBins, keep_events: KeepEvents[RunType], ) -> CorrectedDspacing[RunType]: """ Reduce the pixel-based data to d-spacing and two-theta dimensions. - The two-theta binning does not use :py:class:`TwoThetaBins` but instead - computes the two-theta bins from the 'two_theta' coordinate of the input data. This - is necessary to ensure that we have sufficiently high wavelength resolution when - performing a monitor normalization in a follow-up workflow step. If we were to use - :py:class:`TwoThetaBins` we would be influenced by and limited to the two-theta - binning the user requests for the end result, which may not be sufficient. + The two-theta binning is finer than :py:class:`TwoThetaBins` and covers the full + two-theta range of the detector, not only the requested range. Both are necessary + to have sufficient wavelength resolution when performing a monitor normalization + in a follow-up workflow step. The bins are nevertheless aligned with + :py:class:`TwoThetaBins` such that :func:`group_two_theta` can produce the + requested binning exactly, see :func:`_focussing_two_theta_bins`. Parameters ---------- @@ -48,6 +115,9 @@ def focus_data_dspacing_and_two_theta( 'two_theta' coordinates. dspacing_bins: The bins to use for the d-spacing dimension. + two_theta_bins: + The two-theta bins requested for the final result, or ``None`` if the data + will not be grouped by two-theta. keep_events: Whether to keep the events in the output. If `False`, the output will be histogrammed instead of binned. @@ -57,17 +127,7 @@ def focus_data_dspacing_and_two_theta( : The reduced data with 'dspacing' and 'two_theta' dimensions. """ - ttheta = data.coords['two_theta'] - ttheta_min = ttheta.nanmin() - ttheta_max = ttheta.nanmax() - ttheta_max.value = np.nextafter(ttheta_max.value, np.inf) - twotheta_bins = sc.linspace( - 'two_theta', - start=ttheta_min, - stop=ttheta_max, - num=1024, - unit=ttheta.unit, - ) + twotheta_bins = _focussing_two_theta_bins(data.coords['two_theta'], two_theta_bins) args = {twotheta_bins.dim: twotheta_bins, dspacing_bins.dim: dspacing_bins} if keep_events.value: result = data.bin(args) @@ -102,10 +162,19 @@ def group_two_theta( data: NormalizedDspacing[RunType], two_theta_bins: TwoThetaBins, ) -> FocussedDataDspacingTwoTheta[RunType]: - """Group the data by two-theta bins.""" + """Group the data by two-theta bins. + + ``data`` was focussed onto a finer two-theta grid that is aligned with + ``two_theta_bins``, see :func:`focus_data_dspacing_and_two_theta`. Grouping is + therefore an exact sum over whole sub-bins. + """ + if two_theta_bins is None: + raise ValueError("Cannot group by two-theta, no 'TwoThetaBins' were set.") if 'two_theta' not in data.dims: raise ValueError("Data does not have a 'two_theta' dimension.") - data = data.assign_coords(two_theta=sc.midpoints(data.coords['two_theta'])) + data = data.assign_coords( + two_theta=sc.midpoints(data.coords['two_theta']).to(unit=two_theta_bins.unit) + ) return FocussedDataDspacingTwoTheta[RunType]( data.groupby('two_theta', bins=two_theta_bins).nansum('two_theta') if data.bins is None diff --git a/packages/essdiffraction/src/ess/powder/types.py b/packages/essdiffraction/src/ess/powder/types.py index 81f51435d..b65361cdd 100644 --- a/packages/essdiffraction/src/ess/powder/types.py +++ b/packages/essdiffraction/src/ess/powder/types.py @@ -65,11 +65,15 @@ OutFilename = NewType("OutFilename", str) """Filename of the output.""" -TwoThetaBins = NewType("TwoThetaBins", sc.Variable) +TwoThetaBins = NewType("TwoThetaBins", sc.Variable | None) """Bin edges for grouping in 2theta. This is used by an alternative focussing step that groups detector pixels by scattering angle into bins given by these edges. + +The focussing step aligns its internal, finer 2theta binning with these edges, +so setting them also affects results that are integrated over 2theta. +``None`` means that the data will not be grouped by 2theta. """ UncertaintyBroadcastMode = _UncertaintyBroadcastMode diff --git a/packages/essdiffraction/src/ess/snspowder/powgen/workflow.py b/packages/essdiffraction/src/ess/snspowder/powgen/workflow.py index f6bef6baa..eec6a3e15 100644 --- a/packages/essdiffraction/src/ess/snspowder/powgen/workflow.py +++ b/packages/essdiffraction/src/ess/snspowder/powgen/workflow.py @@ -8,6 +8,7 @@ KeepEvents, NeXusDetectorName, SampleRun, + TwoThetaBins, VanadiumRun, ) @@ -23,6 +24,7 @@ def default_parameters() -> dict: KeepEvents[VanadiumRun]: KeepEvents[VanadiumRun](False), KeepEvents[EmptyCanRun]: KeepEvents[EmptyCanRun](True), NeXusDetectorName: "powgen_detector", + TwoThetaBins: None, } diff --git a/packages/essdiffraction/tests/powder/grouping_test.py b/packages/essdiffraction/tests/powder/grouping_test.py new file mode 100644 index 000000000..f9745b655 --- /dev/null +++ b/packages/essdiffraction/tests/powder/grouping_test.py @@ -0,0 +1,118 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import numpy as np +import pytest +import scipp as sc +from ess.powder.grouping import ( + focus_data_dspacing_and_two_theta, + group_two_theta, + integrate_two_theta, +) +from ess.powder.types import KeepEvents, SampleRun + +DSPACING_BINS = sc.linspace('dspacing', 1.0, 2.0, num=11, unit='angstrom') +"""Coarse d-spacing bins; the tests are only concerned with two-theta.""" + + +@pytest.fixture +def detector() -> sc.DataArray: + """Events spread smoothly over the full two-theta range of a wide detector.""" + rng = np.random.default_rng(90210) + n_pixel = 5000 + n_event = 20 + two_theta = sc.array( + dims=['pixel'], + values=np.deg2rad(np.linspace(10.0, 170.0, n_pixel)), + unit='rad', + ) + events = sc.DataArray( + sc.ones(dims=['event'], shape=[n_pixel * n_event], unit='counts'), + coords={ + 'dspacing': sc.array( + dims=['event'], + values=rng.uniform(1.0, 2.0, n_pixel * n_event), + unit='angstrom', + ), + 'pixel': sc.array( + dims=['event'], values=np.repeat(np.arange(n_pixel), n_event) + ), + }, + ) + return events.group('pixel').drop_coords('pixel').assign_coords(two_theta=two_theta) + + +def _focus_and_group( + detector: sc.DataArray, two_theta_bins: sc.Variable, *, keep_events: bool +) -> sc.DataArray: + focussed = focus_data_dspacing_and_two_theta( + detector, + DSPACING_BINS, + two_theta_bins, + KeepEvents[SampleRun](keep_events), + ) + grouped = group_two_theta(focussed, two_theta_bins) + return grouped if grouped.bins is None else grouped.hist() + + +@pytest.mark.parametrize('keep_events', [True, False]) +@pytest.mark.parametrize( + 'two_theta_bins', + [ + sc.linspace('two_theta', 75.0, 105.0, num=180, unit='deg'), + sc.linspace('two_theta', 0.8, 2.4, num=17, unit='rad'), + sc.array(dims=['two_theta'], values=[0.5, 0.6, 1.5, 1.55, 2.9], unit='rad'), + ], + ids=['many-narrow-deg', 'few-wide-rad', 'non-uniform-rad'], +) +def test_group_two_theta_matches_direct_histogram( + detector, two_theta_bins, keep_events +): + """Focussing must not redistribute counts between requested two-theta bins. + + Focussing bins that are not aligned with the requested bins get assigned to + whichever requested bin contains their center. The number of bins per group + then varies periodically, producing large spikes. + """ + result = _focus_and_group(detector, two_theta_bins, keep_events=keep_events) + expected = detector.hist( + two_theta=two_theta_bins.to(unit=detector.coords['two_theta'].unit), + dspacing=DSPACING_BINS, + ) + assert sc.allclose(result.data, expected.data) + + +def test_group_two_theta_preserves_requested_bins(detector): + two_theta_bins = sc.linspace('two_theta', 75.0, 105.0, num=180, unit='deg') + result = _focus_and_group(detector, two_theta_bins, keep_events=False) + assert sc.identical(result.coords['two_theta'], two_theta_bins) + + +@pytest.mark.parametrize('keep_events', [True, False]) +def test_integrate_two_theta_covers_full_detector_range(detector, keep_events): + """The requested bins may cover only part of the detector.""" + two_theta_bins = sc.linspace('two_theta', 75.0, 105.0, num=180, unit='deg') + focussed = focus_data_dspacing_and_two_theta( + detector, DSPACING_BINS, two_theta_bins, KeepEvents[SampleRun](keep_events) + ) + result = integrate_two_theta(focussed) + if result.bins is not None: + result = result.hist() + expected = detector.hist(dspacing=DSPACING_BINS).sum('pixel') + assert sc.allclose(result.data, expected.data) + + +def test_focussing_without_requested_bins_keeps_all_counts(detector): + focussed = focus_data_dspacing_and_two_theta( + detector, DSPACING_BINS, None, KeepEvents[SampleRun](False) + ) + assert sc.allclose( + focussed.sum().data, detector.hist(dspacing=DSPACING_BINS).sum().data + ) + + +def test_group_two_theta_without_requested_bins_raises(detector): + focussed = focus_data_dspacing_and_two_theta( + detector, DSPACING_BINS, None, KeepEvents[SampleRun](False) + ) + with pytest.raises(ValueError, match='TwoThetaBins'): + group_two_theta(focussed, None) diff --git a/packages/essdiffraction/tests/snspowder/powgen/powgen_reduction_test.py b/packages/essdiffraction/tests/snspowder/powgen/powgen_reduction_test.py index 1e374bd19..3ff09ac1c 100644 --- a/packages/essdiffraction/tests/snspowder/powgen/powgen_reduction_test.py +++ b/packages/essdiffraction/tests/snspowder/powgen/powgen_reduction_test.py @@ -51,6 +51,7 @@ def params(): (x < sc.scalar(0.0, unit="us").to(unit=elem_unit(x))) | (x > sc.scalar(16666.67, unit="us").to(unit=elem_unit(x))) ), + TwoThetaBins: None, TwoThetaMask: None, WavelengthMask: None, GravityVector: sc.vector(value=[0, -1, 0]) * sc.constants.g, From a06f4be602e548a90edf41bf9f2239deda0333ab Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 05:56:53 +0000 Subject: [PATCH 2/7] Document what sets the focussing 2theta resolution The number only matters for histogrammed data normalized by a monitor histogram: the monitor is looked up once per (d-spacing, two-theta) bin and the lookup is piecewise constant on the monitor bins. Data that keeps its events carries its own wavelength and is unaffected. --- .../essdiffraction/src/ess/powder/grouping.py | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index 7b395ad18..7867ea5e7 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -29,10 +29,24 @@ def _reconstruct_wavelength( _BASE_TWO_THETA_RESOLUTION = 1024 """Number of bins used to cover the full two-theta range of a detector. -The two-theta binning used for focussing must be fine enough that the wavelength -reconstructed from bin centers (see :func:`_reconstruct_wavelength`) is accurate -enough for a monitor normalization. This is independent of the two-theta binning +This sets the resolution of the wavelength reconstructed from bin centers, see +:func:`_reconstruct_wavelength`. It is independent of the two-theta binning requested for the final result, see :func:`_focussing_two_theta_bins`. + +It only matters when focussing produces a histogram *and* the run is normalized +by a monitor histogram. The monitor is then looked up once per (d-spacing, +two-theta) bin and the lookup is piecewise constant on the bins of the monitor, +so the wavelength spread within a bin should stay below the width of a monitor +bin. Since + +.. math:: + + \\frac{\\Delta\\lambda}{\\lambda} + = \\frac{\\Delta 2\\theta}{2} \\cot\\theta + +this is most demanding at small scattering angles. If the focussed data keeps +its events, they carry their own wavelength and this number has no influence on +the result. """ From 0d15d4a1eead1138575a1a671ed400cec831dba7 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 06:01:28 +0000 Subject: [PATCH 3/7] Build focussing 2theta bins without a Python loop The bins are recomputed for every chunk in a stream processor, since they depend on the two-theta coordinate of the chunked detector data. Subdividing the requested bins one at a time cost 1.6 ms for 180 requested bins, a third of the cost of the histogramming step it feeds. Computing the sub-bin positions in one go makes the cost independent of the number of requested bins, at around 0.4 ms, half of which is the min/max over the two-theta coordinate. --- .../essdiffraction/src/ess/powder/grouping.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index 7867ea5e7..c927ad231 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -86,18 +86,19 @@ def _focussing_two_theta_bins( unit=two_theta.unit, ) edges = requested.to(unit=two_theta.unit, dtype='float64').values - n_sub = np.ceil(np.diff(edges) / base_width).astype(int) - # Each piece starts at an exact requested edge, so those edges are preserved. - refined = [ - np.linspace(low, high, num=n + 1)[:-1] - for low, high, n in zip(edges[:-1], edges[1:], n_sub, strict=True) - ] + widths = np.diff(edges) + n_sub = np.ceil(widths / base_width).astype(int) + # Index of the requested bin each sub-bin belongs to, and its position within it. + # The first sub-bin of a requested bin reproduces its lower edge exactly. + offset = np.concatenate([[0], np.cumsum(n_sub)]) + index = np.repeat(np.arange(len(n_sub)), n_sub) + position = (np.arange(offset[-1]) - offset[index]) / n_sub[index] return sc.array( dims=['two_theta'], values=np.concatenate( [ _extend_edges(edges[0], -base_width, lo.value), - *refined, + edges[index] + position * widths[index], edges[-1:], _extend_edges(edges[-1], base_width, hi.value), ] From 44f661d481aafb4ceb7eeb725ae6c42785fba061 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 06:39:10 +0000 Subject: [PATCH 4/7] Keep the focussing 2theta bins linearly spaced scipp::bin and scipp::hist have a fast path for linearly spaced bin edges. Continuing the grid beyond the requested range at the base width rather than at the width of the sub-bins broke that, costing 12% of the histogramming step. Two-theta is a per-pixel coordinate, so the binary search is amortized over pixels rather than events; the same loss is a factor of 8-10 where the binning coordinate is an event coordinate. Extending at the width of the outermost sub-bins restores linear spacing whenever the requested bins are linearly spaced. A sub-bin is never narrower than half the base width, so this at most doubles the number of bins outside the requested range. --- packages/essdiffraction/src/ess/powder/grouping.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index c927ad231..917f83472 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -71,6 +71,12 @@ def _focussing_two_theta_bins( this alignment, each focussing bin is assigned as a whole to the requested bin containing its center, so the number of bins per group varies periodically and produces large spikes in the result. + + Beyond the requested range the grid continues at the width of the outermost + sub-bins rather than the base width. This keeps the edges linearly spaced + when ``requested`` is, which :func:`scipp.bin` and :func:`scipp.hist` can + exploit. Since a sub-bin is never narrower than half the base width, it at + most doubles the number of bins needed to reach the ends of the range. """ lo = two_theta.nanmin() hi = two_theta.nanmax() @@ -93,14 +99,15 @@ def _focussing_two_theta_bins( offset = np.concatenate([[0], np.cumsum(n_sub)]) index = np.repeat(np.arange(len(n_sub)), n_sub) position = (np.arange(offset[-1]) - offset[index]) / n_sub[index] + sub_widths = widths / n_sub return sc.array( dims=['two_theta'], values=np.concatenate( [ - _extend_edges(edges[0], -base_width, lo.value), + _extend_edges(edges[0], -sub_widths[0], lo.value), edges[index] + position * widths[index], edges[-1:], - _extend_edges(edges[-1], base_width, hi.value), + _extend_edges(edges[-1], sub_widths[-1], hi.value), ] ), unit=two_theta.unit, From 56c69990f9c92ec5f6d56992205d8ba3bdf08961 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 06:42:30 +0000 Subject: [PATCH 5/7] Say why linearly spaced focussing bins are worth having The docstring stated that scipp can exploit linear spacing without saying what for. --- packages/essdiffraction/src/ess/powder/grouping.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index 917f83472..33d235bc4 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -74,9 +74,11 @@ def _focussing_two_theta_bins( Beyond the requested range the grid continues at the width of the outermost sub-bins rather than the base width. This keeps the edges linearly spaced - when ``requested`` is, which :func:`scipp.bin` and :func:`scipp.hist` can - exploit. Since a sub-bin is never narrower than half the base width, it at - most doubles the number of bins needed to reach the ends of the range. + when ``requested`` is, which is faster to bin into: :func:`scipp.bin` and + :func:`scipp.hist` can then compute bin indices directly instead of + searching for them. Since a sub-bin is never narrower than half the base + width, it at most doubles the number of bins needed to reach the ends of + the range. """ lo = two_theta.nanmin() hi = two_theta.nanmax() From 8509080382e58b75f0a61e71f7ec32e94def69ab Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 09:43:10 +0000 Subject: [PATCH 6/7] Handle degenerate two-theta bins and check grouping alignment Responses to review: - Drop the redundant `max` in `_extend_edges`; `np.arange` already yields no edges when the requested bins extend past the detector range. Order its arguments as start/limit/step. - Subdivide a zero-width requested bin into a single zero-width sub-bin instead of none. It comes out empty, as it would from `scipp.hist`, whereas previously an interior one produced a 0/0 warning and an outermost one a NaN cast error. Descending bins now reach `scipp.bin` and `scipp.hist`, which raise `BinEdgeError`. - Group in the unit of the data and label the result with the requested edges instead of round-tripping the data through the requested unit. The requested edges are then bit-identical to the grid built by focussing, so `group_two_theta` can check exactly that both steps were given the same bins. Without it, calling the two functions directly with different bins silently returns the aliasing this branch removes. - Say accurately what `TwoThetaBins = None` means. --- .../essdiffraction/src/ess/powder/grouping.py | 53 ++++++++++++++----- .../essdiffraction/src/ess/powder/types.py | 3 +- .../tests/powder/grouping_test.py | 34 +++++++++++- 3 files changed, 75 insertions(+), 15 deletions(-) diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index 33d235bc4..af4d7c87d 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -50,9 +50,13 @@ def _reconstruct_wavelength( """ -def _extend_edges(start: float, step: float, limit: float) -> np.ndarray: - """Ascending edges reaching from ``start`` (exclusive) past ``limit``.""" - n = max(int(np.ceil((limit - start) / step)), 0) +def _extend_edges(start: float, limit: float, step: float) -> np.ndarray: + """Ascending edges reaching from ``start`` (exclusive) past ``limit``. + + Empty if ``limit`` has already been passed, which happens when the requested + bins extend beyond the range of the detector. + """ + n = int(np.ceil((limit - start) / step)) edges = start + step * np.arange(1, n + 1) return edges if step > 0 else edges[::-1] @@ -95,21 +99,27 @@ def _focussing_two_theta_bins( ) edges = requested.to(unit=two_theta.unit, dtype='float64').values widths = np.diff(edges) - n_sub = np.ceil(widths / base_width).astype(int) + # Zero-width bins are kept as zero-width sub-bins; like scipp's own binning, + # they end up empty. + n_sub = np.maximum(np.ceil(widths / base_width), 1).astype(int) # Index of the requested bin each sub-bin belongs to, and its position within it. # The first sub-bin of a requested bin reproduces its lower edge exactly. offset = np.concatenate([[0], np.cumsum(n_sub)]) index = np.repeat(np.arange(len(n_sub)), n_sub) position = (np.arange(offset[-1]) - offset[index]) / n_sub[index] sub_widths = widths / n_sub + # A zero-width outermost bin provides no step to continue the grid with. The + # grid is not linearly spaced in that case anyway, so fall back to the base width. + low_step = sub_widths[0] if sub_widths[0] > 0 else base_width + high_step = sub_widths[-1] if sub_widths[-1] > 0 else base_width return sc.array( dims=['two_theta'], values=np.concatenate( [ - _extend_edges(edges[0], -sub_widths[0], lo.value), + _extend_edges(edges[0], lo.value, -low_step), edges[index] + position * widths[index], edges[-1:], - _extend_edges(edges[-1], sub_widths[-1], hi.value), + _extend_edges(edges[-1], hi.value, high_step), ] ), unit=two_theta.unit, @@ -182,6 +192,18 @@ def integrate_two_theta( ) +def _check_aligned(edges: sc.Variable, requested: sc.Variable) -> None: + """Raise unless every requested edge occurs in ``edges``.""" + grid = edges.values + found = np.clip(np.searchsorted(grid, requested.values), 0, len(grid) - 1) + if not np.array_equal(grid[found], requested.values): + raise ValueError( + 'The two-theta binning of the data is not aligned with the requested ' + 'bins. The data must be focussed with the same TwoThetaBins, see ' + 'focus_data_dspacing_and_two_theta.' + ) + + def group_two_theta( data: NormalizedDspacing[RunType], two_theta_bins: TwoThetaBins, @@ -190,19 +212,26 @@ def group_two_theta( ``data`` was focussed onto a finer two-theta grid that is aligned with ``two_theta_bins``, see :func:`focus_data_dspacing_and_two_theta`. Grouping is - therefore an exact sum over whole sub-bins. + therefore an exact sum over whole sub-bins. This requires that both steps were + given the same bins; otherwise, an exception is raised. """ if two_theta_bins is None: raise ValueError("Cannot group by two-theta, no 'TwoThetaBins' were set.") if 'two_theta' not in data.dims: raise ValueError("Data does not have a 'two_theta' dimension.") - data = data.assign_coords( - two_theta=sc.midpoints(data.coords['two_theta']).to(unit=two_theta_bins.unit) + edges = data.coords['two_theta'] + # Grouping in the unit of the data keeps the requested edges bit-identical to the + # grid built by focussing, so that the alignment can be checked exactly. + bins = two_theta_bins.to(unit=edges.unit, dtype='float64') + _check_aligned(edges, bins) + data = data.assign_coords(two_theta=sc.midpoints(edges)) + grouped = ( + data.groupby('two_theta', bins=bins).nansum('two_theta') + if data.bins is None + else data.bin(two_theta=bins) ) return FocussedDataDspacingTwoTheta[RunType]( - data.groupby('two_theta', bins=two_theta_bins).nansum('two_theta') - if data.bins is None - else data.bin(two_theta=two_theta_bins) + grouped.assign_coords(two_theta=two_theta_bins) ) diff --git a/packages/essdiffraction/src/ess/powder/types.py b/packages/essdiffraction/src/ess/powder/types.py index b65361cdd..5cc18dea6 100644 --- a/packages/essdiffraction/src/ess/powder/types.py +++ b/packages/essdiffraction/src/ess/powder/types.py @@ -73,7 +73,8 @@ The focussing step aligns its internal, finer 2theta binning with these edges, so setting them also affects results that are integrated over 2theta. -``None`` means that the data will not be grouped by 2theta. +``None`` means that no 2theta grouping is available: the workflow can then only +produce results integrated over 2theta and raises if grouped data is requested. """ UncertaintyBroadcastMode = _UncertaintyBroadcastMode diff --git a/packages/essdiffraction/tests/powder/grouping_test.py b/packages/essdiffraction/tests/powder/grouping_test.py index f9745b655..eecfd5909 100644 --- a/packages/essdiffraction/tests/powder/grouping_test.py +++ b/packages/essdiffraction/tests/powder/grouping_test.py @@ -51,7 +51,7 @@ def _focus_and_group( KeepEvents[SampleRun](keep_events), ) grouped = group_two_theta(focussed, two_theta_bins) - return grouped if grouped.bins is None else grouped.hist() + return grouped.hist() if grouped.is_binned else grouped @pytest.mark.parametrize('keep_events', [True, False]) @@ -61,8 +61,16 @@ def _focus_and_group( sc.linspace('two_theta', 75.0, 105.0, num=180, unit='deg'), sc.linspace('two_theta', 0.8, 2.4, num=17, unit='rad'), sc.array(dims=['two_theta'], values=[0.5, 0.6, 1.5, 1.55, 2.9], unit='rad'), + sc.linspace('two_theta', 0.0, 180.0, num=19, unit='deg'), + sc.array(dims=['two_theta'], values=[0.5, 0.5, 1.5, 2.9, 2.9], unit='rad'), + ], + ids=[ + 'many-narrow-deg', + 'few-wide-rad', + 'non-uniform-rad', + 'beyond-detector-deg', + 'zero-width-rad', ], - ids=['many-narrow-deg', 'few-wide-rad', 'non-uniform-rad'], ) def test_group_two_theta_matches_direct_histogram( detector, two_theta_bins, keep_events @@ -110,6 +118,28 @@ def test_focussing_without_requested_bins_keeps_all_counts(detector): ) +@pytest.mark.parametrize('keep_events', [True, False]) +def test_focussing_with_descending_bins_raises(detector, keep_events): + two_theta_bins = sc.array(dims=['two_theta'], values=[0.5, 1.5, 1.0], unit='rad') + with pytest.raises(sc.BinEdgeError): + focus_data_dspacing_and_two_theta( + detector, DSPACING_BINS, two_theta_bins, KeepEvents[SampleRun](keep_events) + ) + + +@pytest.mark.parametrize('keep_events', [True, False]) +def test_group_two_theta_with_misaligned_bins_raises(detector, keep_events): + focussed = focus_data_dspacing_and_two_theta( + detector, + DSPACING_BINS, + sc.linspace('two_theta', 0.8, 2.4, num=17, unit='rad'), + KeepEvents[SampleRun](keep_events), + ) + other_bins = sc.linspace('two_theta', 0.8, 2.4, num=23, unit='rad') + with pytest.raises(ValueError, match='aligned'): + group_two_theta(focussed, other_bins) + + def test_group_two_theta_without_requested_bins_raises(detector): focussed = focus_data_dspacing_and_two_theta( detector, DSPACING_BINS, None, KeepEvents[SampleRun](False) From 5e49634137d0096a45ca8a32a0013e789c125baf Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 10:59:47 +0000 Subject: [PATCH 7/7] Test zero-width two-theta bins in every position The previous case had zero-width first and last bins but no interior one, which is the case that produced the 0/0 warning. --- packages/essdiffraction/src/ess/powder/grouping.py | 4 ++-- packages/essdiffraction/tests/powder/grouping_test.py | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/packages/essdiffraction/src/ess/powder/grouping.py b/packages/essdiffraction/src/ess/powder/grouping.py index af4d7c87d..612bf7918 100644 --- a/packages/essdiffraction/src/ess/powder/grouping.py +++ b/packages/essdiffraction/src/ess/powder/grouping.py @@ -99,8 +99,8 @@ def _focussing_two_theta_bins( ) edges = requested.to(unit=two_theta.unit, dtype='float64').values widths = np.diff(edges) - # Zero-width bins are kept as zero-width sub-bins; like scipp's own binning, - # they end up empty. + # A zero-width bin gets a single zero-width sub-bin instead of none, which keeps + # the sub-bin width well-defined. It comes out empty, as it would from scipp. n_sub = np.maximum(np.ceil(widths / base_width), 1).astype(int) # Index of the requested bin each sub-bin belongs to, and its position within it. # The first sub-bin of a requested bin reproduces its lower edge exactly. diff --git a/packages/essdiffraction/tests/powder/grouping_test.py b/packages/essdiffraction/tests/powder/grouping_test.py index eecfd5909..8b03349f8 100644 --- a/packages/essdiffraction/tests/powder/grouping_test.py +++ b/packages/essdiffraction/tests/powder/grouping_test.py @@ -62,6 +62,7 @@ def _focus_and_group( sc.linspace('two_theta', 0.8, 2.4, num=17, unit='rad'), sc.array(dims=['two_theta'], values=[0.5, 0.6, 1.5, 1.55, 2.9], unit='rad'), sc.linspace('two_theta', 0.0, 180.0, num=19, unit='deg'), + sc.array(dims=['two_theta'], values=[0.5, 1.5, 1.5, 2.9], unit='rad'), sc.array(dims=['two_theta'], values=[0.5, 0.5, 1.5, 2.9, 2.9], unit='rad'), ], ids=[ @@ -69,7 +70,8 @@ def _focus_and_group( 'few-wide-rad', 'non-uniform-rad', 'beyond-detector-deg', - 'zero-width-rad', + 'zero-width-interior-rad', + 'zero-width-outermost-rad', ], ) def test_group_two_theta_matches_direct_histogram(