From c7df6c26bd7cec55853ccfe8d7f26c7c48640d4c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 7 Oct 2024 11:38:02 +0200 Subject: [PATCH 1/5] fix: include resolution in stitched curve if present --- src/ess/reflectometry/tools.py | 27 +++++++++++++++++++++++++-- tests/tools_test.py | 23 +++++++++++++++++++++++ 2 files changed, 48 insertions(+), 2 deletions(-) diff --git a/src/ess/reflectometry/tools.py b/src/ess/reflectometry/tools.py index 3409789d..cbe36721 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -229,7 +229,7 @@ def cost(scaling_factors): def combine_curves( curves: Sequence[sc.DataArray], - q_bin_edges: sc.Variable | None = None, + q_bin_edges: sc.Variable, ) -> sc.DataArray: '''Combines the given curves by interpolating them on a 1d grid defined by :code:`q_bin_edges` and averaging @@ -268,7 +268,8 @@ def combine_curves( inv_v = 1.0 / v r_avg = np.nansum(r * inv_v, axis=0) / np.nansum(inv_v, axis=0) v_avg = 1 / np.nansum(inv_v, axis=0) - return sc.DataArray( + + out = sc.DataArray( data=sc.array( dims='Q', values=r_avg, @@ -277,6 +278,28 @@ def combine_curves( ), coords={'Q': q_bin_edges}, ) + if any('Q_resolution' in c.coords for c in curves): + # This might need to be revisited. The question about how to combine curves + # with different Q-resolution is not completely resolved. + # However, in practice the difference in Q-resolution between different curves + # is small so it's not likely to make a big difference. + q_res = ( + sc.DataArray( + data=c.coords.get( + 'Q_resolution', sc.scalar(float('nan')) * sc.values(c.data.copy()) + ), + coords={'Q': c.coords['Q']}, + ) + for c in curves + ) + qs = _interpolate_on_qgrid(q_res, q_bin_edges).values + qs_avg = np.nansum(qs * inv_v, axis=0) / np.nansum( + ~np.isnan(qs) * inv_v, axis=0 + ) + out.coords['Q_resolution'] = sc.array( + dims='Q', values=qs_avg, unit=next(iter(curves)).coords['Q_resolution'].unit + ) + return out def orso_datasets_from_measurements( diff --git a/tests/tools_test.py b/tests/tools_test.py index 03c07aaf..345d233c 100644 --- a/tests/tools_test.py +++ b/tests/tools_test.py @@ -146,6 +146,29 @@ def test_combined_curves(): ) +@pytest.mark.filterwarnings("ignore:invalid value encountered in divide") +def test_combined_curves_resolution(): + qgrid = sc.linspace('Q', 0, 1, 26) + data = sc.concat( + ( + sc.ones(dims=['Q'], shape=[10], with_variances=True), + 0.5 * sc.ones(dims=['Q'], shape=[15], with_variances=True), + ), + dim='Q', + ) + data.variances[:] = 0.1 + curves = ( + curve(data, 0, 0.3), + curve(0.5 * data, 0.2, 0.7), + curve(0.25 * data, 0.6, 1.0), + ) + curves[0].coords['Q_resolution'] = sc.midpoints(curves[0].coords['Q']) / 5 + combined = combine_curves(curves, qgrid) + assert 'Q_resolution' in combined.coords + assert combined.coords['Q_resolution'][0] == curves[0].coords['Q_resolution'][1] + assert sc.isnan(combined.coords['Q_resolution'][-1]) + + def test_linlogspace_linear(): q_lin = linlogspace( dim='qz', edges=[0.008, 0.08], scale='linear', num=50, unit='1/angstrom' From 56852e6c87ccd76a13ce0b88e9966ac9df962ba3 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 10 Oct 2024 15:01:43 +0200 Subject: [PATCH 2/5] fix: use scipp nansum --- src/ess/reflectometry/tools.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/src/ess/reflectometry/tools.py b/src/ess/reflectometry/tools.py index cbe36721..529e4060 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -261,19 +261,19 @@ def combine_curves( if len({c.coords['Q'].unit for c in curves}) != 1: raise ValueError('The Q-coordinates must have the same unit for each curve') - r = _interpolate_on_qgrid(map(sc.values, curves), q_bin_edges).values - v = _interpolate_on_qgrid(map(sc.variances, curves), q_bin_edges).values + r = _interpolate_on_qgrid(map(sc.values, curves), q_bin_edges) + v = _interpolate_on_qgrid(map(sc.variances, curves), q_bin_edges) - v[v == 0] = np.nan + v = sc.where(v == 0, sc.scalar(np.nan, unit=v.unit), v) inv_v = 1.0 / v - r_avg = np.nansum(r * inv_v, axis=0) / np.nansum(inv_v, axis=0) - v_avg = 1 / np.nansum(inv_v, axis=0) + r_avg = sc.nansum(r * inv_v, dim='curves') / sc.nansum(inv_v, dim='curves') + v_avg = 1 / sc.nansum(inv_v, dim='curves') out = sc.DataArray( data=sc.array( dims='Q', - values=r_avg, - variances=v_avg, + values=r_avg.values, + variances=v_avg.values, unit=next(iter(curves)).data.unit, ), coords={'Q': q_bin_edges}, @@ -286,18 +286,15 @@ def combine_curves( q_res = ( sc.DataArray( data=c.coords.get( - 'Q_resolution', sc.scalar(float('nan')) * sc.values(c.data.copy()) + 'Q_resolution', sc.full_like(c.coords['Q'], value=np.nan) ), coords={'Q': c.coords['Q']}, ) for c in curves ) - qs = _interpolate_on_qgrid(q_res, q_bin_edges).values - qs_avg = np.nansum(qs * inv_v, axis=0) / np.nansum( - ~np.isnan(qs) * inv_v, axis=0 - ) - out.coords['Q_resolution'] = sc.array( - dims='Q', values=qs_avg, unit=next(iter(curves)).coords['Q_resolution'].unit + qs = _interpolate_on_qgrid(q_res, q_bin_edges) + out.coords['Q_resolution'] = sc.nansum(qs * inv_v, dim='curves') / sc.nansum( + sc.where(sc.isnan(qs), sc.scalar(0.0, unit=inv_v.unit), inv_v), dim='curves' ) return out From ec567785c0843da25bbc4b4f5988e0a9a1a07e71 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 18 Dec 2025 09:37:43 +0100 Subject: [PATCH 3/5] fix: remove moved test --- tests/tools_test.py | 272 -------------------------------------------- 1 file changed, 272 deletions(-) delete mode 100644 tests/tools_test.py diff --git a/tests/tools_test.py b/tests/tools_test.py deleted file mode 100644 index 345d233c..00000000 --- a/tests/tools_test.py +++ /dev/null @@ -1,272 +0,0 @@ -# SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -import numpy as np -import pytest -import sciline as sl -import scipp as sc -from numpy.testing import assert_allclose as np_assert_allclose -from orsopy.fileio import Orso, OrsoDataset -from scipp.testing import assert_allclose - -from ess.reflectometry.orso import OrsoIofQDataset -from ess.reflectometry.tools import ( - combine_curves, - linlogspace, - orso_datasets_from_measurements, - scale_reflectivity_curves_to_overlap, -) -from ess.reflectometry.types import Filename, ReflectivityOverQ, SampleRun - - -def curve(d, qmin, qmax): - return sc.DataArray(data=d, coords={'Q': sc.linspace('Q', qmin, qmax, len(d) + 1)}) - - -def test_reflectivity_curve_scaling(): - data = sc.concat( - ( - sc.ones(dims=['Q'], shape=[10], with_variances=True), - 0.5 * sc.ones(dims=['Q'], shape=[15], with_variances=True), - ), - dim='Q', - ) - data.variances[:] = 0.1 - - curves, factors = scale_reflectivity_curves_to_overlap( - (curve(data, 0, 0.3), curve(0.8 * data, 0.2, 0.7), curve(0.1 * data, 0.6, 1.0)), - ) - - assert_allclose(curves[0].data, data, rtol=sc.scalar(1e-5)) - assert_allclose(curves[1].data, 0.5 * data, rtol=sc.scalar(1e-5)) - assert_allclose(curves[2].data, 0.25 * data, rtol=sc.scalar(1e-5)) - np_assert_allclose((1, 0.5 / 0.8, 0.25 / 0.1), factors, 1e-4) - - -def test_reflectivity_curve_scaling_with_critical_edge(): - data = sc.concat( - ( - sc.ones(dims=['Q'], shape=[10], with_variances=True), - 0.5 * sc.ones(dims=['Q'], shape=[15], with_variances=True), - ), - dim='Q', - ) - data.variances[:] = 0.1 - - curves, factors = scale_reflectivity_curves_to_overlap( - ( - 2 * curve(data, 0, 0.3), - curve(0.8 * data, 0.2, 0.7), - curve(0.1 * data, 0.6, 1.0), - ), - critical_edge_interval=(sc.scalar(0.01), sc.scalar(0.05)), - ) - - assert_allclose(curves[0].data, data, rtol=sc.scalar(1e-5)) - assert_allclose(curves[1].data, 0.5 * data, rtol=sc.scalar(1e-5)) - assert_allclose(curves[2].data, 0.25 * data, rtol=sc.scalar(1e-5)) - np_assert_allclose((0.5, 0.5 / 0.8, 0.25 / 0.1), factors, 1e-4) - - -def test_combined_curves(): - qgrid = sc.linspace('Q', 0, 1, 26) - data = sc.concat( - ( - sc.ones(dims=['Q'], shape=[10], with_variances=True), - 0.5 * sc.ones(dims=['Q'], shape=[15], with_variances=True), - ), - dim='Q', - ) - data.variances[:] = 0.1 - curves = ( - curve(data, 0, 0.3), - curve(0.5 * data, 0.2, 0.7), - curve(0.25 * data, 0.6, 1.0), - ) - - combined = combine_curves(curves, qgrid) - assert_allclose( - combined.data, - sc.array( - dims='Q', - values=[ - 1.0, - 1, - 1, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.25, - 0.25, - 0.25, - 0.25, - 0.25, - 0.25, - 0.25, - 0.25, - 0.25, - 0.125, - 0.125, - 0.125, - 0.125, - 0.125, - 0.125, - ], - variances=[ - 0.1, - 0.1, - 0.1, - 0.1, - 0.1, - 0.02, - 0.02, - 0.025, - 0.025, - 0.025, - 0.025, - 0.025, - 0.025, - 0.025, - 0.025, - 0.005, - 0.005, - 0.00625, - 0.00625, - 0.00625, - 0.00625, - 0.00625, - 0.00625, - 0.00625, - 0.00625, - ], - ), - ) - - -@pytest.mark.filterwarnings("ignore:invalid value encountered in divide") -def test_combined_curves_resolution(): - qgrid = sc.linspace('Q', 0, 1, 26) - data = sc.concat( - ( - sc.ones(dims=['Q'], shape=[10], with_variances=True), - 0.5 * sc.ones(dims=['Q'], shape=[15], with_variances=True), - ), - dim='Q', - ) - data.variances[:] = 0.1 - curves = ( - curve(data, 0, 0.3), - curve(0.5 * data, 0.2, 0.7), - curve(0.25 * data, 0.6, 1.0), - ) - curves[0].coords['Q_resolution'] = sc.midpoints(curves[0].coords['Q']) / 5 - combined = combine_curves(curves, qgrid) - assert 'Q_resolution' in combined.coords - assert combined.coords['Q_resolution'][0] == curves[0].coords['Q_resolution'][1] - assert sc.isnan(combined.coords['Q_resolution'][-1]) - - -def test_linlogspace_linear(): - q_lin = linlogspace( - dim='qz', edges=[0.008, 0.08], scale='linear', num=50, unit='1/angstrom' - ) - expected = sc.linspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom') - assert sc.allclose(q_lin, expected) - - -def test_linlogspace_linear_list_input(): - q_lin = linlogspace( - dim='qz', edges=[0.008, 0.08], unit='1/angstrom', scale=['linear'], num=[50] - ) - expected = sc.linspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom') - assert sc.allclose(q_lin, expected) - - -def test_linlogspace_log(): - q_log = linlogspace( - dim='qz', edges=[0.008, 0.08], unit='1/angstrom', scale='log', num=50 - ) - expected = sc.geomspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom') - assert sc.allclose(q_log, expected) - - -def test_linlogspace_linear_log(): - q_linlog = linlogspace( - dim='qz', - edges=[0.008, 0.03, 0.08], - unit='1/angstrom', - scale=['linear', 'log'], - num=[16, 20], - ) - exp_lin = sc.linspace(dim='qz', start=0.008, stop=0.03, num=16, unit='1/angstrom') - exp_log = sc.geomspace(dim='qz', start=0.03, stop=0.08, num=21, unit='1/angstrom') - expected = sc.concat([exp_lin, exp_log['qz', 1:]], 'qz') - assert sc.allclose(q_linlog, expected) - - -def test_linlogspace_log_linear(): - q_loglin = linlogspace( - dim='qz', - edges=[0.008, 0.03, 0.08], - unit='1/angstrom', - scale=['log', 'linear'], - num=[16, 20], - ) - exp_log = sc.geomspace(dim='qz', start=0.008, stop=0.03, num=16, unit='1/angstrom') - exp_lin = sc.linspace(dim='qz', start=0.03, stop=0.08, num=21, unit='1/angstrom') - expected = sc.concat([exp_log, exp_lin['qz', 1:]], 'qz') - assert sc.allclose(q_loglin, expected) - - -def test_linlogspace_linear_log_linear(): - q_linloglin = linlogspace( - dim='qz', - edges=[0.008, 0.03, 0.08, 0.12], - unit='1/angstrom', - scale=['linear', 'log', 'linear'], - num=[16, 20, 10], - ) - exp_lin = sc.linspace(dim='qz', start=0.008, stop=0.03, num=16, unit='1/angstrom') - exp_log = sc.geomspace(dim='qz', start=0.03, stop=0.08, num=21, unit='1/angstrom') - exp_lin2 = sc.linspace(dim='qz', start=0.08, stop=0.12, num=11, unit='1/angstrom') - expected = sc.concat([exp_lin, exp_log['qz', 1:], exp_lin2['qz', 1:]], 'qz') - assert sc.allclose(q_linloglin, expected) - - -def test_linlogspace_bad_input(): - with pytest.raises(ValueError, match="Sizes do not match"): - _ = linlogspace( - dim='qz', - edges=[0.008, 0.03, 0.08, 0.12], - unit='1/angstrom', - scale=['linear', 'log'], - num=[16, 20], - ) - - -@pytest.mark.filterwarnings("ignore:No suitable") -def test_orso_datasets_tool(): - def normalized_ioq(filename: Filename[SampleRun]) -> ReflectivityOverQ: - return filename - - def orso_dataset(filename: Filename[SampleRun]) -> OrsoIofQDataset: - class Reduction: - corrections = [] # noqa: RUF012 - - return OrsoDataset( - Orso({}, Reduction, [], name=f'{filename}.orso'), np.ones((0, 0)) - ) - - workflow = sl.Pipeline( - [normalized_ioq, orso_dataset], params={Filename[SampleRun]: 'default'} - ) - datasets = orso_datasets_from_measurements( - workflow, - [{}, {Filename[SampleRun]: 'special'}], - scale_to_overlap=False, - ) - assert len(datasets) == 2 - assert tuple(d.info.name for d in datasets) == ('default.orso', 'special.orso') From 07eef7aa1334279db1b892aff13bf87481920245 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 18 Dec 2025 09:44:04 +0100 Subject: [PATCH 4/5] fix: use scalar --- src/ess/reflectometry/tools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ess/reflectometry/tools.py b/src/ess/reflectometry/tools.py index 355b43c6..b7077cae 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -443,7 +443,7 @@ def combine_curves( r = _interpolate_on_qgrid(map(sc.values, curves), q_bin_edges) v = _interpolate_on_qgrid(map(sc.variances, curves), q_bin_edges) - v = sc.where(v == 0, sc.scalar(np.nan, unit=v.unit), v) + v = sc.where(v == sc.scalar(0.0, unit=v.unit), sc.scalar(np.nan, unit=v.unit), v) inv_v = 1.0 / v r_avg = sc.nansum(r * inv_v, dim='curves') / sc.nansum(inv_v, dim='curves') v_avg = 1 / sc.nansum(inv_v, dim='curves') From 86225efc54b5d3c47a48fced31153ec03a10e7e4 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 18 Dec 2025 13:44:01 +0100 Subject: [PATCH 5/5] feat: add Estia resolution --- src/ess/estia/normalization.py | 6 +++--- src/ess/estia/resolution.py | 30 +++++++++++++++++++++--------- 2 files changed, 24 insertions(+), 12 deletions(-) diff --git a/src/ess/estia/normalization.py b/src/ess/estia/normalization.py index 82c99def..5140d44b 100644 --- a/src/ess/estia/normalization.py +++ b/src/ess/estia/normalization.py @@ -23,6 +23,7 @@ angular_resolution, q_resolution, sample_size_resolution, + wavelength_resolution, ) @@ -101,9 +102,8 @@ def evaluate_reference( ), { **graph, - # The wavelength resolution of Estia is not implemented yet - # when it is this will be replaced by the correct value. - "wavelength_resolution": lambda: sc.scalar(0.05), + "source_pulse_length": lambda: sc.scalar(3.0, unit='ms').to(unit='ns'), + "wavelength_resolution": wavelength_resolution, "sample_size_resolution": sample_size_resolution, "angular_resolution": angular_resolution, "Q_resolution": q_resolution, diff --git a/src/ess/estia/resolution.py b/src/ess/estia/resolution.py index 75205945..63f5f48b 100644 --- a/src/ess/estia/resolution.py +++ b/src/ess/estia/resolution.py @@ -1,12 +1,15 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) import scipp as sc +from scipp import constants from ..reflectometry.tools import fwhm_to_std def wavelength_resolution( - # What parameters are needed? + source_pulse_length: sc.Variable, + wavelength: sc.Variable, + Ltotal: sc.Variable, ): """ Find the wavelength resolution contribution of the ESTIA instrument. @@ -14,18 +17,27 @@ def wavelength_resolution( Parameters ---------- - L1: - Distance from midpoint between choppers to sample. - L2: - Distance from sample to detector. + source_pulse_length: + The length of the ESS source pulse. + wavelength: + The wavelength of the neutron. + Ltotal: + The distance from source to detector. Returns ------- : - The wavelength resolution variable, as standard deviation. + The wavelength resolution, + ratio of the standard deviation of wavelength and the wavelength. """ - # Don't yet know how to do this - raise NotImplementedError() + # The exact factor depends on the shape of the source pulse. + # The factor here assumes a rectangular source pulse: + # import sympy as sp + # x, D = sp.symbols('x, D', positive=True) + # sp.sqrt(sp.integrate(1/D * (x-D/2)**2, (x, 0, D))) + standard_deviation_of_time_at_source = (1 / 12**0.5) * source_pulse_length + tof = wavelength * Ltotal * constants.m_n / constants.h + return (standard_deviation_of_time_at_source / tof).to(unit='dimensionless') def sample_size_resolution( @@ -72,7 +84,7 @@ def angular_resolution( Returns ------- : - Angular resolution standard deviation + Angular resolution standard deviation over the reflection angle. """ return ( fwhm_to_std(