Add function which calculates Polarization & propagating uncertainties #97

esmith1729 · 2025-01-30T16:10:31Z

Polarization/asymmetry #58

…g, refactoring of function.

Tom-Willemsen

Most of these things are relatively minor - core implementation looks good.

src/ibex_bluesky_core/devices/simpledae/reducers.py

Tom-Willemsen · 2025-02-07T14:12:29Z

tests/devices/simpledae/test_reducers.py

+    ("a", "b", "variance_a", "variance_b", "expected_value", "expected_uncertainty"),
+    [
+        # Case 1: Symmetric case with equal uncertainties
+        (5.0, 3.0, 0.1, 0.1, 0.25, 0.05762215285808055),  # comment where tf this number came from


comment where tf this number came from

good idea

tests/devices/simpledae/test_reducers.py

Tom-Willemsen · 2025-02-07T14:13:52Z

tests/devices/simpledae/test_reducers.py

+        # Case 2: Asymmetric case with different uncertainties
+        (10.0, 6.0, 0.2, 0.3, 0.25, 0.04764984588117782),
+        # Case 3: Case with larger values and different uncertainty magnitudes
+        (100.0, 60.0, 1.0, 2.0, 0.25, 0.0120017902310447),


In our chat we also had some "simple" example cases - the sort which could be worked out easily "in your head". Those would make good test cases to add here.

Tom-Willemsen · 2025-02-07T14:15:30Z

tests/devices/simpledae/test_reducers.py

+    var_a = sc.Variable(dims=[], values=a, variances=variance_a, unit="", dtype="float64")
+    var_b = sc.Variable(dims=[], values=b, variances=variance_b, unit="", dtype="float64")


From https://scipp.github.io/generated/classes/scipp.Variable.html#scipp.Variable:

This constructor is meant primarily for internal use. Use one of the Specialized creation functions instead. See in particular scipp.array() and scipp.scalar().

(scalar is what you want here) - it still gives you a sc.Variable object in the end though.

I added sc.scalar in the tests where it had Variable... there was an error arising in the assert result_value == pytest.approx(expected_value) lines, but I changed it to assert result_value.value == pytest.approx(expected_value) instead, and that seems to have silenced Pylance at least

Tom-Willemsen · 2025-02-07T14:17:50Z

tests/devices/simpledae/test_reducers.py

+        polarization(var_a, var_b)
+
+
+# test cases uncert = 0


Tom-Willemsen · 2025-02-07T14:18:29Z

tests/devices/simpledae/test_reducers.py

+@pytest.mark.parametrize(
+    ("a", "b", "variances_a", "variances_b", "expected_values", "expected_uncertainties"),
+    [
+        # Case 1: Symmetric case with equal uncertainties


Not sure what "symmetric" means here - could you clarify?

Also the equal uncertainties seems not true.

Copy/paste?

symmetric means that the uncertainty values for both a and b are the same. It seems not true in that the value is wrong? I need to check, yes this was copy/paste but I think it might have been correlated uncertainty.

Tom-Willemsen · 2025-02-07T14:30:24Z

doc/devices/dae.md

+Polarization refers to the property of transverse waves which specifies the geometrical orientation of the 
+oscillations. 
+
+The polarization funtion provided will calculate the polarization between two values, A and B, which 
+have different definitions based on the instrument context.
+
+Instrument-Specific Interpretations
+SANS Instruments (e.g., LARMOR)
+A: Intensity in DAE period before switching a flipper.
+B: Intensity in DAE period after switching a flipper.
+
+Reflectometry Instruments (e.g., POLREF)
+Similar to LARMOR, A and B represent intensities before and after flipper switching.
+
+Muon Instruments
+A and B refer to Measurements from different detector banks.


I think we need to re-word this section a bit. Not your fault, this is kind of specific science knowledge here - which I'm no expert on myself but have picked up bits & pieces from talking to scientists.

How about something like:

Polarization/asymmetry

ibex_bluesky_core provides a helper method,
{py:obj}ibex_bluesky_core.devices.simpledae.reducers.polarization, for calculating the quantity (a-b)/(a+b). This quantity is used, for example, in neutron polarization measurements, and in calculating asymmetry for muon measurements.

For this expression, scipp's default uncertainty propagation rules cannot be used as the uncertainties on (a-b) are correlated with those of (a+b) in the division step - but scipp assumes uncorrelated data. This helper method calculates the uncertainties following linear error propagation theory, using the partial derivatives of the above expression.

The partial derivatives are:
[INSERT MATHS HERE - and see the fitting models docs pages for how to render nice maths equations in docs]
Which then means the variances computed by this helper function are:
[INSERT MATHS HERE]

Co-authored-by: Tom Willemsen <tom.willemsen@stfc.ac.uk>

…etry

Tom-Willemsen

Couple of tiny comments and then we're there I think...

Tom-Willemsen · 2025-06-17T15:13:14Z

doc/devices/dae.md

+Muon Instruments
+A and B refer to Measurements from different detector banks.
+
+[`polarization`](ibex_bluesky_core.devices.simpledae.polarization)


Not sure this is needed - there's a link to the function in the first sentence?

Tom-Willemsen · 2025-06-17T15:25:53Z

tests/devices/simpledae/test_reducers.py

+    a_ufloat = ufloat(a, variance_a)
+    b_ufloat = ufloat(b, variance_b)
+
+    # polarization value, i.e. (a - b) / (a + b)
+    polarization_ufloat = (a_ufloat.n - b_ufloat.n) / (a_ufloat.n + b_ufloat.n)
+
+    # the partial derivatives of a and b, calculated with 'uncertainties' library's ufloat type
+    partial_a = (2 * b_ufloat.n) / ((a_ufloat.n + b_ufloat.n) ** 2)
+    partial_b = (-2 * a_ufloat.n) / ((a_ufloat.n + b_ufloat.n) ** 2)
+
+    # variance calculated with 'uncertainties' library
+    variance = (partial_a**2 * a_ufloat.s) + (partial_b**2 * b_ufloat.s)
+    uncertainty = variance**0.5  # uncertainty is sqrt of variance
+
+    # Two scipp scalars, to test our polarization function
+    var_a = sc.scalar(value=a, variance=variance_a, unit="", dtype="float64")
+    var_b = sc.scalar(value=b, variance=variance_b, unit="", dtype="float64")
+    result_value = polarization(var_a, var_b)
+    result_uncertainy = (result_value.variance) ** 0.5  # uncertainty is sqrt of variance


I don't quite understand what's going on in this test; we shouldn't need partial derivatives with the uncertainties library as that's what the library is doing for us. I think conceptually it should be:

a_ufloat = ... b_ufloat = ... a_scipp = ... b_scipp = ... pol_ufloat = (a-b)/(a+b) pol_scipp = polarization(a, b) assert pol_ufloat.nominal == pol_scipp.value assert pol_ufloat.stddev == sqrt(pol_scipp.variance)

i.e. similar to the test below

whoops yes sorry I meant to remove the unnecessary bits and forgot, will do it now.

Actually sorry, pretty sure we need those; when calculating the polarization from the nominal values, the variances do not get carried into that value, so we need to calculate it in long form. It was easier with the arrays, because it was just adding an array of ufloats. But in this test we are calculating with the ufloat's nominal values themselves.

for whatever reason, ufloat library really didn't like doing this

polarization_ufloat = (a_ufloat - b_ufloat) / (a_ufloat + b_ufloat)

And needs it to be nominal values.... which then creates polarization_ufloat as a "float" and not a ufloat. We could possibly get around this by instead making two arrays of one element, sort of copying the array test, but it might be more confusing in future for the next person looking at this code to do it that way.

Initial polarization func and test

88202fa

pull-request-size bot added the size/M label Jan 30, 2025

esmith1729 added 5 commits January 30, 2025 16:13

Merge branch 'main' into 58_polarization_asymmetry

46cee18

Clean up polarization function

bce1bbe

parametrize and write test class for polarization

60904c0

Add polarization doc comment, some changes to edge case error checkin…

a3c5fbf

…g, refactoring of function.

Add more unit tests

b1b2e7e

pull-request-size bot added size/L and removed size/M labels Feb 5, 2025

esmith1729 added 2 commits February 6, 2025 11:54

Add some more explanation to code doc

c3e53d8

Add polarization to API reference

6bcfac1

esmith1729 linked an issue Feb 6, 2025 that may be closed by this pull request

Polarization/asymmetry #58

Closed

5 tasks

esmith1729 changed the title ~~Initial polarization func and test~~ Polarization/asymmetry function Feb 6, 2025

esmith1729 added 2 commits February 6, 2025 17:26

make strings raw to satisfy ruff complaints

1c2e292

Merge branch 'main' into 58_polarization_asymmetry

bbfd749

Tom-Willemsen reviewed Feb 7, 2025

View reviewed changes

esmith1729 and others added 14 commits March 17, 2025 10:02

Update src/ibex_bluesky_core/devices/simpledae/reducers.py

b4ff8ed

Co-authored-by: Tom Willemsen <tom.willemsen@stfc.ac.uk>

Update src/ibex_bluesky_core/devices/simpledae/reducers.py

68be42e

Co-authored-by: Tom Willemsen <tom.willemsen@stfc.ac.uk>

Merge branch 'main' into 58_polarization_asymmetry

8cd8a33

Change type hint to include DataArray for polarization

a9ee9a2

Fix polarization info in docs

bc707fb

Update polarization info to include DataArray

897d12f

Add mathematical formula for variance calculation

e2a42aa

Add uncertainties as dev dependency

968f77d

Change polarization signature type hint and description

f363136

refactor polarization tests to compare against uncertainties library

2fc7803

Fix imports

68d1744

Merge remote-tracking branch 'origin/main' into 58_polarization_asymm…

7ba2d5a

…etry

add polarization to import

7bd51d0

reformat API docs to link to scipp and scipp neutron docs

9e2b122

esmith1729 added 5 commits June 13, 2025 15:13

fix merge conflict with imports

e7d10f4

add sum_spectra to API

8632f3a

ruff format

08855ef

reformat function doc

93acddd

remove pyobject DaeSpectra

97628ec

Tom-Willemsen reviewed Jun 17, 2025

View reviewed changes

esmith1729 and others added 4 commits June 19, 2025 09:53

Remove unnecessary link to polarization

2f5db04

remove unnecessary calculations in test

67e90cb

revert changes to polarization calculation function

cbb43ae

Merge branch 'main' into 58_polarization_asymmetry

24fb6d0

jackbdoughty approved these changes Jun 19, 2025

View reviewed changes

jackbdoughty added the bluesky-Semver-Minor label Jun 19, 2025

jackbdoughty changed the title ~~Polarization/asymmetry function~~ Add function which calculates Polarization & propagating uncertainties Jun 19, 2025

jackbdoughty merged commit 240ba51 into main Jun 19, 2025
13 checks passed

jackbdoughty deleted the 58_polarization_asymmetry branch June 19, 2025 14:26

		var_a = sc.Variable(dims=[], values=a, variances=variance_a, unit="", dtype="float64")
		var_b = sc.Variable(dims=[], values=b, variances=variance_b, unit="", dtype="float64")

Add function which calculates Polarization & propagating uncertainties #97

Add function which calculates Polarization & propagating uncertainties #97

Uh oh!

Conversation

esmith1729 commented Jan 30, 2025

Uh oh!

Tom-Willemsen left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

esmith1729 Jun 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Polarization/asymmetry

Uh oh!

Tom-Willemsen left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

esmith1729 Jun 5, 2025 •

edited

Loading