Polarization correction #122

jokasimr · 2025-02-27T14:04:30Z

Polarization correction for ESTIA

Only important differences from #108 are in src/ess/estia/calibration.py: https://github.com/scipp/essreflectometry/pull/122/files#diff-df676f2a67f8251976ea2fcc161ea759d49deef6ef66c5adbf01d10665230077

SimonHeybrock · 2025-02-28T04:59:31Z

src/ess/estia/calibration.py

Can we use esspolarization to perform the correction?

Right now we can't because of scipp/esspolarization#84

When that is fixed we should be able to use esspolarization to do the correction.

…sking

src/ess/estia/calibration.py

SimonHeybrock · 2025-04-25T08:55:12Z

src/ess/estia/calibration.py

+)
+
+
+def solve_for_calibration_parameters(Io, Is):


Missing type hints everywhere.

SimonHeybrock · 2025-04-25T08:57:01Z

src/ess/estia/calibration.py

+    Pp = sc.sqrt(
+        (Ipp + Iaa - Ipa - Iap)
+        * (alp * (Ipp - Iap) - Iaa + Ipa)
+        / (
+            (rho * alp * Ipp + Iaa + rho * Ipa + alp * Iap)
+            * (rho * (Ipp - Ipa) - Iaa + Iap)
+        )
+    )
+    Ap = sc.sqrt(
+        (Ipp + Iaa - Ipa - Iap)
+        * (rho * (Ipp - Ipa) - Iaa + Iap)
+        / (
+            (rho * alp * Ipp + Iaa + rho * Ipa + alp * Iap)
+            * (alp * (Ipp - Iap) - Iaa + Ipa)
+        )
+    )
+    Rs = sc.sqrt(


A lot of duplication between this equations. Can we share some helper vars?

I think it's easier to read this way, the same way it's written as in the article.
Some shared vars would improve performance, but that would be a small improvement relative to the runtime of the workflow.

Well, I find it harder to read.

There's a practical reason to keep it similar to the description in the article. The code is just a solution to an equation system, it's pretty meaningless to read it without comparing to what is in the article, it is verified by comparison to the article, and keeping it similar to the description there makes it easier to verify.

If you still think it should be rewritten for readability I'm of course open to review any suggestions.

As suggested, you can define a variable, e.g., for long repeated bits such as (rho * alp * Ipp + Iaa + rho * Ipa + alp * Iap). While you correctly argue that the math equations do repeat it, those are typically way more compact than code and easier to "parse".

Please rewrite it in a way that you find easier to read, make a commit suggestion, and I'll approve that. It's a waste of time for both of us that I try to rewrite it the way you prefer it.

SimonHeybrock · 2025-04-25T08:58:34Z

src/ess/estia/calibration.py

+    return I0, Pp, Pa, Ap, Aa, Rspp, Rsaa
+
+
+def correction_matrix(Pp, Pa, Ap, Aa):


Many params in return and signature. Looks error prone. Maybe we can use keyword args, or a simple dataclass?

SimonHeybrock · 2025-04-25T09:02:34Z

src/ess/estia/calibration.py

+    return [
+        [
+            (1 + Pp) * (1 + Ap) / 4,
+            (1 + Pp) * (1 - Ap) / 4,
+            (1 - Pp) * (1 + Ap) / 4,
+            (1 - Pp) * (1 - Ap) / 4,
+        ],
+        [
+            (1 + Pp) * (1 + Aa) / 4,
+            (1 + Pp) * (1 - Aa) / 4,
+            (1 - Pp) * (1 + Aa) / 4,
+            (1 - Pp) * (1 - Aa) / 4,
+        ],
+        [
+            (1 + Pa) * (1 + Ap) / 4,
+            (1 + Pa) * (1 - Ap) / 4,
+            (1 - Pa) * (1 + Ap) / 4,
+            (1 - Pa) * (1 - Ap) / 4,
+        ],
+        [
+            (1 + Pa) * (1 + Aa) / 4,
+            (1 + Pa) * (1 - Aa) / 4,
+            (1 - Pa) * (1 + Aa) / 4,
+            (1 - Pa) * (1 - Aa) / 4,
+        ],
+    ]


Would this be more readable if we define two 2x2 matrices and use something like https://numpy.org/doc/stable/reference/generated/numpy.kron.html?

SimonHeybrock · 2025-04-25T09:05:34Z

src/ess/estia/calibration.py

+def compute_reflectivity_calibrate_on_q(
+    reference_supermirror,
+    reference_polarized_supermirror,
+    sample,
+    qbins,
+):
+    """Reduces reference and sample to Q before applying
+    the polarization correction and normalization."""
+    reference_supermirror = [reduce_to_q(i, qbins) for i in reference_supermirror]
+    reference_polarized_supermirror = [
+        reduce_to_q(i, qbins) for i in reference_polarized_supermirror
+    ]
+    I0, C = calibration_factors_from_reference_measurements(
+        reference_supermirror, reference_polarized_supermirror
+    )


Here and in quite a few other functions: Please use (and enforce, if possible) function args to be keyword args. There are many simple mess-ups than can be avoided that way.

SimonHeybrock · 2025-04-25T09:08:05Z

src/ess/estia/types.py

+
+
+class Intensity(
+    sciline.ScopeTwoParams[PolarizedRunType, FlipperSetting, sc.DataArray], sc.DataArray


Can now use Scope if you update sciline.

SimonHeybrock · 2025-04-25T09:14:06Z

src/ess/estia/types.py

+OffOff = NewType("OffOff", str)
+OffOn = NewType("OffOn", str)
+OnOff = NewType("OnOff", str)
+OnOn = NewType("OnOn", str)
+FlipperSetting = TypeVar("FlipperSetting", OffOff, OffOn, OnOff, OnOn)


Is it the flipper setting that matters to the workflow, or the polarization (which, as far as I understand, by be the reverse, depending on how the supermirror is setup)?

I think either works, it's just labels, but in my opinion it's more natural to talk about the flipper setting to distinguish those labels from the labels on the resulting reflectivity curves.

But overall I'm not sure this is the nicest way to define the domain types in the workflow.
We probably need to get some experience using it to see how convenient/cumbersome it is to work with in practice.

Can you use the types defined in esspolarization?

SimonHeybrock · 2025-04-25T09:16:20Z

src/ess/reflectometry/normalization.py

@@ -20,6 +20,19 @@
 )


+def reduce_to_q(da, qbins):
+    if da.bins:
+        return da.bins.concat().bin(Q=qbins)


Is the bins.concat() necessary? Might be better to use the newer dim arg to bin instead.

…mple

SimonHeybrock · 2025-04-30T02:33:46Z

src/ess/estia/calibration.py

+def _kronecker_product(A, B):
+    return [
+        [A[ia][ja] * B[ib][jb] for ja in range(2) for jb in range(2)]
+        for ia in range(2)
+        for ib in range(2)
+    ]


Can you clarify why we cannot use np.kron?

The 2D arrays containing DataArray are treated by np.kron as 3D arrays, and that makes the result different from what we want.

SimonHeybrock · 2025-04-30T02:35:56Z

src/ess/estia/calibration.py

+        '''Corrects the intensities for imperfections of polarizing components'''
+        return _linsolve(self.polarization_matrix, Is)
+
+    def correct_intensities_and_normalize_by_reference(self, Is):


Still missing type hints in many places.

Co-authored-by: Simon Heybrock <[email protected]>

SimonHeybrock reviewed Feb 28, 2025

View reviewed changes

jokasimr force-pushed the polcal branch from 5310d7e to 1e2baa9 Compare March 18, 2025 08:13

jokasimr added 9 commits March 18, 2025 15:04

docs: polarization calibration prototype

8768696

docs: include figures

f78be6c

docs: add header

0c686af

changes

a25ef55

feat: implement calibration and reduction

d22df0f

refactor: move test code to tests, add docstrings

e3bb2c1

feat: draft sciline domain types

2e50d3b

test: solve recovers expected parameters

7b70b11

test: stacking logic in linsolve

ab78efb

jokasimr force-pushed the polcal branch from 9755352 to ab78efb Compare March 18, 2025 14:04

jokasimr added 18 commits March 19, 2025 14:11

fix: remove unused

fc20fd0

fix: add position already in reduction to lz

893abdb

fix: calibration on lz grid, compute Q on lz grid for sample + add ma…

0fb609a

…sking

docs: add esspolarization + tox -e deps

727bf44

docs: add example reducing polarized mcstas data

ec48ec6

fix: move factor from normalization to polarization matrix

e9acb7c

fix: use pandas to load csv

6e11c7d

fix: tweak center of incident angle for mcstas data

381227e

fix: old value was better

0329a22

fix: tests

179aaef

feat: make normalization able to handle different kinds of references

261ac50

fix: no need to have two reduce to Q functions

cca0ce5

docs: stap template

4930b00

fix: remove unnecessary

725ee1d

deps: add ipympl to docs

91fb8cd

Merge branch 'polcal' of github.com:scipp/essreflectometry into polcal

8d9ba1d

docs: remove polarization notebook

6089d41

docs: remove polarization notebook from docs

f8083b2

jokasimr added 2 commits April 9, 2025 09:50

docs: remove polarization correction experiment notebook

d7d31c0

Merge branch 'main' into polcal

b187fbd

jokasimr marked this pull request as ready for review April 9, 2025 11:24

jokasimr added 2 commits April 9, 2025 13:24

docs: remove figures

6e1a83f

Merge branch 'polcal' of github.com:scipp/essreflectometry into polcal

1c92913

SimonHeybrock reviewed Apr 25, 2025

View reviewed changes

jokasimr added 10 commits April 25, 2025 15:04

fix: use np.kron to construct calibration matrix

c5610b7

deps: use latest sciline for new Scope

0913455

fix: use dims arg to do concatenation

4099486

refactor: implement container for polarization calibration parameters

f6764bb

fix: better doc comment + fix typo

be37702

Merge remote-tracking branch 'origin/main' into polcal

76bcf3c

Merge branch 'main' into polcal

7d1753c

fix: reimplement kron, numpy and scipp doesn't play nicely in this case

ca28f99

fix: don't assume only dimension of 1d reference is Q

e0db8f0

docs: add notebook and data for test case with simulated spin flip sa…

e72eb85

…mple

SimonHeybrock reviewed Apr 30, 2025

View reviewed changes

jokasimr and others added 6 commits April 30, 2025 08:11

Update src/ess/estia/calibration.py

dc36161

Co-authored-by: Simon Heybrock <[email protected]>

Apply automatic formatting

3f7182b

docs: add new example notebook to index

9aaf79d

Merge branch 'polcal' of github.com:scipp/essreflectometry into polcal

e801737

docs: fix hader

889ecc5

docs: fix sample names

8d071c7

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Polarization correction #122

Polarization correction #122

jokasimr commented Feb 27, 2025 •

edited

Loading

SimonHeybrock Feb 28, 2025

jokasimr Apr 9, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 25, 2025

jokasimr Apr 25, 2025

SimonHeybrock Apr 25, 2025

jokasimr Apr 25, 2025 •

edited

Loading

SimonHeybrock Apr 30, 2025

jokasimr Apr 30, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 25, 2025

jokasimr Apr 25, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 25, 2025

SimonHeybrock Apr 30, 2025

jokasimr Apr 30, 2025

SimonHeybrock Apr 30, 2025

		return I0, Pp, Pa, Ap, Aa, Rspp, Rsaa


		def correction_matrix(Pp, Pa, Ap, Aa):



		class Intensity(
		sciline.ScopeTwoParams[PolarizedRunType, FlipperSetting, sc.DataArray], sc.DataArray

Polarization correction #122

Are you sure you want to change the base?

Polarization correction #122

Conversation

jokasimr commented Feb 27, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jokasimr Apr 25, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jokasimr commented Feb 27, 2025 •

edited

Loading

jokasimr Apr 25, 2025 •

edited

Loading