From 92050871e2b06b9b8612465b5af1e63a8508a08a Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 14 Oct 2024 15:18:00 +0200 Subject: [PATCH 1/6] feat: scale critical edge to 1 --- src/ess/reflectometry/tools.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/ess/reflectometry/tools.py b/src/ess/reflectometry/tools.py index 341fa4f0..a8f57105 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -171,6 +171,7 @@ def _interpolate_on_qgrid(curves, grid): def scale_reflectivity_curves_to_overlap( curves: Sequence[sc.DataArray], return_scaling_factors=False, + critical_edge_interval=None, ) -> list[sc.DataArray] | list[sc.scalar]: '''Make the curves overlap by scaling all except the first by a factor. The scaling factors are determined by a maximum likelihood estimate @@ -192,6 +193,20 @@ def scale_reflectivity_curves_to_overlap( : A list of scaled reflectivity curves or a list of scaling factors. ''' + if critical_edge_interval is not None: + q = next(iter(curves)).coords['Q'] + N = ( + ((q > critical_edge_interval[0]) & (q < critical_edge_interval[1])) + .sum() + .value + ) + edge = sc.DataArray( + data=sc.ones(dims=('Q',), shape=(N,), with_variances=True), + coords={'Q': sc.linspace('Q', *critical_edge_interval, N + 1)}, + ) + return scale_reflectivity_curves_to_overlap( + [edge, *curves], return_scaling_factors=return_scaling_factors + )[1:] if len({c.data.unit for c in curves}) != 1: raise ValueError('The reflectivity curves must have the same unit') if len({c.coords['Q'].unit for c in curves}) != 1: From 6773d82b87eb5f0f07a1fd82b9a1e2ef93b1379d Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 18 Oct 2024 13:44:04 +0200 Subject: [PATCH 2/6] test --- tests/tools_test.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/tools_test.py b/tests/tools_test.py index b447fcf4..1c9745d4 100644 --- a/tests/tools_test.py +++ b/tests/tools_test.py @@ -29,6 +29,30 @@ def test_reflectivity_curve_scaling(): assert_allclose(curves[2].data, 0.25 * data, rtol=sc.scalar(1e-5)) +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 = 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)) + + def test_reflectivity_curve_scaling_return_factors(): data = sc.concat( ( From 6b85d5c0d894d077daa1ba9eed69d97d20f4cfeb Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 18 Oct 2024 14:01:45 +0200 Subject: [PATCH 3/6] fix: return curves and edges --- docs/user-guide/amor/amor-reduction.ipynb | 2 +- src/ess/reflectometry/tools.py | 28 +++++++++++------------ tests/tools_test.py | 27 ++++------------------ 3 files changed, 19 insertions(+), 38 deletions(-) diff --git a/docs/user-guide/amor/amor-reduction.ipynb b/docs/user-guide/amor/amor-reduction.ipynb index 01ffaca2..04115d19 100644 --- a/docs/user-guide/amor/amor-reduction.ipynb +++ b/docs/user-guide/amor/amor-reduction.ipynb @@ -172,7 +172,7 @@ "from ess.reflectometry.tools import scale_reflectivity_curves_to_overlap\n", "results_scaled = dict(zip(\n", " results.keys(),\n", - " scale_reflectivity_curves_to_overlap(results.values()),\n", + " scale_reflectivity_curves_to_overlap(results.values())[0],\n", " strict=True\n", "))\n", "sc.plot(results_scaled, norm='log', vmin=1e-5)" diff --git a/src/ess/reflectometry/tools.py b/src/ess/reflectometry/tools.py index a8f57105..2b4893c8 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -170,28 +170,29 @@ def _interpolate_on_qgrid(curves, grid): def scale_reflectivity_curves_to_overlap( curves: Sequence[sc.DataArray], - return_scaling_factors=False, - critical_edge_interval=None, -) -> list[sc.DataArray] | list[sc.scalar]: + critical_edge_interval: tuple[sc.Variable, sc.Variable] | None = None, +) -> tuple[list[sc.DataArray], list[sc.scalar]]: '''Make the curves overlap by scaling all except the first by a factor. The scaling factors are determined by a maximum likelihood estimate (assuming the errors are normal distributed). + If :code:`critical_edge_interval` is provided then all curves are scaled. + All curves must be have the same unit for data and the Q-coordinate. Parameters --------- curves: the reflectivity curves that should be scaled together - return_scaling_factor: - If True the return value of the function - is a list of the scaling factors that should be applied. - If False (default) the function returns the scaled curves. + critical_edge_interval: optional, + a tuple denoting an interval that is known to belong + to the critical edge, i.e. where the reflectivity is + known to be 1. Returns --------- : - A list of scaled reflectivity curves or a list of scaling factors. + A list of scaled reflectivity curves and a list of the scaling factors. ''' if critical_edge_interval is not None: q = next(iter(curves)).coords['Q'] @@ -204,9 +205,8 @@ def scale_reflectivity_curves_to_overlap( data=sc.ones(dims=('Q',), shape=(N,), with_variances=True), coords={'Q': sc.linspace('Q', *critical_edge_interval, N + 1)}, ) - return scale_reflectivity_curves_to_overlap( - [edge, *curves], return_scaling_factors=return_scaling_factors - )[1:] + curves, factors = scale_reflectivity_curves_to_overlap([edge, *curves]) + return curves[1:], factors[1:] if len({c.data.unit for c in curves}) != 1: raise ValueError('The reflectivity curves must have the same unit') if len({c.coords['Q'].unit for c in curves}) != 1: @@ -229,13 +229,11 @@ def cost(scaling_factors): return np.nansum((r_scaled - r_avg) ** 2 * inv_v_scaled) sol = opt.minimize(cost, [1.0] * (len(curves) - 1)) - scaling_factors = (1.0, *sol.x) - if return_scaling_factors: - return [sc.scalar(x) for x in scaling_factors] + scaling_factors = (1.0, *map(float, sol.x)) return [ scaling_factor * curve for scaling_factor, curve in zip(scaling_factors, curves, strict=True) - ] + ], scaling_factors def combine_curves( diff --git a/tests/tools_test.py b/tests/tools_test.py index 1c9745d4..193e9245 100644 --- a/tests/tools_test.py +++ b/tests/tools_test.py @@ -1,6 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) import scipp as sc +from numpy.testing import assert_allclose as np_assert_allclose from scipp.testing import assert_allclose from ess.reflectometry.tools import combine_curves, scale_reflectivity_curves_to_overlap @@ -20,13 +21,14 @@ def test_reflectivity_curve_scaling(): ) data.variances[:] = 0.1 - curves = scale_reflectivity_curves_to_overlap( + 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(): @@ -39,7 +41,7 @@ def test_reflectivity_curve_scaling_with_critical_edge(): ) data.variances[:] = 0.1 - curves = scale_reflectivity_curves_to_overlap( + curves, factors = scale_reflectivity_curves_to_overlap( ( 2 * curve(data, 0, 0.3), curve(0.8 * data, 0.2, 0.7), @@ -51,26 +53,7 @@ def test_reflectivity_curve_scaling_with_critical_edge(): 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)) - - -def test_reflectivity_curve_scaling_return_factors(): - 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 - - 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)), - return_scaling_factors=True, - ) - - assert_allclose(factors[0], sc.scalar(1.0), rtol=sc.scalar(1e-5)) - assert_allclose(factors[1], sc.scalar(0.5 / 0.8), rtol=sc.scalar(1e-5)) - assert_allclose(factors[2], sc.scalar(0.25 / 0.1), 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(): From af1eb590dfa995d6850849eed67678343ba78234 Mon Sep 17 00:00:00 2001 From: jokasimr Date: Mon, 21 Oct 2024 09:17:10 +0200 Subject: [PATCH 4/6] Update src/ess/reflectometry/tools.py Co-authored-by: Simon Heybrock <12912489+SimonHeybrock@users.noreply.github.com> --- 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 2b4893c8..564309db 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -197,7 +197,7 @@ def scale_reflectivity_curves_to_overlap( if critical_edge_interval is not None: q = next(iter(curves)).coords['Q'] N = ( - ((q > critical_edge_interval[0]) & (q < critical_edge_interval[1])) + ((q >= critical_edge_interval[0]) & (q < critical_edge_interval[1])) .sum() .value ) From e2743f7af251f52b16cc84846070a158acf0900b Mon Sep 17 00:00:00 2001 From: jokasimr Date: Mon, 21 Oct 2024 09:17:21 +0200 Subject: [PATCH 5/6] Update src/ess/reflectometry/tools.py Co-authored-by: Simon Heybrock <12912489+SimonHeybrock@users.noreply.github.com> --- 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 564309db..8843f89f 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -171,7 +171,7 @@ def _interpolate_on_qgrid(curves, grid): def scale_reflectivity_curves_to_overlap( curves: Sequence[sc.DataArray], critical_edge_interval: tuple[sc.Variable, sc.Variable] | None = None, -) -> tuple[list[sc.DataArray], list[sc.scalar]]: +) -> tuple[list[sc.DataArray], list[sc.Variable]]: '''Make the curves overlap by scaling all except the first by a factor. The scaling factors are determined by a maximum likelihood estimate (assuming the errors are normal distributed). From 8c05477aec7b4ac1866adf4a1a002c3db2ed798b Mon Sep 17 00:00:00 2001 From: jokasimr Date: Mon, 21 Oct 2024 09:17:30 +0200 Subject: [PATCH 6/6] Update src/ess/reflectometry/tools.py Co-authored-by: Simon Heybrock <12912489+SimonHeybrock@users.noreply.github.com> --- 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 8843f89f..86823efc 100644 --- a/src/ess/reflectometry/tools.py +++ b/src/ess/reflectometry/tools.py @@ -184,7 +184,7 @@ def scale_reflectivity_curves_to_overlap( --------- curves: the reflectivity curves that should be scaled together - critical_edge_interval: optional, + critical_edge_interval: a tuple denoting an interval that is known to belong to the critical edge, i.e. where the reflectivity is known to be 1.