From e18d5eec7d2adcae783c1d6d881af7bedafdfc7d Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 25 Aug 2026 07:34:17 +0000 Subject: [PATCH] fix(sans): snap wavelength band edges onto WavelengthBins Band boundaries are realized by slicing the wavelength axis, which can only cut at bin boundaries. Event data was re-binned onto the requested band edges and was therefore exact, while the dense denominator was sliced and snapped to whole bins. Numerator and denominator then covered different wavelength ranges, biasing I(Q) per band by several percent for bands only a few bins wide. The dense parts are inconsistent among themselves too: the denominator carries wavelength midpoints and so assigns each bin to the band containing its midpoint, whereas the numerator that `direct_beam` histograms carries bin edges and so picks up every overlapping bin. On real Larmor data the existing overlapping-bands direct-beam configuration was off by up to 22%. Snapping the band edges onto WavelengthBins makes all three representations select identical wavelength ranges. Aligned bands, which is everything in the docs and tests bar the overlapping-bands case, are unchanged bit-for-bit. Bands falling outside WavelengthBins used to be clamped to whatever range existed; they now raise, as do reversed bands and NaN edges, so that the remaining silent adjustment is bounded by half a bin width. --- .../esssans/src/ess/sans/normalization.py | 64 +++++++- packages/esssans/src/ess/sans/types.py | 8 +- packages/esssans/tests/normalization_test.py | 140 ++++++++++++++++++ 3 files changed, 210 insertions(+), 2 deletions(-) diff --git a/packages/esssans/src/ess/sans/normalization.py b/packages/esssans/src/ess/sans/normalization.py index cc32bdcbe..3720c6bd6 100644 --- a/packages/esssans/src/ess/sans/normalization.py +++ b/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() diff --git a/packages/esssans/src/ess/sans/types.py b/packages/esssans/src/ess/sans/types.py index 1c7152070..fd824dfd2 100644 --- a/packages/esssans/src/ess/sans/types.py +++ b/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) diff --git a/packages/esssans/tests/normalization_test.py b/packages/esssans/tests/normalization_test.py index 4ea4e7a57..0608b9bde 100644 --- a/packages/esssans/tests/normalization_test.py +++ b/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)