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)