[ESSSANS] fix(sans): snap wavelength band edges onto WavelengthBins

packages/esssans/src/ess/sans/normalization.py

@@ -1,5 +1,6 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
+import numpy as np
 import scipp as sc
 import scippnexus as snx
 from scipp.core import concepts
@@ -281,6 +282,20 @@ def process_wavelength_bands(
 
     The final bands must have a size of 2 in the wavelength dimension, defining a start
     and an end wavelength.
+
+    The band edges are snapped onto ``wavelength_bins``, moving each edge by at most
+    half a bin width. :func:`_reduce` selects a band from event data by re-binning onto
+    the band edges but from dense data by slicing whole bins, and only bands aligned
+    with the bins make the two select the same wavelength range. Unaligned bands leave
+    the I(Q) numerator and denominator covering different ranges, biasing each band by
+    several percent when bands are only a few bins wide. See
+    :py:class:`ess.sans.types.ProcessedWavelengthBands`.
+
+    Raises
+    ------
+    ValueError
+        If the bands are not two wavelength values per band, fall outside
+        ``wavelength_bins``, are not ordered, or are narrower than a wavelength bin.
     """
     if wavelength_bands is None:
         wavelength_bands = sc.concat(
@@ -301,7 +316,49 @@ def process_wavelength_bands(
             'defining a start and an end wavelength, '
             f'got {wavelength_bands.sizes["wavelength"]}.'
         )
-    return wavelength_bands
+    wavelength_bands = wavelength_bands.to(
+        unit=wavelength_bins.unit, dtype='float64', copy=False
+    )
+    lo = wavelength_bins.min()
+    hi = wavelength_bins.max()
+    # A NaN edge compares False against both bounds, so it is rejected here rather
+    # than snapping to an arbitrary bin.
+    if not sc.all((wavelength_bands >= lo) & (wavelength_bands <= hi)).value:
+        raise ValueError(
+            f'Wavelength bands must lie within the wavelength bins [{lo:c}, {hi:c}], '
+            f'got {wavelength_bands}'
+        )
+    if not sc.all(_widths(wavelength_bands) > sc.scalar(0.0, unit=lo.unit)).value:
+        raise ValueError(
+            f'Wavelength bands must start before they end, got {wavelength_bands}'
+        )
+    bands = _snap_to_bins(wavelength_bands, wavelength_bins)
+    widths = _widths(bands)
+    collapsed = int((widths <= sc.scalar(0.0, unit=lo.unit)).sum().value)
+    if collapsed:
+        raise ValueError(
+            f'{collapsed} of {widths.size} wavelength bands are narrower than a '
+            'wavelength bin and collapse to zero width when snapped onto it. Widen the '
+            'bands or refine WavelengthBins.'
+        )
+    return bands
+
+
+def _widths(bands: sc.Variable) -> sc.Variable:
+    return bands['wavelength', 1] - bands['wavelength', 0]
+
+
+def _snap_to_bins(bands: sc.Variable, bins: sc.Variable) -> sc.Variable:
+    """Move every band edge to the nearest wavelength-bin edge."""
+    edges = bins.values
+    upper = np.clip(np.searchsorted(edges, bands.values), 1, len(edges) - 1)
+    lower = upper - 1
+    nearest = np.where(
+        bands.values - edges[lower] <= edges[upper] - bands.values, lower, upper
+    )
+    # Values are taken from ``bins`` rather than rounded, so they are bit-identical to
+    # the bin edges. Slicing rules that differ in how they treat a boundary then agree.
+    return sc.array(dims=bands.dims, values=edges[nearest], unit=bins.unit)
 
 
 def _normalize(
@@ -413,6 +470,11 @@ def _reduce(part: sc.DataArray, /, *, bands: ProcessedWavelengthBands) -> sc.Dat
         # If in event mode the desired wavelength binning has not been applied, we need
         # it for splitting by bands, or restricting the range in case of a single band.
         part = part.bin(wavelength=sc.sort(bands.flatten(to=wav), wav))
+    # Selection by label slicing means different things depending on the wavelength
+    # coord: bin edges select every overlapping bin, midpoints (as carried by the
+    # denominator, see `norm_detector_term_denominator`) select by nearest bin. The
+    # three cases coincide only because `process_wavelength_bands` aligned the bands
+    # with the bins.
     parts = [
         _do_reduce(part[wav, wav_range[0] : wav_range[1]])
         for wav_range in sc.collapse(bands, keep=wav).values()

packages/esssans/src/ess/sans/types.py

@@ -116,7 +116,13 @@ class TransmissionRun(Generic[ScatteringRunType]):
 ProcessedWavelengthBands = NewType('ProcessedWavelengthBands', sc.Variable)
 """Processed wavelength bands, as a two-dimensional variable, with the first dimension
 being the band index and the second dimension being the wavelength. For each band, there
-must be two wavelength values defining the start and end wavelength of the band."""
+must be two wavelength values defining the start and end wavelength of the band.
+
+Band edges coincide exactly with values of :py:obj:`WavelengthBins`. Reducing over a
+band relies on this, so that event data, dense data with a bin-edge wavelength coord and
+dense data with a midpoint wavelength coord all select the same wavelength range.
+Construct this via :py:func:`ess.sans.normalization.process_wavelength_bands`, which
+establishes the alignment; hand-built bands silently bias I(Q) per band."""
 
 
 QBins = NewType('QBins', sc.Variable)

packages/esssans/tests/normalization_test.py

@@ -188,3 +188,143 @@ def test_transmission_fraction():
             direct_transmission_monitor=direct_transmission_monitor,
         ).data,
     )
+
+
+@pytest.fixture
+def wavelength_bins() -> sc.Variable:
+    return sc.linspace('wavelength', 1.0, 13.0, num=51, unit='angstrom')
+
+
+@pytest.fixture
+def q_bins() -> sc.Variable:
+    return sc.linspace('Q', 0.0, 1.0, num=2, unit='1/angstrom')
+
+
+def _flat_density(bins: sc.Variable, q_bins: sc.Variable) -> sc.DataArray:
+    """Events with a uniform density of one count per angstrom.
+
+    Reducing this over a band yields the width of the band in angstrom, so any
+    discrepancy between two representations is the discrepancy in the wavelength range
+    they select. The event count must stay large compared to the number of bands, else
+    the discretization of the uniform placement dominates the comparison tolerance.
+    """
+    n = 120_000
+    lo, hi = bins.min().value, bins.max().value
+    return sc.DataArray(
+        sc.full(dims=['event'], shape=[n], value=(hi - lo) / n, unit='counts'),
+        coords={
+            'wavelength': sc.array(
+                dims=['event'],
+                values=np.linspace(lo, hi, n, endpoint=False),
+                unit='angstrom',
+            ),
+            'Q': sc.full(dims=['event'], shape=[n], value=0.5, unit='1/angstrom'),
+        },
+    ).bin(Q=q_bins, wavelength=bins)
+
+
+def _as_midpoints(histogram: sc.DataArray) -> sc.DataArray:
+    """Replace the wavelength bin edges by midpoints.
+
+    The I(Q) denominator is dense in this form, because computing Q requires one
+    wavelength value per bin. See :func:`normalization.norm_detector_term_denominator`.
+    """
+    return histogram.assign_coords(
+        wavelength=sc.midpoints(histogram.coords['wavelength'])
+    )
+
+
+@pytest.mark.parametrize('nbands', [7, 10, 13])
+def test_reduce_q_selects_same_wavelength_range_for_all_representations(
+    wavelength_bins, q_bins, nbands
+):
+    bands = normalization.process_wavelength_bands(
+        sc.linspace(
+            'wavelength',
+            wavelength_bins.min().value,
+            wavelength_bins.max().value,
+            num=nbands + 1,
+            unit='angstrom',
+        ),
+        wavelength_bins,
+    )
+    events = _flat_density(wavelength_bins, q_bins)
+    representations = {
+        'events': events,
+        'bin_edges': events.hist(),
+        'midpoints': _as_midpoints(events.hist()),
+    }
+    reduced = {
+        name: normalization.reduce_q(data, bands=bands)
+        for name, data in representations.items()
+    }
+    reduced['events'] = reduced['events'].hist()
+    for name, result in reduced.items():
+        assert sc.allclose(result.data, reduced['events'].data, rtol=sc.scalar(1e-4)), (
+            name
+        )
+
+
+def test_process_wavelength_bands_returns_exact_bin_edges(wavelength_bins):
+    bands = sc.linspace('wavelength', 1.0, 13.0, num=11, unit='angstrom')
+    processed = normalization.process_wavelength_bands(bands, wavelength_bins)
+    assert sc.identical(
+        processed,
+        sc.concat([bands[:-1], bands[1:]], dim='x').rename(
+            x='wavelength', wavelength='band'
+        ),
+    )
+
+
+def test_process_wavelength_bands_snaps_unaligned_bands_onto_bins(wavelength_bins):
+    bands = sc.linspace('wavelength', 1.0, 13.0, num=8, unit='angstrom')
+    processed = normalization.process_wavelength_bands(bands, wavelength_bins)
+    assert set(np.unique(processed.values)) <= set(wavelength_bins.values)
+    half_width = 0.5 * (wavelength_bins[1] - wavelength_bins[0]).value
+    assert np.all(np.abs(np.unique(processed.values) - bands.values) <= half_width)
+
+
+def test_process_wavelength_bands_is_idempotent(wavelength_bins):
+    """`direct_beam` feeds already-processed bands back in as `WavelengthBands`."""
+    bands = sc.linspace('wavelength', 1.0, 13.0, num=8, unit='angstrom')
+    once = normalization.process_wavelength_bands(bands, wavelength_bins)
+    assert sc.identical(
+        normalization.process_wavelength_bands(once, wavelength_bins), once
+    )
+
+
+def test_process_wavelength_bands_snaps_overlapping_bands(wavelength_bins):
+    edges = sc.linspace('band', 1.0, 13.0, num=12, unit='angstrom')
+    bands = sc.concat([edges[:-2], edges[2:]], dim='wavelength').transpose()
+    processed = normalization.process_wavelength_bands(bands, wavelength_bins)
+    assert processed.dims == bands.dims
+    assert set(np.unique(processed.values)) <= set(wavelength_bins.values)
+    bin_width = (wavelength_bins[1] - wavelength_bins[0]).value
+    processed_hi = processed['wavelength', 1]['band', :-1]
+    processed_lo = processed['wavelength', 0]['band', 1:]
+    overlap = (processed_hi - processed_lo).values
+    bands_hi = bands['wavelength', 1]['band', :-1]
+    bands_lo = bands['wavelength', 0]['band', 1:]
+    expected = (bands_hi - bands_lo).values
+    assert np.all(np.abs(overlap - expected) <= bin_width)
+
+
+def test_process_wavelength_bands_raises_if_bands_narrower_than_bins(wavelength_bins):
+    bands = sc.linspace('wavelength', 1.0, 13.0, num=201, unit='angstrom')
+    with pytest.raises(ValueError, match='collapse to zero width'):
+        normalization.process_wavelength_bands(bands, wavelength_bins)
+
+
+@pytest.mark.parametrize(
+    'values', [[0.5, 5.0], [5.0, 20.0], [-3.0, -1.0], [2.0, float('nan')]]
+)
+def test_process_wavelength_bands_raises_if_bands_outside_bins(wavelength_bins, values):
+    bands = sc.array(dims=['wavelength'], values=values, unit='angstrom')
+    with pytest.raises(ValueError, match='must lie within'):
+        normalization.process_wavelength_bands(bands, wavelength_bins)
+
+
+def test_process_wavelength_bands_raises_if_band_is_reversed(wavelength_bins):
+    bands = sc.array(dims=['band', 'wavelength'], values=[[6.0, 3.0]], unit='angstrom')
+    with pytest.raises(ValueError, match='start before they end'):
+        normalization.process_wavelength_bands(bands, wavelength_bins)