Frame unwrapping and WFM #472

SimonHeybrock · 2023-11-30T08:45:05Z

This represents the latest approach to frame unwrapping and WFM (see old #415).

This clearly makes a bunch of assumptions which may turn out to be wrong in practice. At this point I have little confidence that this will be anything close to the final solution (or even that there is a solution at all — in the end we might still be forced to fall back to defining bounds more manually).

Testing is obviously challenging, and given that it is still not clear how all the choppers will be handled, I believe we must iterate on the whole solution (including what @jl-wynen has kicked of with #466) in the coming months.

Notes:

The old tof.frames was removed.
Adds tof.fakes, which I hope we can develop further for these and similar tests. This is currently my only idea on how to move forward: Build more realistic fakes, fix the chopper and unwrapping code, iterate.
Unlike the old unwrapping code, this can also handle event_time_zero. This will be essential for experiments with high time-resolution requirements.
Relies on the recently added scippneutron.tof.chopper_cascade module
My initial benchmarks indicate that using scipp.lookup over scipp.where and scipp.bin during unwrapping and WFM stitching (respectively) has equivalent or better performance. I think it is also clearer.
I added Sciline to the docs builds (so we can show the graphs), but not as a Scippneutron requirement. That is, we will provide providers (:grimacing: ) and domain types, but not more.

Fixes scipp/ess#165.

Follow-up work:

nvaytet · 2023-11-30T15:08:29Z

I had a short look but I think it would be much easier for me if we looked at the graphs together and you can explain the different steps, as well as the differences between the different graphs.

nvaytet · 2023-12-08T09:23:35Z

src/scippneutron/tof/unwrap.py

+    # Find last frame before ltotal
+    frame_before_detector = None
+    for frame in frames:
+        if frame.distance > ltotal:


Just checking: I guess we are guaranteed that the frames are ordered according to distance?

Yes, chopper_toolkit sorts the choppers before applying them.

nvaytet · 2023-12-08T09:27:02Z

src/scippneutron/tof/unwrap.py

-    choppers: Choppers, source_chopper_name: SourceChopperName
+    choppers: Choppers,
+    source_time_range: SourceTimeRange,
+    source_chopper_name: Optional[SourceChopperName],


Should source_chopper_name have a default None value?

Sciline does not use default values.

nvaytet · 2023-12-08T09:27:17Z

src/scippneutron/tof/unwrap.py

+    are used to define a fake chopper.
+    """
+    if source_chopper_name is not None:
+        return choppers[source_chopper_name]


So if the source_chopper_name exists, then source_time_range is unused. Should that also be an optional input?

I don't think so, because then the pipeline would fail at runtime if neither is given?

Say you know the source chopper name, but not the source time range, you'd still have to supply a source time range as input.
If you don't know it, you could always make it to be nonsense because it won't be used, but it seems a little strange to require it?

Without a source time range the whole chopper-cascade math won't work. This is required anyway for other parts of the workflow.

Should we then have a check to make sure that the source_time_range is consistent with the source chopper if source_chopper_name is given?

How would that be possible?

Couldn't we check that the chopper open and close times match the source_time_range? Or maybe I misunderstood something and in some part of the pipeline you need the source chopper, and in others you need the source_time_range and they can in principle be different? (e.g. if the source chopper is not at a distance of 0.0?)

The "source chopper" may define a different time range, i.e., it is a pulse-shaping chopper to increase the resolution. It will typically have a tighter time range.

SimonHeybrock · 2024-01-31T14:34:50Z

src/scippneutron/tof/unwrap.py

+
+
+def offset_to_time_of_flight_wfm(
+    time_offset: OffsetFromWrapped,


I think this may actually be wrong. I think we need the unwrapped time offset, not just the offset from wrapped to perform the WFM TOF offsetting.

Converting to draft until verified.

jokasimr · 2024-02-06T15:33:04Z

src/scippneutron/tof/unwrap.py

+        wavelength_max=source_wavelength_range[-1],
+    )
+    frames = frames.chop(choppers.values())
+    return FrameAtDetector(frames[ltotal])


Doesn't this assume:

that the neutrons move in straight lines to the detector elements?

that there is no inelastic scattering?

Correct. See my comment in the PR:

This clearly makes a bunch of assumptions which may turn out to be wrong in practice. At this point I have little confidence that this will be anything close to the final solution (or even that there is a solution at all — in the end we might still be forced to fall back to defining bounds more manually).

SimonHeybrock · 2024-02-12T06:21:44Z

@nvaytet @jl-wynen @jokasimr I'd like some input here, on the latest changes.

jl-wynen

I don't have anything enlightening to say. This just looks like what we talked about.

jl-wynen · 2024-02-12T08:17:38Z

docs/api-reference/frame-unwrapping.ipynb

Typo in l73: "For, e.g., detector 1 when then observe:\n", should be "For, e.g., detector 1 we then observe:\n",. And 'then' can be left out, too.

jl-wynen · 2024-02-12T08:20:39Z

docs/api-reference/frame-unwrapping.ipynb

-    "   $$\n",
-    "   \n",
-    "Note that there are other valid definitions of $t_\\text{pivot}$, differing only in how neutrons in the intermediate \"background\" region (region between two pulses) are mapped to frames."
+    "As illustrated in the figure, the pivot time $t_\\text{pivot}$ depends on the detector or rather the distance of the detector (or monitor) from the scattering position.\n"


In the diagram, the pivot for detector 1 is at the black dashed line, i.e., the beginning of the orange frame. But in detector 2 it is in the middle between the frames.

Not sure what you are describing. Doesn't it look like in the docs for you?

docs/api-reference/frame-unwrapping.ipynb

jl-wynen · 2024-02-12T08:32:56Z

docs/api-reference/frame-unwrapping.ipynb

-   "metadata": {},
-   "source": [
-    "At the sample or a detector we would only see this for a particular `Ltotal`:"
+    "Choppers may be used to skip pulses, for the purpose of a simultaneous study of a wider wavelength range.\n",


What do you mean by simultaneous here?

You can use every pulse, measure 1-4 Angstrom, change the choppers, measure 4-8 Angstrom, and combine the results. That would be serial, as opposed to simultaneous (I think this is the term I have heard scientists use).

src/scippneutron/tof/unwrap.py

nvaytet

A couple of questions.

Also, maybe I've asked this before, but why is this notebook in the API reference? It makes it pretty difficult to find in the docs I feel.

nvaytet · 2024-02-12T10:06:25Z

src/scippneutron/tof/unwrap.py

@@ -378,7 +388,9 @@ def time_of_flight_origin_from_chopper(
    source_time_open = source_chopper.time_open[0]
    source_time_close = source_chopper.time_close[0]
    time_zero = 0.5 * (source_time_open + source_time_close)
-    return TimeOfFlightOrigin(time_zero)
+    return TimeOfFlightOrigin(
+        time=time_zero, L1=sample_distance - source_chopper.distance


Is the TimeOfFlightOrigin ill-defined if there is no sample in the run, or if you are trying to unwrap times for a monitor?

It does not matter if there is no sample, this is about the definition of L1. This (or sample_position) has to be defined even if there is no "sample", since all the coord transforms in ScippNeutron rely on it.

But you have actually stumbled upon a related problem: If we are unwrapping a monitor then the coord transforms defined in ScippNeutron do not use L1 but Ltotal, so we should be correcting that.

@nvaytet Please have a look at the latest commit, this should fix the aforementioned issue.

👍 Remember to add a test that checks for this when you get round to adding tests.

nvaytet · 2024-02-12T10:08:02Z

src/scippneutron/tof/unwrap.py

-    return TimeOfFlightOrigin(sc.DataArray(neg_shift, coords={'subframe': times}))
+    return TimeOfFlightOrigin(
+        time=sc.DataArray(neg_shift, coords={'subframe': times}),
+        L1=sample_distance - source_chopper.distance,


Should they be distances or should they be defined as vector positions in space and take the norm of the difference?

The chopper_cascade using a simple distance, so this is what we have to use here. This will hopefully allow for hiding details such as curved guides (we do not want to chopper_cascade toolkit to deal with such complications, if possible).

docs/api-reference/frame-unwrapping.ipynb

jokasimr

I still don't really understand the WFM case and why it's necessary to treat it separately. My understanding was that given any chopper configuration we can propagate the source pulse through the cascade and compute the wavelength-time diagram at every detector(element), and from it + the event time data we can compute the time of flight. From that perspective the WFM case is just one particular chopper configuration and it's not clear to me why it needs a separate code path.

But this shouldn't be an obstacle to merge this, just a heads up that this is probably something we should talk about again.

jokasimr · 2024-02-12T10:57:53Z

src/scippneutron/tof/unwrap.py

+TimeOffset = NewType('TimeOffset', sc.Variable)
+"""Unwrapped time offset relative to the pulse time."""
+
+DeltaFromWrapped = NewType('DeltaFromWrapped', sc.Variable)


Maybe add a docstring to describe what DeltaFromWrapped represents

Did you review an outdated diff?

jokasimr · 2024-02-12T12:29:05Z

src/scippneutron/tof/unwrap.py

@@ -311,7 +314,7 @@ def offset_from_wrapped(
        + sc.concat([frame_period, sc.zeros_like(frame_period)], dim),
        coords={dim: time.to(unit=elem_unit(wrapped_time_offset))},
    )
-    return OffsetFromWrapped(sc.lookup(offset, dim=dim)[wrapped_time_offset])
+    return DeltaFromWrapped(sc.lookup(offset, dim=dim)[wrapped_time_offset])


Would it be possible to use interpolate here? If so, is there a reason why lookup is advantageous

scipp.lookup can handle event data, interpolate cannot, and may be slower.

jokasimr · 2024-02-12T12:37:09Z

src/scippneutron/tof/unwrap.py

-    time_offset: OffsetFromWrapped,
-    source_chopper: SourceChopper,
-) -> OffsetFromTimeOfFlight:
+def time_of_flight_origin(source_chopper: SourceChopper) -> TimeOfFlightOrigin:


Why do the user have to define the SourceChopper? Can't that be automatically inferred from the description of the chopper setup?

No, not really, it is not necessarily the first chopper, and it might not even be an actual chopper, but a "virtual" one that one might define (e.g., to use the source pulse as "source chopper").

Co-authored-by: Jan-Lukas Wynen <[email protected]> Co-authored-by: Neil Vaytet <[email protected]>

Co-authored-by: Jan-Lukas Wynen <[email protected]>

SimonHeybrock · 2024-02-13T05:28:05Z

Also, maybe I've asked this before, but why is this notebook in the API reference? It makes it pretty difficult to find in the docs I feel.

It used to be more visible, but was moved there during the switch to the new theme. @jl-wynen can comment on why it is here or where we can move it.

SimonHeybrock · 2024-02-13T05:37:42Z

From that perspective the WFM case is just one particular chopper configuration and it's not clear to me why it needs a separate code path.

It does not a different code part, for the most part. But the time-of-flight origin needs to be defined differently, in particular different neutrons within a frame have a different origin.

Setting L1 would have been broken for monitors. Furthermore, it might have been ignored if Ltotal was computed previously.

nvaytet · 2024-02-13T13:11:36Z

Green light from me for this round of review.

Do I remember correctly that you said you were going to add tests after this round of review?

SimonHeybrock · 2024-02-20T09:04:57Z

@nvaytet Tests updated!

nvaytet · 2024-02-23T15:26:28Z

src/scippneutron/tof/unwrap.py

@@ -505,14 +508,21 @@ def to_time_of_flight(
    if da.bins is not None:
        da = da.copy(deep=False)
        da.data = sc.bins(**da.bins.constituents)
-        da.bins.coords['tof'] = time_offset + delta
-        da.bins.coords['time_zero'] = da.bins.coords['pulse_time'] - delta.to(
+        da.bins.coords['tof'] = time_offset - delta


SimonHeybrock added 8 commits November 27, 2023 15:22

Experiment with new frame-unwrapping workflow

60ecf33

Simpler WFM frame offsetting

243729c

Move code to module

ce8ef99

Cleanup

6179931

Add provider setup helper

c8dcebb

Add docstrings

69f7078

Start graph refactor so we can also offset event_time_zero

218c45e

Compute OffsetToTimeOfFlight

9522b82

SimonHeybrock requested a review from nvaytet November 30, 2023 08:45

SimonHeybrock mentioned this pull request Dec 6, 2023

Implement frame unwrapping (and pulse skipping) into an actual utility #402

Closed

SimonHeybrock added 5 commits December 7, 2023 10:56

Support beamlines using full pulse

712244c

Document

6da49f3

Handle inter-subframe events using NAN

99c0478

Generalize graph so we can handle monitors

23b9306

Sketch approach for non-event data

13c908f

nvaytet reviewed Dec 8, 2023

View reviewed changes

SimonHeybrock added 13 commits December 12, 2023 12:51

Fix subframe merge in chopper_cascade module

266141e

WIP tests [no ci]

c1c5b7f

Support getting frame at any distance

f5bcb93

Basic test of unwrapping code

172b1c9

Use choppers from fakes for docs

6be1439

Begin updating frame-unwrapping notebook

42f11bb

Remove merged notebook

ae1f4f5

Remove old tof.frames module

8016270

Reject overlapping frames

ac58634

Module docstring

aeda208

Test frame period

20f51e3

More tests

68fca82

Test frame skipping

1da5272

Extend docstrings

62158fd

SimonHeybrock requested a review from nvaytet January 24, 2024 14:04

Merge branch 'main' into unwrap

938486e

SimonHeybrock commented Jan 31, 2024

View reviewed changes

SimonHeybrock marked this pull request as draft January 31, 2024 14:35

jokasimr reviewed Feb 6, 2024

View reviewed changes

SimonHeybrock added 3 commits February 7, 2024 09:25

Fix WFM offset computation and restructure logic

96e5de0

Improve graph structure

5d60e14

Add L1 and improve docs

bb79688

jl-wynen reviewed Feb 12, 2024

View reviewed changes

nvaytet reviewed Feb 12, 2024

View reviewed changes

docs/api-reference/frame-unwrapping.ipynb Outdated Show resolved Hide resolved

jokasimr reviewed Feb 12, 2024

View reviewed changes

SimonHeybrock and others added 3 commits February 13, 2024 06:26

Apply suggestions from code review

ec8b9fe

Co-authored-by: Jan-Lukas Wynen <[email protected]> Co-authored-by: Neil Vaytet <[email protected]>

Update src/scippneutron/tof/unwrap.py

c5ccebe

Co-authored-by: Jan-Lukas Wynen <[email protected]>

Apply automatic formatting

c67e32c

SimonHeybrock added 2 commits February 13, 2024 06:42

Grammar

efd7ed4

Correct Ltotal instead of setting L1

f41165e

Setting L1 would have been broken for monitors. Furthermore, it might have been ignored if Ltotal was computed previously.

jl-wynen mentioned this pull request Feb 14, 2024

[Requirement] WFM stitching scipp/essdiffraction#26

Open

SimonHeybrock added 2 commits February 20, 2024 09:19

Fix some tests

91a84f5

Test to_time_of_flight

d24ac9f

SimonHeybrock marked this pull request as ready for review February 20, 2024 09:04

nvaytet approved these changes Feb 23, 2024

View reviewed changes

SimonHeybrock merged commit 22ee6eb into main Feb 26, 2024
6 checks passed

SimonHeybrock deleted the unwrap branch February 26, 2024 03:48



		def offset_to_time_of_flight_wfm(
		time_offset: OffsetFromWrapped,

Frame unwrapping and WFM #472

Frame unwrapping and WFM #472

Conversation

SimonHeybrock commented Nov 30, 2023 • edited Loading

nvaytet commented Nov 30, 2023

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

SimonHeybrock Jan 31, 2024 • edited Loading

Choose a reason for hiding this comment

jokasimr Feb 6, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SimonHeybrock commented Feb 12, 2024

jl-wynen left a 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

nvaytet left a 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 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jokasimr Feb 12, 2024 • edited Loading

Choose a reason for hiding this comment

SimonHeybrock Feb 13, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SimonHeybrock commented Feb 13, 2024

SimonHeybrock commented Feb 13, 2024

nvaytet commented Feb 13, 2024

SimonHeybrock commented Feb 20, 2024

Choose a reason for hiding this comment

SimonHeybrock commented Nov 30, 2023 •

edited

Loading

SimonHeybrock Jan 31, 2024 •

edited

Loading

jokasimr Feb 6, 2024 •

edited

Loading

jokasimr Feb 12, 2024 •

edited

Loading

SimonHeybrock Feb 13, 2024 •

edited

Loading