From c17c64d1aadd2d0395f02780eaa325cd3abb231f Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 1 Feb 2024 08:40:28 +0100 Subject: [PATCH 01/11] Support event-mode I(Q) This uses an upper-bound error estimate. --- docs/examples/loki-direct-beam.ipynb | 4 +++- docs/examples/sans2d.ipynb | 31 ++++++++++++++++++++++------ docs/examples/zoom.ipynb | 3 ++- src/esssans/__init__.py | 2 ++ src/esssans/i_of_q.py | 7 ++++++- src/esssans/io.py | 7 +++++-- src/esssans/normalization.py | 16 +++++++++++--- src/esssans/types.py | 3 +++ src/esssans/uncertainty.py | 29 ++++++++++++++++++++++++++ 9 files changed, 88 insertions(+), 14 deletions(-) diff --git a/docs/examples/loki-direct-beam.ipynb b/docs/examples/loki-direct-beam.ipynb index e1ded591..fab51b3f 100644 --- a/docs/examples/loki-direct-beam.ipynb +++ b/docs/examples/loki-direct-beam.ipynb @@ -98,6 +98,7 @@ "\n", "params[CorrectForGravity] = True\n", "params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", + "params[sans.ReturnEvents] = False\n", "\n", "params[QBins] = sc.linspace(dim='Q', start=0.01, stop=0.3, num=101, unit='1/angstrom')" ] @@ -413,7 +414,8 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3" + "pygments_lexer": "ipython3", + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/examples/sans2d.ipynb b/docs/examples/sans2d.ipynb index 05334917..4d4e17a5 100644 --- a/docs/examples/sans2d.ipynb +++ b/docs/examples/sans2d.ipynb @@ -98,7 +98,8 @@ " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", ")\n", "pipeline[CorrectForGravity] = True\n", - "pipeline[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound" + "pipeline[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", + "pipeline[sans.ReturnEvents] = True" ] }, { @@ -158,7 +159,25 @@ "outputs": [], "source": [ "result = iofq.compute()\n", - "result.plot()" + "result.hist().plot(scale={'Q': 'log'}, norm='log')" + ] + }, + { + "cell_type": "markdown", + "id": "28532aa7", + "metadata": {}, + "source": [ + "As the result was computed in event-mode, we can also use a different $Q$-binning, without re-reducing the data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0748a1a", + "metadata": {}, + "outputs": [], + "source": [ + "result.hist(Q=60).plot(scale={'Q': 'log'}, norm='log')" ] }, { @@ -182,7 +201,7 @@ "result_drop = pipeline.compute(BackgroundSubtractedIofQ)\n", "# Reset the UnsertaintyBroadcastMode to the old value\n", "pipeline[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", - "sc.DataGroup(upper_bound=result, dropped=result_drop).plot(norm='log')" + "sc.DataGroup(upper_bound=result, dropped=result_drop).hist().plot(norm='log')" ] }, { @@ -291,7 +310,7 @@ ")\n", "pipeline.set_param_table(param_table)\n", "results = pipeline.compute(sciline.Series[Mode, BackgroundSubtractedIofQ])\n", - "sc.DataGroup(results).plot(norm='log')" + "sc.DataGroup(results).hist().plot(norm='log')" ] }, { @@ -357,7 +376,7 @@ "metadata": {}, "outputs": [], "source": [ - "pp.plot(sc.collapse(result, keep='Q'), norm='log')" + "pp.plot(sc.collapse(result.hist(), keep='Q'), norm='log')" ] } ], @@ -377,7 +396,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/examples/zoom.ipynb b/docs/examples/zoom.ipynb index 9e466c03..982a65fe 100644 --- a/docs/examples/zoom.ipynb +++ b/docs/examples/zoom.ipynb @@ -109,7 +109,8 @@ " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", ")\n", "params[CorrectForGravity] = True\n", - "params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound" + "params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", + "params[sans.ReturnEvents] = False" ] }, { diff --git a/src/esssans/__init__.py b/src/esssans/__init__.py index 87afcfb0..49021e69 100644 --- a/src/esssans/__init__.py +++ b/src/esssans/__init__.py @@ -21,6 +21,7 @@ ) from .common import transmission_from_background_run, transmission_from_sample_run from .direct_beam import direct_beam +from .types import ReturnEvents providers = ( *conversions.providers, @@ -41,6 +42,7 @@ del importlib __all__ = [ + 'ReturnEvents', 'beam_center_finder', 'common', 'conversions', diff --git a/src/esssans/i_of_q.py b/src/esssans/i_of_q.py index d7f4c817..6785ca8b 100644 --- a/src/esssans/i_of_q.py +++ b/src/esssans/i_of_q.py @@ -24,6 +24,7 @@ MonitorType, NonBackgroundWavelengthRange, QBins, + ReturnEvents, RunType, SampleRun, UncertaintyBroadcastMode, @@ -314,8 +315,12 @@ def _dense_merge_spectra( def subtract_background( - sample: IofQ[SampleRun], background: IofQ[BackgroundRun] + sample: IofQ[SampleRun], + background: IofQ[BackgroundRun], + return_events: ReturnEvents, ) -> BackgroundSubtractedIofQ: + if return_events and sample.bins is not None and background.bins is not None: + return sample.bins.concatenate(-background) if sample.bins is not None: sample = sample.bins.sum() if background.bins is not None: diff --git a/src/esssans/io.py b/src/esssans/io.py index 03954d1b..85a5323f 100644 --- a/src/esssans/io.py +++ b/src/esssans/io.py @@ -15,8 +15,11 @@ def save_background_subtracted_iofq( run_number: RunNumber, run_title: RunTitle, ) -> None: - """Save background subtracted IofQ as an NXcanSAS file.""" - da = iofq.copy(deep=False) + """Save background-subtracted I(Q) histogram as an NXcanSAS file.""" + if iofq.bins is None: + da = iofq.copy(deep=False) + else: + da = iofq.hist() if da.coords.is_edges('Q'): da.coords['Q'] = sc.midpoints(da.coords['Q']) with snx.File(out_filename, 'w') as f: diff --git a/src/esssans/normalization.py b/src/esssans/normalization.py index ce9682d2..b538a49b 100644 --- a/src/esssans/normalization.py +++ b/src/esssans/normalization.py @@ -20,6 +20,7 @@ LabFrameTransform, NormWavelengthTerm, Numerator, + ReturnEvents, RunType, SolidAngle, Transmission, @@ -28,6 +29,7 @@ UncertaintyBroadcastMode, ) from .uncertainty import ( + broadcast_to_events_with_upper_bound_variances, broadcast_with_upper_bound_variances, drop_variances_if_broadcast, ) @@ -298,6 +300,7 @@ def iofq_denominator( def normalize( numerator: CleanSummedQ[RunType, Numerator], denominator: CleanSummedQ[RunType, Denominator], + return_events: ReturnEvents, ) -> IofQ[RunType]: """ Perform normalization of counts as a function of Q. @@ -312,17 +315,24 @@ def normalize( denominator: The divisor for the normalization operation. This cannot be event data, it must contain histogrammed data. + return_events: + Whether to return the result as event data or histogrammed data. Returns ------- : The input data normalized by the supplied denominator. """ - if denominator.variances is not None and numerator.bins is not None: - # Event-mode normalization is not correct of norm-term has variances. + if return_events and numerator.bins is not None: + # Naive event-mode normalization is not correct if norm-term has variances. # See https://doi.org/10.3233/JNR-220049 for context. + if denominator.variances is not None: + denominator = broadcast_to_events_with_upper_bound_variances( + denominator, events=numerator + ) + else: numerator = numerator.hist() - if numerator.bins is not None: + if numerator.bins is not None and denominator.bins is None: da = numerator.bins / sc.lookup(func=denominator, dim='Q') else: da = numerator / denominator diff --git a/src/esssans/types.py b/src/esssans/types.py index 572a04b0..5da0e231 100644 --- a/src/esssans/types.py +++ b/src/esssans/types.py @@ -82,6 +82,9 @@ class TransmissionRun(sciline.Scope[RunType, int], int): See https://doi.org/10.3233/JNR-220049 for context. """ +ReturnEvents = NewType('ReturnEvents', bool) +"""Whether to return events in the output I(Q)""" + WavelengthBins = NewType('WavelengthBins', sc.Variable) """Wavelength binning""" diff --git a/src/esssans/uncertainty.py b/src/esssans/uncertainty.py index a0fa3018..f708fe86 100644 --- a/src/esssans/uncertainty.py +++ b/src/esssans/uncertainty.py @@ -51,3 +51,32 @@ def _no_variance_broadcast( return (data.variances is None) or all( data.sizes.get(dim) == size for dim, size in sizes.items() ) + + +def broadcast_to_events_with_upper_bound_variances( + da: sc.DataArray, *, events: sc.DataArray +) -> sc.DataArray: + """ + Upper-bound estimate for errors from normalization in event-mode. + + Count the number of events in each bin of the input data array. Then scale the + variances by the number of events in each bin. An explicit broadcast is performed + to bypass Scipp's safety check on broadcasting variances. + + Details will be published in an upcoming publication by Simon Heybrock et al. + """ + if da.variances is None: + return da + constituents = events.bins.constituents + Q = constituents['data'].coords['Q'] + constituents['data'] = sc.DataArray( + sc.ones(sizes=Q.sizes, dtype='float32'), coords={'Q': Q} + ) + # Combine all bins of the events that correspond to the same bin in the input data + to_concat = set(events.dims) - set(da.dims) + binned = sc.DataArray(sc.bins(**constituents).bins.concat(to_concat)) + counts = binned.hist(Q=da.coords['Q']) + da = da.copy() + da.variances *= counts.values + da.data = sc.bins_like(events, da.data) + return da From be912398fda2ab80a40bf0fb630a7f80bd822565 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 1 Feb 2024 08:44:45 +0100 Subject: [PATCH 02/11] Add "drop" broadcast mode. --- src/esssans/normalization.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/esssans/normalization.py b/src/esssans/normalization.py index b538a49b..294a893d 100644 --- a/src/esssans/normalization.py +++ b/src/esssans/normalization.py @@ -301,6 +301,7 @@ def normalize( numerator: CleanSummedQ[RunType, Numerator], denominator: CleanSummedQ[RunType, Denominator], return_events: ReturnEvents, + uncertainties: UncertaintyBroadcastMode, ) -> IofQ[RunType]: """ Perform normalization of counts as a function of Q. @@ -327,9 +328,12 @@ def normalize( # Naive event-mode normalization is not correct if norm-term has variances. # See https://doi.org/10.3233/JNR-220049 for context. if denominator.variances is not None: - denominator = broadcast_to_events_with_upper_bound_variances( - denominator, events=numerator - ) + if uncertainties == UncertaintyBroadcastMode.drop: + denominator = sc.values(denominator) + else: + denominator = broadcast_to_events_with_upper_bound_variances( + denominator, events=numerator + ) else: numerator = numerator.hist() if numerator.bins is not None and denominator.bins is None: From c0c105bc7dd1353013b3df3e65aa0fadf97c915a Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 1 Feb 2024 09:30:51 +0100 Subject: [PATCH 03/11] Add variances to event data loaded with ScippNexus --- src/esssans/loki/io.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/esssans/loki/io.py b/src/esssans/loki/io.py index b4395b66..07818090 100644 --- a/src/esssans/loki/io.py +++ b/src/esssans/loki/io.py @@ -37,6 +37,10 @@ def _patch_data( da: sc.DataArray, sample_position: sc.Variable, source_position: sc.Variable ) -> sc.DataArray: out = da.copy(deep=False) + if out.bins is not None: + content = out.bins.constituents['data'] + if content.variances is None: + content.variances = content.values out.coords['sample_position'] = sample_position out.coords['source_position'] = source_position out.coords['gravity'] = gravity_vector() From 54dfa06c79be9ee0b3d46cd2d497761efd491334 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 1 Feb 2024 09:39:39 +0100 Subject: [PATCH 04/11] Fix direct-beam iteration to work with data that has variances --- src/esssans/direct_beam.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/esssans/direct_beam.py b/src/esssans/direct_beam.py index 621e298d..b07f52e3 100644 --- a/src/esssans/direct_beam.py +++ b/src/esssans/direct_beam.py @@ -132,7 +132,7 @@ def direct_beam(pipeline: Pipeline, I0: sc.Variable, niter: int = 5) -> List[dic # parameters, nor given by any providers, so it will be considered flat. # TODO: Should we have a check that DirectBeam cannot be computed from the # pipeline? - iofq = pipeline.compute(BackgroundSubtractedIofQ) + iofq = sc.values(pipeline.compute(BackgroundSubtractedIofQ)) iofq_full = iofq['band', -1] iofq_bands = iofq['band', :-1] From 2f82e335b4092deed9df375481852d68481348ee Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 1 Feb 2024 09:40:28 +0100 Subject: [PATCH 05/11] Add and update tests --- src/esssans/__init__.py | 2 +- src/esssans/beam_center_finder.py | 2 ++ tests/loki/common.py | 2 ++ tests/loki/iofq_test.py | 38 +++++++++++++++++++++++++++++++ tests/sans2d/reduction_test.py | 4 +++- 5 files changed, 46 insertions(+), 2 deletions(-) diff --git a/src/esssans/__init__.py b/src/esssans/__init__.py index 49021e69..4ca9387e 100644 --- a/src/esssans/__init__.py +++ b/src/esssans/__init__.py @@ -21,7 +21,7 @@ ) from .common import transmission_from_background_run, transmission_from_sample_run from .direct_beam import direct_beam -from .types import ReturnEvents +from .types import BackgroundSubtractedIofQ, IofQ, ReturnEvents, SampleRun providers = ( *conversions.providers, diff --git a/src/esssans/beam_center_finder.py b/src/esssans/beam_center_finder.py index e6d482d5..1ce338e5 100644 --- a/src/esssans/beam_center_finder.py +++ b/src/esssans/beam_center_finder.py @@ -27,6 +27,7 @@ MaskedData, NormWavelengthTerm, QBins, + ReturnEvents, SampleRun, UncertaintyBroadcastMode, WavelengthBins, @@ -176,6 +177,7 @@ def _iofq_in_quadrants( ] params = {} params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound + params[ReturnEvents] = False params[WavelengthBins] = wavelength_bins params[QBins] = q_bins params[DetectorPixelShape[SampleRun]] = pixel_shape diff --git a/tests/loki/common.py b/tests/loki/common.py index a4d4a371..0fa0f54d 100644 --- a/tests/loki/common.py +++ b/tests/loki/common.py @@ -13,6 +13,7 @@ EmptyBeamRun, FileList, QBins, + ReturnEvents, SampleRun, TransmissionRun, UncertaintyBroadcastMode, @@ -65,6 +66,7 @@ def make_params( params[BeamStopRadius] = sc.scalar(0.042, unit='m') params[CorrectForGravity] = True params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound + params[ReturnEvents] = False params[QBins] = sc.linspace( dim='Q', start=0.01, stop=0.3, num=101, unit='1/angstrom' diff --git a/tests/loki/iofq_test.py b/tests/loki/iofq_test.py index 4cc8675f..395827a2 100644 --- a/tests/loki/iofq_test.py +++ b/tests/loki/iofq_test.py @@ -13,6 +13,7 @@ BackgroundSubtractedIofQ, BeamCenter, DimsToKeep, + ReturnEvents, UncertaintyBroadcastMode, WavelengthBands, ) @@ -39,6 +40,43 @@ def test_pipeline_can_compute_IofQ(uncertainties): assert result.dims == ('Q',) +@pytest.mark.parametrize( + 'uncertainties', + [UncertaintyBroadcastMode.drop, UncertaintyBroadcastMode.upper_bound], +) +@pytest.mark.parametrize( + 'target', [sans.IofQ[sans.SampleRun], sans.BackgroundSubtractedIofQ] +) +def test_pipeline_can_compute_IofQ_in_event_mode(uncertainties, target): + params = make_params() + params[UncertaintyBroadcastMode] = uncertainties + pipeline = sciline.Pipeline(loki_providers(), params=params) + reference = pipeline.compute(target) + pipeline[ReturnEvents] = True + result = pipeline.compute(target) + assert result.bins is not None + assert sc.allclose( + sc.values(reference.data), + sc.values(result.hist().data), + rtol=sc.scalar(1e-11), + atol=sc.scalar(1e-11), + ) + if uncertainties == UncertaintyBroadcastMode.drop: + assert sc.allclose( + sc.variances(reference).data, + sc.variances(result.hist()).data, + rtol=sc.scalar(1e-14), + atol=sc.scalar(1e-14), + ) + else: + assert sc.allclose( + sc.variances(reference).data, + sc.variances(result.hist()).data, + rtol=sc.scalar(1e-9), + atol=sc.scalar(1e-9), + ) + + def test_pipeline_can_compute_IofQ_in_wavelength_slices(): params = make_params() band = np.linspace(1.0, 13.0, num=11) diff --git a/tests/sans2d/reduction_test.py b/tests/sans2d/reduction_test.py index b9e31210..dbff8c63 100644 --- a/tests/sans2d/reduction_test.py +++ b/tests/sans2d/reduction_test.py @@ -22,6 +22,7 @@ NonBackgroundWavelengthRange, QBins, RawData, + ReturnEvents, SampleRun, SolidAngle, TransmissionRun, @@ -61,6 +62,7 @@ def make_params() -> dict: ) params[CorrectForGravity] = True params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound + params[ReturnEvents] = False return params @@ -140,7 +142,7 @@ def test_uncertainty_broadcast_mode_drop_yields_smaller_variances(): ) params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop pipeline = sciline.Pipeline(sans2d_providers(), params=params) - drop = pipeline.compute(IofQ[SampleRun]).hist().data + drop = pipeline.compute(IofQ[SampleRun]).data params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound pipeline = sciline.Pipeline(sans2d_providers(), params=params) upper_bound = pipeline.compute(IofQ[SampleRun]).data From 0b4880add1f0ccb37e69d6ddcf12128639feb2f1 Mon Sep 17 00:00:00 2001 From: Simon Heybrock <12912489+SimonHeybrock@users.noreply.github.com> Date: Mon, 5 Feb 2024 06:53:26 +0100 Subject: [PATCH 06/11] Update tests/loki/iofq_test.py Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com> --- tests/loki/iofq_test.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/tests/loki/iofq_test.py b/tests/loki/iofq_test.py index 395827a2..15959518 100644 --- a/tests/loki/iofq_test.py +++ b/tests/loki/iofq_test.py @@ -62,19 +62,15 @@ def test_pipeline_can_compute_IofQ_in_event_mode(uncertainties, target): atol=sc.scalar(1e-11), ) if uncertainties == UncertaintyBroadcastMode.drop: - assert sc.allclose( - sc.variances(reference).data, - sc.variances(result.hist()).data, - rtol=sc.scalar(1e-14), - atol=sc.scalar(1e-14), - ) + tol = sc.scalar(1e-14) else: - assert sc.allclose( - sc.variances(reference).data, - sc.variances(result.hist()).data, - rtol=sc.scalar(1e-9), - atol=sc.scalar(1e-9), - ) + tol = sc.scalar(1e-9) + assert sc.allclose( + sc.variances(reference).data, + sc.variances(result.hist()).data, + rtol=tol, + atol=tol, + ) def test_pipeline_can_compute_IofQ_in_wavelength_slices(): From ca4eb37e1c696626d39a736b065ba0ce634789fe Mon Sep 17 00:00:00 2001 From: Simon Heybrock <12912489+SimonHeybrock@users.noreply.github.com> Date: Mon, 5 Feb 2024 08:34:03 +0100 Subject: [PATCH 07/11] Update src/esssans/__init__.py --- src/esssans/__init__.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/esssans/__init__.py b/src/esssans/__init__.py index 4ca9387e..18f929fe 100644 --- a/src/esssans/__init__.py +++ b/src/esssans/__init__.py @@ -42,6 +42,9 @@ del importlib __all__ = [ + 'BackgroundSubtractedIofQ', + 'IofQ', + 'SampleRun', 'ReturnEvents', 'beam_center_finder', 'common', From cb97c6e6fc8a259381ae81b9128d18b844a04639 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 5 Feb 2024 13:32:33 +0100 Subject: [PATCH 08/11] Refactor to handle wavelength bands during normalization step --- docs/examples/sans2d.ipynb | 11 ++++- src/esssans/conversions.py | 5 +- src/esssans/i_of_q.py | 93 ++++++------------------------------ src/esssans/normalization.py | 32 +++++++++++++ 4 files changed, 60 insertions(+), 81 deletions(-) diff --git a/docs/examples/sans2d.ipynb b/docs/examples/sans2d.ipynb index a0f6d011..9fdb8ce1 100644 --- a/docs/examples/sans2d.ipynb +++ b/docs/examples/sans2d.ipynb @@ -265,8 +265,15 @@ "\n", "display(plot_flat_detector_xy(results[MaskedData[SampleRun]].sum('tof'), norm='log'))\n", "\n", - "parts = {str(key): results[key] for key in parts}\n", - "parts = {key: val if val.bins is None else val.hist() for key, val in parts.items()}\n", + "wavelength = pipeline.compute(WavelengthBins)\n", + "display(\n", + " results[CleanSummedQ[SampleRun, Numerator]]\n", + " .hist(wavelength=wavelength)\n", + " .transpose()\n", + " .plot(norm='log')\n", + ")\n", + "display(results[CleanSummedQ[SampleRun, Denominator]].plot(norm='log'))\n", + "parts = {str(key): results[key].sum('wavelength') for key in parts}\n", "display(sc.plot(parts, norm='log'))\n", "\n", "iofqs = {str(key): results[key] for key in iofqs}\n", diff --git a/src/esssans/conversions.py b/src/esssans/conversions.py index 1e30ecb8..b30f7648 100644 --- a/src/esssans/conversions.py +++ b/src/esssans/conversions.py @@ -228,7 +228,10 @@ def to_Q( """ Convert a data array from wavelength to Q. """ - return CleanQ[RunType, IofQPart](data.transform_coords('Q', graph=graph)) + # Keep naming of wavelength dim, subsequent steps use a (Q, wavelength) binning. + return CleanQ[RunType, IofQPart]( + data.transform_coords('Q', graph=graph, rename_dims=False) + ) providers = ( diff --git a/src/esssans/i_of_q.py b/src/esssans/i_of_q.py index 0782cddf..147a6f5b 100644 --- a/src/esssans/i_of_q.py +++ b/src/esssans/i_of_q.py @@ -2,10 +2,8 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) from typing import Dict, List, Optional, Union -from uuid import uuid4 import scipp as sc -from scipp.core.concepts import irreducible_mask from scipp.scipy.interpolate import interp1d from .common import mask_range @@ -170,11 +168,7 @@ def process_wavelength_bands( def merge_spectra( - data: CleanQ[RunType, IofQPart], - q_bins: QBins, - wavelength_bins: WavelengthBins, - wavelength_bands: ProcessedWavelengthBands, - dims_to_keep: Optional[DimsToKeep], + data: CleanQ[RunType, IofQPart], q_bins: QBins, dims_to_keep: Optional[DimsToKeep] ) -> CleanSummedQ[RunType, IofQPart]: """ Merges all spectra: @@ -188,13 +182,6 @@ def merge_spectra( A DataArray containing the data that is to be converted to Q. q_bins: The binning in Q to be used. - wavelength_bins: - The binning in wavelength to be used. - wavelength_bands: - Defines bands in wavelength that can be used to separate different wavelength - ranges that contribute to different regions in Q space. Note that this needs to - be defined, so if all wavelengths should be used, this should simply be a start - and end edges that encompass the entire wavelength range. dims_to_keep: Dimensions that should not be reduced and thus still be present in the final I(Q) result (this is typically the layer dimension). @@ -204,23 +191,17 @@ def merge_spectra( : The input data converted to Q and then summed over all detector pixels. """ - dims_to_reduce = set(data.dims) - {'Q'} + dims_to_reduce = set(data.dims) - {'wavelength'} if dims_to_keep is not None: dims_to_reduce -= set(dims_to_keep) if data.bins is not None: out = _events_merge_spectra( - data_q=data, - q_bins=q_bins, - wavelength_bands=wavelength_bands, - dims_to_reduce=dims_to_reduce, + data_q=data, q_bins=q_bins, dims_to_reduce=dims_to_reduce ) else: out = _dense_merge_spectra( - data_q=data, - q_bins=q_bins, - wavelength_bands=wavelength_bands, - dims_to_reduce=dims_to_reduce, + data_q=data, q_bins=q_bins, dims_to_reduce=dims_to_reduce ) return CleanSummedQ[RunType, IofQPart](out.squeeze()) @@ -236,87 +217,43 @@ def _to_q_bins(q_bins: Union[int, sc.Variable]) -> Dict[str, Union[int, sc.Varia def _events_merge_spectra( - data_q: sc.DataArray, - q_bins: Union[int, sc.Variable], - dims_to_reduce: List[str], - wavelength_bands: sc.Variable, + data_q: sc.DataArray, q_bins: Union[int, sc.Variable], dims_to_reduce: List[str] ) -> sc.DataArray: """ Merge spectra of event data """ - q_all_pixels = data_q.bins.concat(dims_to_reduce) edges = _to_q_bins(q_bins) - q_binned = q_all_pixels.bin(**edges) - dim = 'wavelength' - wav_binned = q_binned.bin({dim: sc.sort(wavelength_bands.flatten(to=dim), dim)}) - # At this point we kind of already have what we need, would be cheapest to just - # return, if follow up providers can work with the result. - # Otherwise we need to duplicate events: - sections = [] - for bounds in sc.collapse(wavelength_bands, keep=dim).values(): - # The extra concat can probably be avoided if we insert some dummy edges for - # first and last band, but we would need to know how many edges to insert, as - # the bands can be very wide and overlap by more than one bin. - sections.append(wav_binned[dim, bounds[0] : bounds[1]].bins.concat(dim)) - band_dim = (set(wavelength_bands.dims) - {'wavelength'}).pop() - out = sc.concat(sections, band_dim) - out.coords[dim] = wavelength_bands - return out + q_all_pixels = data_q.bins.concat(dims_to_reduce) + return q_all_pixels.bin(**edges) def _dense_merge_spectra( - data_q: sc.DataArray, - q_bins: Union[int, sc.Variable], - dims_to_reduce: List[str], - wavelength_bands: sc.Variable, + data_q: sc.DataArray, q_bins: Union[int, sc.Variable], dims_to_reduce: set[str] ) -> sc.DataArray: """ Merge spectra of dense data """ edges = _to_q_bins(q_bins) - bands = [] - band_dim = (set(wavelength_bands.dims) - {'wavelength'}).pop() # We want to flatten data to make histogramming cheaper (avoiding allocation of # large output before summing). We strip unnecessary content since it makes # flattening more expensive. stripped = data_q.copy(deep=False) for name, coord in data_q.coords.items(): - if name not in ['Q', 'wavelength'] and any( - [dim in dims_to_reduce for dim in coord.dims] - ): + if name not in {'Q', 'wavelength'} and set(coord.dims) & dims_to_reduce: del stripped.coords[name] to_flatten = [dim for dim in data_q.dims if dim in dims_to_reduce] - dummy_dim = str(uuid4()) - # Make sure that dims to flatten are contiguous, and that Q is the last dim + # Make dims to flatten contiguous, keep wavelength is the last dim data_dims = list(stripped.dims) - rearranged_dims = to_flatten + ['Q'] - for dim in rearranged_dims: + for dim in to_flatten + ['wavelength']: data_dims.remove(dim) - for dim in rearranged_dims: data_dims.append(dim) stripped = stripped.transpose(data_dims) - flat = stripped.flatten(dims=to_flatten, to=dummy_dim) - - # Apply masks once, to avoid repeated work when iterating over bands - mask = irreducible_mask(flat, dummy_dim) - # When not all dims are reduced there may be extra dims in the mask and it is not - # possible to select data based on it. In this case the masks will be applied - # in the loop below, which is slightly slower. - if mask.ndim == 1: - flat = flat.drop_masks( - [name for name, mask in flat.masks.items() if dummy_dim in mask.dims] - ) - flat = flat[~mask] - - dims_to_reduce = tuple(dim for dim in dims_to_reduce if dim not in to_flatten) - for wav_range in sc.collapse(wavelength_bands, keep='wavelength').values(): - band = flat['wavelength', wav_range[0] : wav_range[1]] - # By flattening before histogramming we avoid allocating a large output array, - # which would then require summing over all pixels. - bands.append(band.flatten(dims=(dummy_dim, 'Q'), to='Q').hist(**edges)) - return sc.concat(bands, band_dim) + # Flatten to Q such that `hist` below will turn this into the new Q dim. + flat = stripped.flatten(dims=to_flatten, to='Q') + + return flat.hist(**edges) def subtract_background( diff --git a/src/esssans/normalization.py b/src/esssans/normalization.py index 294a893d..bc45266b 100644 --- a/src/esssans/normalization.py +++ b/src/esssans/normalization.py @@ -20,6 +20,7 @@ LabFrameTransform, NormWavelengthTerm, Numerator, + ProcessedWavelengthBands, ReturnEvents, RunType, SolidAngle, @@ -302,6 +303,7 @@ def normalize( denominator: CleanSummedQ[RunType, Denominator], return_events: ReturnEvents, uncertainties: UncertaintyBroadcastMode, + wavelength_bands: ProcessedWavelengthBands, ) -> IofQ[RunType]: """ Perform normalization of counts as a function of Q. @@ -318,12 +320,42 @@ def normalize( contain histogrammed data. return_events: Whether to return the result as event data or histogrammed data. + wavelength_bands: + Defines bands in wavelength that can be used to separate different wavelength + ranges that contribute to different regions in Q space. Note that this needs to + be defined, so if all wavelengths should be used, this should simply be a start + and end edges that encompass the entire wavelength range. Returns ------- : The input data normalized by the supplied denominator. """ + wav = 'wavelength' + wavelength_bounds = sc.sort(wavelength_bands.flatten(to=wav), wav) + if numerator.bins is not None: + # 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. + numerator = numerator.bin(wavelength=wavelength_bounds) + + def _reduce(da: sc.DataArray) -> sc.DataArray: + return da.sum(wav) if da.bins is None else da.bins.concat(wav) + + num_parts = [] + denom_parts = [] + for wav_range in sc.collapse(wavelength_bands, keep=wav).values(): + num_parts.append(_reduce(numerator[wav, wav_range[0] : wav_range[1]])) + denom_parts.append(_reduce(denominator[wav, wav_range[0] : wav_range[1]])) + band_dim = (set(wavelength_bands.dims) - {'wavelength'}).pop() + if len(num_parts) == 1: + numerator = num_parts[0] + denominator = denom_parts[0] + else: + numerator = sc.concat(num_parts, band_dim) + denominator = sc.concat(denom_parts, band_dim) + numerator.coords[wav] = wavelength_bands.squeeze() + denominator.coords[wav] = wavelength_bands.squeeze() + if return_events and numerator.bins is not None: # Naive event-mode normalization is not correct if norm-term has variances. # See https://doi.org/10.3233/JNR-220049 for context. From 50122c343dfc1053ba5b1b0c68cbb3441c13024f Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 5 Feb 2024 13:53:39 +0100 Subject: [PATCH 09/11] Update notebook --- docs/examples/zoom.ipynb | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/docs/examples/zoom.ipynb b/docs/examples/zoom.ipynb index 9d102a9b..e77946ee 100644 --- a/docs/examples/zoom.ipynb +++ b/docs/examples/zoom.ipynb @@ -247,8 +247,19 @@ " )\n", ")\n", "\n", + "wavelength = pipeline.compute(WavelengthBins)\n", + "display(\n", + " results[CleanSummedQ[SampleRun, Numerator]]\n", + " .hist(wavelength=wavelength)\n", + " .transpose()\n", + " .plot(norm='log')\n", + ")\n", + "display(results[CleanSummedQ[SampleRun, Denominator]].plot(norm='log'))\n", "parts = {str(key): results[key] for key in parts}\n", - "parts = {key: val if val.bins is None else val.hist() for key, val in parts.items()}\n", + "parts = {\n", + " key: val.sum('wavelength') if val.bins is None else val.hist()\n", + " for key, val in parts.items()\n", + "}\n", "display(sc.plot(parts, norm='log', scale={'Q': 'log'}))\n", "\n", "iofqs = {str(key): results[key] for key in iofqs}\n", From 6268fbed96e5aecb419da3f880f56c84a438c3e1 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 6 Feb 2024 05:26:01 +0100 Subject: [PATCH 10/11] Inline helper functions --- src/esssans/i_of_q.py | 81 +++++++++++++------------------------------ 1 file changed, 24 insertions(+), 57 deletions(-) diff --git a/src/esssans/i_of_q.py b/src/esssans/i_of_q.py index 147a6f5b..476e18ea 100644 --- a/src/esssans/i_of_q.py +++ b/src/esssans/i_of_q.py @@ -1,7 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -from typing import Dict, List, Optional, Union +from typing import Optional import scipp as sc from scipp.scipy.interpolate import interp1d @@ -195,67 +195,34 @@ def merge_spectra( if dims_to_keep is not None: dims_to_reduce -= set(dims_to_keep) + edges = {'Q': q_bins} if isinstance(q_bins, int) else {q_bins.dim: q_bins} + if data.bins is not None: - out = _events_merge_spectra( - data_q=data, q_bins=q_bins, dims_to_reduce=dims_to_reduce - ) + q_all_pixels = data.bins.concat(dims_to_reduce) + out = q_all_pixels.bin(**edges) else: - out = _dense_merge_spectra( - data_q=data, q_bins=q_bins, dims_to_reduce=dims_to_reduce - ) + # We want to flatten data to make histogramming cheaper (avoiding allocation of + # large output before summing). We strip unnecessary content since it makes + # flattening more expensive. + stripped = data.copy(deep=False) + for name, coord in data.coords.items(): + if name not in {'Q', 'wavelength'} and set(coord.dims) & dims_to_reduce: + del stripped.coords[name] + to_flatten = [dim for dim in data.dims if dim in dims_to_reduce] + + # Make dims to flatten contiguous, keep wavelength as the last dim + data_dims = list(stripped.dims) + for dim in to_flatten + ['wavelength']: + data_dims.remove(dim) + data_dims.append(dim) + stripped = stripped.transpose(data_dims) + # Flatten to Q such that `hist` below will turn this into the new Q dim. + flat = stripped.flatten(dims=to_flatten, to='Q') + + out = flat.hist(**edges) return CleanSummedQ[RunType, IofQPart](out.squeeze()) -def _to_q_bins(q_bins: Union[int, sc.Variable]) -> Dict[str, Union[int, sc.Variable]]: - """ - If the input bins are an integer, convert them to a dictionary that can be used - to bin a DataArray. - """ - if isinstance(q_bins, int): - return {'Q': q_bins} - return {q_bins.dim: q_bins} - - -def _events_merge_spectra( - data_q: sc.DataArray, q_bins: Union[int, sc.Variable], dims_to_reduce: List[str] -) -> sc.DataArray: - """ - Merge spectra of event data - """ - edges = _to_q_bins(q_bins) - q_all_pixels = data_q.bins.concat(dims_to_reduce) - return q_all_pixels.bin(**edges) - - -def _dense_merge_spectra( - data_q: sc.DataArray, q_bins: Union[int, sc.Variable], dims_to_reduce: set[str] -) -> sc.DataArray: - """ - Merge spectra of dense data - """ - edges = _to_q_bins(q_bins) - - # We want to flatten data to make histogramming cheaper (avoiding allocation of - # large output before summing). We strip unnecessary content since it makes - # flattening more expensive. - stripped = data_q.copy(deep=False) - for name, coord in data_q.coords.items(): - if name not in {'Q', 'wavelength'} and set(coord.dims) & dims_to_reduce: - del stripped.coords[name] - to_flatten = [dim for dim in data_q.dims if dim in dims_to_reduce] - - # Make dims to flatten contiguous, keep wavelength is the last dim - data_dims = list(stripped.dims) - for dim in to_flatten + ['wavelength']: - data_dims.remove(dim) - data_dims.append(dim) - stripped = stripped.transpose(data_dims) - # Flatten to Q such that `hist` below will turn this into the new Q dim. - flat = stripped.flatten(dims=to_flatten, to='Q') - - return flat.hist(**edges) - - def subtract_background( sample: IofQ[SampleRun], background: IofQ[BackgroundRun], From 2a728d73973e69348033476afbf2373c4a3dc2af Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 6 Feb 2024 05:29:54 +0100 Subject: [PATCH 11/11] Move function to related file --- src/esssans/beam_center_finder.py | 9 +++++-- src/esssans/i_of_q.py | 39 ------------------------------- src/esssans/normalization.py | 39 +++++++++++++++++++++++++++++++ 3 files changed, 46 insertions(+), 41 deletions(-) diff --git a/src/esssans/beam_center_finder.py b/src/esssans/beam_center_finder.py index c7b8cb6d..0fba6f6a 100644 --- a/src/esssans/beam_center_finder.py +++ b/src/esssans/beam_center_finder.py @@ -16,9 +16,14 @@ mask_wavelength, to_Q, ) -from .i_of_q import merge_spectra, process_wavelength_bands +from .i_of_q import merge_spectra from .logging import get_logger -from .normalization import iofq_denominator, normalize, solid_angle +from .normalization import ( + iofq_denominator, + normalize, + process_wavelength_bands, + solid_angle, +) from .types import ( BeamCenter, DetectorPixelShape, diff --git a/src/esssans/i_of_q.py b/src/esssans/i_of_q.py index 476e18ea..79938f4d 100644 --- a/src/esssans/i_of_q.py +++ b/src/esssans/i_of_q.py @@ -21,13 +21,11 @@ IofQPart, MonitorType, NonBackgroundWavelengthRange, - ProcessedWavelengthBands, QBins, ReturnEvents, RunType, SampleRun, UncertaintyBroadcastMode, - WavelengthBands, WavelengthBins, WavelengthMonitor, ) @@ -131,42 +129,6 @@ def resample_direct_beam( return CleanDirectBeam(func(wavelength_bins, midpoints=True)) -def process_wavelength_bands( - wavelength_bands: Optional[WavelengthBands], - wavelength_bins: WavelengthBins, -) -> ProcessedWavelengthBands: - """ - Perform some checks and potential reshaping on the wavelength bands. - - The wavelength bands must be either one- or two-dimensional. - If the wavelength bands are defined as a one-dimensional array, convert them to a - two-dimensional array with start and end wavelengths. - - The final bands must have a size of 2 in the wavelength dimension, defining a start - and an end wavelength. - """ - if wavelength_bands is None: - wavelength_bands = sc.concat( - [wavelength_bins.min(), wavelength_bins.max()], dim='wavelength' - ) - if wavelength_bands.ndim == 1: - wavelength_bands = sc.concat( - [wavelength_bands[:-1], wavelength_bands[1:]], dim='x' - ).rename(x='wavelength', wavelength='band') - if wavelength_bands.ndim != 2: - raise ValueError( - 'Wavelength_bands must be one- or two-dimensional, ' - f'got {wavelength_bands.ndim}.' - ) - if wavelength_bands.sizes['wavelength'] != 2: - raise ValueError( - 'Wavelength_bands must have a size of 2 in the wavelength dimension, ' - 'defining a start and an end wavelength, ' - f'got {wavelength_bands.sizes["wavelength"]}.' - ) - return wavelength_bands - - def merge_spectra( data: CleanQ[RunType, IofQPart], q_bins: QBins, dims_to_keep: Optional[DimsToKeep] ) -> CleanSummedQ[RunType, IofQPart]: @@ -242,5 +204,4 @@ def subtract_background( resample_direct_beam, merge_spectra, subtract_background, - process_wavelength_bands, ) diff --git a/src/esssans/normalization.py b/src/esssans/normalization.py index bc45266b..a1a4a063 100644 --- a/src/esssans/normalization.py +++ b/src/esssans/normalization.py @@ -28,6 +28,8 @@ TransmissionFraction, TransmissionRun, UncertaintyBroadcastMode, + WavelengthBands, + WavelengthBins, ) from .uncertainty import ( broadcast_to_events_with_upper_bound_variances, @@ -298,6 +300,42 @@ def iofq_denominator( return CleanWavelength[RunType, Denominator](denominator) +def process_wavelength_bands( + wavelength_bands: Optional[WavelengthBands], + wavelength_bins: WavelengthBins, +) -> ProcessedWavelengthBands: + """ + Perform some checks and potential reshaping on the wavelength bands. + + The wavelength bands must be either one- or two-dimensional. + If the wavelength bands are defined as a one-dimensional array, convert them to a + two-dimensional array with start and end wavelengths. + + The final bands must have a size of 2 in the wavelength dimension, defining a start + and an end wavelength. + """ + if wavelength_bands is None: + wavelength_bands = sc.concat( + [wavelength_bins.min(), wavelength_bins.max()], dim='wavelength' + ) + if wavelength_bands.ndim == 1: + wavelength_bands = sc.concat( + [wavelength_bands[:-1], wavelength_bands[1:]], dim='x' + ).rename(x='wavelength', wavelength='band') + if wavelength_bands.ndim != 2: + raise ValueError( + 'Wavelength_bands must be one- or two-dimensional, ' + f'got {wavelength_bands.ndim}.' + ) + if wavelength_bands.sizes['wavelength'] != 2: + raise ValueError( + 'Wavelength_bands must have a size of 2 in the wavelength dimension, ' + 'defining a start and an end wavelength, ' + f'got {wavelength_bands.sizes["wavelength"]}.' + ) + return wavelength_bands + + def normalize( numerator: CleanSummedQ[RunType, Numerator], denominator: CleanSummedQ[RunType, Denominator], @@ -380,5 +418,6 @@ def _reduce(da: sc.DataArray) -> sc.DataArray: iofq_norm_wavelength_term, iofq_denominator, normalize, + process_wavelength_bands, solid_angle, )