Simplify AMOR/reflectometry workflow #107
Conversation
There was a problem hiding this comment.
- There's one or two hardcoded constants that I haven't gotten around to move and name properly. Those should of course get proper names and be documented.
Please do this. But not here, instead:
This PR is very large and after more than 1h and 13 comments, I am still no where near reviewing it all. Can you please split it into smaller PRs.
It is particularly difficult to tell for a change whether something was simply renamed or moved or modified in a meaningful way. Also, this PR seems to contain changes related to several aspects of the workflow such as loading, saving, coord transformations, normalisations, etc.
| @@ -111,11 +93,6 @@ class SampleSize(sciline.Scope[RunType, sc.Variable], sc.Variable): | |||
| """Size of the sample. If None it is assumed to be the same as the reference.""" | |||
|
|
|||
|
|
|||
| Gravity = NewType("Gravity", sc.Variable) | |||
| """This parameter determines if gravity is taken into account | |||
| when computing the scattering angle and momentum transfer.""" | |||
There was a problem hiding this comment.
Why did you remove this? Now the user can no longer control the gravity correction.
| BeamDivergenceLimits: ( | ||
| sc.scalar(-0.75, unit='deg').to(unit='rad'), | ||
| sc.scalar(0.75, unit='deg').to(unit='rad'), | ||
| ), |
There was a problem hiding this comment.
This is not very user friendly. Ideally, the code that uses these should convert to the unit it needs.
There was a problem hiding this comment.
On the other hand it's messy and somewhat error prone to have unit conversion code all over the package (forgetting a unit conversion somewhere might lead to UnitErrors and having extraneous unit conversions in every provider obscures the intention of the code, which is also bad user experience) 🤔
It's easiest to have unit conversion code on the boundaries of the application and only work with a set of predetermined units in the "internal" code.
One way to achieve that would be to create a separate provider just for converting the unit of the input to the common unit. But doing that for every input would clutter the workflow graph and make it considerably less readable.
There was a problem hiding this comment.
It's easiest to have unit conversion code on the boundaries of the application and only work with a set of predetermined units in the "internal" code.
I disagree. It is often very difficult to understand which units are required by any given piece of code. And we don't encode them in our interfaces.
forgetting a unit conversion somewhere might lead to UnitErrors
Yes, the same as using an incorrect unit in a parameter. Except that the latter UnitErrors are user errors and it is typically very difficulty to figure out which unit was wrong because those errors can happen after many processing steps.
having extraneous unit conversions in every provider obscures the intention of the code, which is also bad user experience
True. That is a tradeoff. But we accepted that tradeoff in ScippNeutron in favour of reducing the number of user-facing UnitErrors. In fact, we took great care to make coord transform kernels work with any units and do conversions in an optimal way.
I don't expect the same amount of care in all code. But I do think that users should not have to know which units are required. Especially given that the requirements are not document anywhere.
| def _pixel_coordinate_in_detector_system( | ||
| pixelID: sc.Variable, | ||
| ) -> tuple[sc.Variable, sc.Variable]: | ||
| def pixel_coordinates_in_detector_system() -> tuple[sc.Variable, sc.Variable]: |
There was a problem hiding this comment.
The return annotation is incorrect now. And I am not sure the function name still makes sense either.
There was a problem hiding this comment.
👍 I'll fix the annotation. What name would you have?
src/ess/amor/load.py
Outdated
| fp: Filename[RunType], | ||
| ) -> AngleCenterOfIncomingToHorizon[RunType]: | ||
| (kad,) = load_nx(fp, "NXentry/NXinstrument/master_parameters/kad") | ||
| return sc.scalar(kad['value'].data['dim_1', 0]['time', 0].value, unit='deg') |
There was a problem hiding this comment.
Why does it access the 0th element in dim_1 and time? What are the other elements?
There was a problem hiding this comment.
It's stored as an NXLog and when we read that it ends up this shape, but it doesn't make sense for this to be anything other than a scalar, so we only read the first value.
There was a problem hiding this comment.
Ok. Can you check that there is indeed a single entry? Otherwise, it may silently mess up when there is a time-dependent log. You could just use squeeze.
There was a problem hiding this comment.
It's not unusual that there is more than one entry. But it has been determined that it's correct to read the first entry. See the other parameters that are read the same way, sample_rotation and detector_rotation.
I don't think it makes sense to complicate the processing or raise errors just because there are more than a single entry here.
There was a problem hiding this comment.
In that case, please add a comment that explains the choice and why it is ok to use the first element.
src/ess/amor/conversions.py
Outdated
| "divergence_angle": "pixel_divergence_angle", | ||
| "wavelength": wavelength, | ||
| "theta": theta, | ||
| "angle_of_divergence": angle_of_divergence, |
There was a problem hiding this comment.
What is the difference of 'divergence_angle' and 'angle_of_divergence'? If they are different, can you rename one? The two are synonyms in English.
There was a problem hiding this comment.
That's a very good point. It's confusing. divergence_angle is the divergence angle on the detector, angle_of_divergence is the divergence angle of the incoming beam on the sample.
Without gravity they are identical.
I'll try to come up with better names.
I understand it's a lot, the changes are large. |
|
I'll close this and split it up into smaller PRs. |
Co-authored-by: Jan-Lukas Wynen <[email protected]>
Co-authored-by: Jan-Lukas Wynen <[email protected]>
There was a problem hiding this comment.
Unused AngleCenterOfIncomingToHorizon in this file? Where does the correction happen?
There was a problem hiding this comment.
The correction happens in normalize.py/reduce_reference.
Not sure what you mean by AngleCenterOfIncomingToHorizon is unused? Do you mean that it should be used, or do you mean that it's imported but unused? I can't see that it occurs in the file.
There was a problem hiding this comment.
Do you mean that it should be used, or do you mean that it's imported but unused?
Yes, this is what I meant. It should have been caught by ruff, I think. But it looks like the setup in ESSreflectometry is rather old.
There was a problem hiding this comment.
AngleCenterOfIncomingToHorizon does not occur in the file src/ess/reflectometry/supermirror/__init__.py
jokasimr
force-pushed
the
simplify-amor
branch
from
84e0a55 to
afd5ce9
Compare
src/ess/amor/load.py
Outdated
| data.coords["chopper_frequency"] = chopper_frequency | ||
| data.coords["chopper_separation"] = sc.abs(chopper_separation) | ||
| data.coords["L1"] = sc.abs(chopper_distance) | ||
| data.coords["L2"] = data.coords['distance_in_detector'] + Detector.distance |
There was a problem hiding this comment.
As a reminder, this should be done by adapting the bemaline graph for coord transforms from ScippNeutron.
There was a problem hiding this comment.
The beamline graph defines the following quantities:
- incident_beam is the vector of incoming neutrons on the sample
- L1
- scattered_beam is the vector of neutrons that were scattered off the sample
- L2
- Ltotal
But some of those are not used in this workflow, and some are to be replaced anyway:
- incident_beam - this quantity is not needed in the workflow
- L1 - should be replaced
- scattered_beam - this quantity is not used in the workflow
- L2 - should be replaced
- Ltotal - this is used in the workflow
In the end, the only useful quantity we get from the graph is Ltotal. But that is neither difficult to compute or likely to change.
For this reason I don't think it's necessary to use that graph here. Adapting it to fit here is likely more code than saved by using it.
I think it's good to take into consideration that it becomes harder for a user to follow what transformations were applied if they have to go and look in a different package to find the definitions. And while that's sometimes unavoidable, it doesn't seem like it is here.
There was a problem hiding this comment.
Take a look at the LoKI workflow: https://scipp.github.io/esssans/user-guide/loki/loki-iofq.html It has a provider for its coord transform graphs. This way, we can have the same code do the actual transforms for multiple instruments. We only need to replace the graphs.
| @@ -54,7 +54,7 @@ def mask_events_where_supermirror_does_not_cover( | |||
| theta( | |||
| sam.bins.coords["wavelength"], | |||
| sam.coords["pixel_divergence_angle"], | |||
| ref.coords["L2"], | |||
| sam.coords["L2"], | |||
There was a problem hiding this comment.
Is this a bugfix? If yes, please add a test!
There was a problem hiding this comment.
Not a bugfix, they are equivalent, I just found the second version more clear.
Co-authored-by: Jan-Lukas Wynen <[email protected]>
| sam.bins.coords["wavelength"], | ||
| sam.coords["pixel_divergence_angle"], | ||
| sam.coords["L2"], | ||
| ref.coords["sample_rotation"], |
There was a problem hiding this comment.
It would be good to use the AmorCoordinates transformation graph to compute supermirror reflectivity here, because as it is now it looks like the coordinate transformations can be changed everywhere by replacing AmorCoordinates but they can't because here they are hardcoded.
However, I wasn't able to use transform_coords here in a way that doesn't rely on knowing the input arguments to the coordinate transformation functions in the graph.
If you have any suggestions about how to do this, please tell me.
There was a problem hiding this comment.
What do you mean by
However, I wasn't able to use transform_coords here in a way that doesn't rely on knowing the input arguments to the coordinate transformation functions in the graph.
What about the inputs would you have to know that you don't have to know otherwise?
There was a problem hiding this comment.
Sorry that was unclear. You have to know it either way. It's just that since you have to know it either way we don't get any benefit from using transform_coords here and then I thought this version was more explicit and easier to understand.
| ) | ||
| out = sc.asin(out, out=out) | ||
| out -= sample_rotation | ||
| out -= sample_rotation.to(unit='rad') |
There was a problem hiding this comment.
Are the variables you are converting all scalar? For arrays, I would use copy=False.
src/ess/amor/types.py
Outdated
| @@ -9,6 +9,8 @@ | |||
| AngularResolution = NewType("AngularResolution", sc.Variable) | |||
| SampleSizeResolution = NewType("SampleSizeResolution", sc.Variable) | |||
|
|
|||
| AmorCoordinates = NewType("AmorCoordinates", dict) | |||
There was a problem hiding this comment.
| AmorCoordinates = NewType("AmorCoordinates", dict) | |
| CoordTransformationGraph = NewType("CoordTransformationGraph", dict) |
to match the name in ESSsans.
Why did you include the word 'Amor' in the name?
There was a problem hiding this comment.
Because it's specific to the Amor beamline
There was a problem hiding this comment.
But you didn't prefix the other domain types with 'Amor'. The concrete value is specific, but the concept of a coord transform graph is generic.
| sam.bins.coords["wavelength"], | ||
| sam.coords["pixel_divergence_angle"], | ||
| sam.coords["L2"], | ||
| ref.coords["sample_rotation"], |
There was a problem hiding this comment.
What do you mean by
However, I wasn't able to use transform_coords here in a way that doesn't rely on knowing the input arguments to the coordinate transformation functions in the graph.
What about the inputs would you have to know that you don't have to know otherwise?
A while back I realized the reflectometry workflow can be made a lot simpler than it is now.
I think it was worthwhile to do this refactoring even if Amor is not an ESS instrument because the data reduction is very similar to Estia, and this new version serves as a better template for the Estia data reduction.
There are some things to fix here that I am aware of:
ReducibleDataReducedReferenceSampleReference) Name suggestions are appreciated!Documentation for some of the new coordinate transformations.Done. This version does not use the coordinate transformations in scippneutron, mainly because most of them just are not in scippneutron and are only used here, and partly because it was just simpler to keep the transformations directly in the package.