Skip to content

[essreduce] On-the-fly wavelength lookup table in generic workflow#603

Merged
nvaytet merged 66 commits into
mainfrom
lut-in-generic-wf
Jun 3, 2026
Merged

[essreduce] On-the-fly wavelength lookup table in generic workflow#603
nvaytet merged 66 commits into
mainfrom
lut-in-generic-wf

Conversation

@nvaytet
Copy link
Copy Markdown
Member

@nvaytet nvaytet commented May 22, 2026

Insert the providers from the LUT workflow into the GenericUnwrapWorkflow to compute the lut on the fly with chopper cascade.
Breaking change: we compute a LUT per component (detector bank, monitor) instead of one table for all.

Building small tables is faster, and interpolating inside a small table is also faster (most probably due to cache misses in larger tables).

By default, the workflow tries to load the lookup table from file, which behaves like the old code.
The workflow now has 3 modes it can operate in:

  • wavelength_from='file': the default = read pre-computed table from file
  • wavelength_from='simulation': on-the-fly computed table using tof simulation from chopper params (slow)
  • wavelength_from='analytical': on-the-fly computed table using chopper_cascade -> what we want to default to in the future.

Suggestions for a better name than 'analytical' are welcome.

The idea is that each package would then move over to analytical mode when they wish to do so (in follow-up PRs).


Additional notes:

I tried to reduce the number of parameters we need to set on the workflow:

  • LtotalRange is now determined from the pixel positions of the detector bank, or the monitor position.
  • PulseStride is now guessed from the chopper frequencies (see comment here)

Finally, the LookupTableWorkflow that was used to create lookup tables using tof is now just an alias for the GenericUnwrapWorkflow that uses tof by default, so that the behaviour remains the same as before. We plan to remove the alias in the future.

@nvaytet nvaytet added the essreduce Issues for essreduce. label May 22, 2026
Comment on lines +353 to 356
"lut_wf[LtotalRange[SampleRun, snx.NXdetector]] = (\n",
" sc.scalar(25.0, unit=\"m\"),\n",
" sc.scalar(80.0, unit=\"m\"),\n",
")\n",
Copy link
Copy Markdown
Member

@SimonHeybrock SimonHeybrock Jun 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we chatted about this, but remind me why this needs to be set. It seems the range is not per-component, so we is it not computed?

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Jun 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It no longer needs to be set.

It is however useful to keep as an entry point in the workflow in case we wish to compute a lookup table that covers a wide range of distances with an expensive calculation using tof or McStas.

In the present example, instead I will set the DetectorLtotal (which is defined further up in the notebooks).

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that my plan is to probably remove the dream.ipynb and wfm.ipynb notebooks here; the new analytical-unwrap.ipynb is meant to replace them.


@dataclass
class SimulationResults:
class SimulationResultsBaseClass:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do you need both instead of making the dataclass generic?

Comment on lines +131 to +134
Note also that the range of total flight paths is supplied manually to the workflow
instead of being read from the input data, as it allows us to compute the expensive
part of the workflow in advance (the lookup table) and does not need to be repeated
for each run, or for new data coming in in the case of live data collection.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this answers my comment above... partially. But live data collection caches static nodes (using StreamProcessor), so this is not a reason, I think?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But you still need to repeat this for every RunType which usually have the same detector positions. Or does the workflow still need to do that anyway?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We cannot have two runs at the same time in live data. If we need to subtract vanadium etc. then that must come from a previously stored intermediate result.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment is actually a remnant of when we were only using tables from disk or computed using tof.
I think I will just remove the comment.

- table.coords["distance"][0],
time_resolution=table.coords["event_time_offset"][1]
- table.coords["event_time_offset"][0],
# TODO: Do we still want to store the chopper info in the lookup table?
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Planned to resolve here or later?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Later

DistanceResolution: sc.scalar(0.1, unit="m"),
TimeResolution: sc.scalar(250.0, unit='us'),
SourceBounds: SourceBounds(
time=(sc.scalar(0.0, unit='ms'), sc.scalar(5.0, unit='ms')),
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where does 5.0 come from?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a comment.

Comment on lines -572 to +574
for compress_mode in (Compression.GZIP, Compression.BITSHUFFLE_LZ4):
for compress_mode in (Compression.NONE, Compression.BITSHUFFLE_LZ4):
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did this one removed?

Comment thread packages/essdiffraction/src/ess/dream/workflows.py Outdated
analytical calculations to propagate and chop a pulse through the chopper
cascade and build the lookup table. The 'simulation' mode uses ``tof`` to trace
individual neutrons through the chopper system and build the table.
The 'file' mode loads a pre-computed table from a file.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of explaining this here, can you link to a docs page that explains the 3 modes in detail?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok I will do that.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See lut-building-methods.ipynb.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

return lut_wf.compute(unwrap.LookupTable[AnyRun, NXdetector])


def test_pipeline_can_compute_dspacing_result_using_custom_built_tof_lookup(
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now we should also have tests with a table that is computed on the fly. Is it easy to parametrize a test?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I parametrized most tests for essreduce.unwrap.
Do you want me to parametrize all such tests in the child packages?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is some custom code to handle lut filenames, so ideally, they should be parametrized. But that would lead to additional reference data because the workflows don't produce exactly the same result for each lut, right?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I parametrized the dream geant4 tests.
I suggest we update other tests once packages individually move over to using the analytical mode?

Comment on lines +31 to +38
for run_type in (SampleRun, OpenBeamRun):
wf[LookupTableFilename[run_type, NXdetector]] = (
odin.data.odin_wavelength_lookup_table()
)
# Cache the lookup table
wf[LookupTable[run_type, NXdetector]] = wf.compute(
LookupTable[run_type, NXdetector]
)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this possible? I.e., can we set a lut for all runs without looping through them?

Suggested change
for run_type in (SampleRun, OpenBeamRun):
wf[LookupTableFilename[run_type, NXdetector]] = (
odin.data.odin_wavelength_lookup_table()
)
# Cache the lookup table
wf[LookupTable[run_type, NXdetector]] = wf.compute(
LookupTable[run_type, NXdetector]
)
wf[LookupTableFilename[RunType, NXdetector]] = (
odin.data.odin_wavelength_lookup_table()
)
# Cache the lookup table
wf[LookupTable[RunType, NXdetector]] = wf.compute(
LookupTable[RunType, NXdetector]
)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this is the feature of sciline I stumbled upon this morning that I didn't know existed.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm doesn't seem to work, I get UnsatisfiedRequirement when I try to compute LookupTable[RunType, NXdetector]

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, that can't work. Maybe

    wf[LookupTable[RunType, NXdetector]] = wf.compute(
        LookupTable[SampleRun, NXdetector]
    )

if you are sure the tables are all the same.

But this is only in a test. I was after a general approach to simplify code. In real workflows

    wf[LookupTableFilename[RunType, NXdetector]] = (
        odin.data.odin_wavelength_lookup_table()
    )

should be all you need. And that should work.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe even

wf[LookupTableFilename] = odin.data.odin_wavelength_lookup_table()

?

Comment thread packages/essreduce/docs/user-guide/unwrap/frame-unwrapping.ipynb
Comment thread packages/essreduce/src/ess/reduce/unwrap/types.py Outdated
Comment on lines +131 to +134
Note also that the range of total flight paths is supplied manually to the workflow
instead of being read from the input data, as it allows us to compute the expensive
part of the workflow in advance (the lookup table) and does not need to be repeated
for each run, or for new data coming in in the case of live data collection.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But you still need to repeat this for every RunType which usually have the same detector positions. Or does the workflow still need to do that anyway?

def make_wavelength_lookup_table(
simulation: SimulationResults,
ltotal_range: LtotalRange,
def make_wavelength_lut_from_simulation(
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about 'monte_carlo' instead of 'simulation' to be more precise? Not sure what to do about 'analytical', though.

@SimonHeybrock
Copy link
Copy Markdown
Member

Found a crash in the new analytical mode while integrating it into ESSlivedata: make_wavelength_lut_from_polygons raises IndexError on a degenerate subframe when a chopper has a positive frequency.

File ".../ess/reduce/unwrap/lut.py", line 554, in _polygon_intersections
    if b[0, 0] == b[1, 0]:
                  ~^^^^^^
IndexError: index 1 is out of bounds for axis 0 with size 1

A polygon bound decomposes to a single vertex (the subframe collapses), so indexing b[1] is out of range.

Minimal reproducer (essreduce only)

import scipp as sc
import scippnexus as snx
from scippneutron.chopper import DiskChopper
from ess.reduce import unwrap
from ess.reduce.unwrap import GenericUnwrapWorkflow
from ess.reduce.nexus.types import AnyRun, Position


def chopper(z, freq):
    return DiskChopper(
        axle_position=sc.vector([0, 0, z], unit='m'),
        frequency=sc.scalar(freq, unit='Hz'),
        beam_position=sc.scalar(0.0, unit='deg'),
        phase=sc.scalar(0.0, unit='rad'),
        slit_begin=sc.array(dims=['cutout'], values=[0.0], unit='deg'),
        slit_end=sc.array(dims=['cutout'], values=[90.0], unit='deg'),
        slit_height=None,
        radius=sc.scalar(0.35, unit='m'),
    )


def run(freq):
    wf = GenericUnwrapWorkflow(
        run_types=[AnyRun], monitor_types=[], wavelength_from='analytical'
    )
    wf[unwrap.DiskChoppers[AnyRun]] = sc.DataGroup(
        {'chopper1': chopper(-15.0, freq), 'chopper2': chopper(-13.0, freq)}
    )
    wf[Position[snx.NXsource, AnyRun]] = sc.vector([0, 0, -25.0], unit='m')
    wf[unwrap.PulseStride[AnyRun]] = 1
    wf[unwrap.LtotalRange[AnyRun, snx.NXdetector]] = (
        sc.scalar(5.0, unit='m'),
        sc.scalar(30.0, unit='m'),
    )
    wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit='m')
    wf[unwrap.TimeResolution] = sc.scalar(250.0, unit='us')
    return wf.compute(unwrap.LookupTable[AnyRun, snx.NXdetector])


run(-14.0)  # OK
run(+14.0)  # IndexError

Trigger

The sign of frequency is what flips it: -14 Hz builds the table fine, +14 Hz crashes. Everything else identical. The simulation mode tolerated this geometry (it just yields sparse neutrons), so it only surfaced on the analytical path.

Real beamline choppers use the signed (negative, anti-clockwise) convention, so this does not bite in practice — our LOKI geometry computes fine. We hit it only with a synthetic test fixture that fed +14 Hz. Reporting because a positive frequency arguably shouldn't produce an unguarded IndexError; a degenerate-subframe guard or an explicit error would be friendlier. Happy to leave the fix to you — on our side we'll just fix the fixture's sign.

SimonHeybrock added a commit to scipp/esslivedata that referenced this pull request Jun 2, 2026
Build the LUT pipeline via GenericUnwrapWorkflow(wavelength_from='analytical')
instead of the simulation-based LookupTableWorkflow, following the essreduce
rewrite that moved the wavelength lookup table into the generic workflow.
Analytical mode derives the table from chopper-cascade polygon geometry: far
faster than the tof neutron simulation and needs no simulated-neutron count.

- SourcePosition -> Position[NXsource, AnyRun]; LtotalRange / PulseStride /
  LookupTable re-parametrised with [AnyRun, NXdetector].
- Drop NumberOfSimulatedNeutrons and the Simulation params model.
- Synthetic chopper test fixture uses negative (anti-clockwise) chopper
  frequencies; a positive frequency drives essreduce's analytical chopper
  cascade into a degenerate frame (reported on scipp/ess#603).

Depends on an unreleased ess.reduce: the analytical wavelength-LUT API lives on
the essreduce lut-in-generic-wf branch and is absent from the latest release
(26.5.0). The essreduce dependency pin must be bumped once that work is
released before this can merge.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@nvaytet
Copy link
Copy Markdown
Member Author

nvaytet commented Jun 2, 2026

I added

if len(b) < 2:
    continue

inside the loop. Does that fix it for you?

LookupTable,
LookupTableFilename,
Lut,
# Lut,
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove?

Comment on lines +550 to +553
if len(b) < 2:
# This can happen if the polygon is degenerate (collapsed to a single
# vertex).
continue
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you implying that the code below that also uses b and the bounds works nevertheless?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it does?
But maybe you think it's safer to just discard the polygon if it has collapsed?

Can you clarify what a collapsed polygon means? Should there just be no polygon at all, or is there a single infinitessimally small point where neutrons could get through? (maybe it amounts to the same anyway)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question. Does scippneutron keep around such degenerate polygons? Should they be dropped?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From private conversation:

So I looked at the subframes in your exampl (Minimal reproducer).
the subframe that is still there is not just a single vertex, it's made up of 3 vertices and they are all the same: [0, 0, 0] in both time and wavelength.

It's a peculiarity that happens because we allow a 0 wavelength (the default bounds of the pulse) which is infinite speed.
The choppers are open just before t=0 and close exactly at t=0.
I'm thinking that the neutrons with infinite speed can make it through just at t=0, but neutrons with infinite speed that started at t=1ms don't make it through because they were born too late.
So we are left with a single point at the corner of the pulse that can make it through.

If I change the pulse bounds to start at 0.001 angstrom the error goes away.
I get a new error which crashes when the subframes are empty, but I also fixed that now.

@nvaytet nvaytet added this pull request to the merge queue Jun 3, 2026
Merged via the queue into main with commit ee27de7 Jun 3, 2026
24 checks passed
@nvaytet nvaytet deleted the lut-in-generic-wf branch June 3, 2026 14:19
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

essreduce Issues for essreduce.

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants