Essreduce + Essimaging: tofectomy

[nvaytet]

Replacing time-of-flight with wavelength in the GenericTofWorkflow, which is now called GenericUnwrapWorkflow.

Note: once this is merged, every package using the GenericTofWorkflow not in the monorepo will break...

[jokasimr]

About the name, what about WavelengthWorkflow or EstimateWavelengthWorkflow?

[jokasimr]

What about WavelenghtLookupTable? It's a bit more specific.

I'm sure the technique packages will have other lookup tables so we will have to rename this type in the technique packages anyway. Is it easier to directly give it a specific name here instead?

[jokasimr]

After thinking a bit more, LookupTable is probably better in this module. We can rename it in the instrument packages.

[nvaytet]

Yes, I think that when you take the module name into account, unwrap.Lookuptable becomes specific enough.

[jokasimr]

Would it make sense to name the module wavelength instead of unwrap?

I get that wavelength is such a common variable name, that there will likely be clashes if the module is imported like from ess.reduce import wavelength but I think it is still better than unwrap because that term is a bit overloaded.

[nvaytet]

The name unwrap was suggested by both JL & Sun, and I also liked it, so I'd like to go with that.
We had an issue with time_of_flight which turned out to be too specific because when I changed things to work with wavelengths, I had to rename everywhere.

If we decide in the future to use something else than wavelength (I don't think we will but we never know), we won't have to rename, because unwrapping is the essence of what we are doing here, no matter what underlying quantity we operate on.

[jokasimr]

Suggested change

# tofs = distance / simulation.speed

[jokasimr]

There's one thing I've been thinking about related to this change, and it is that there might be a bit higher interpolation error when we interpolate wavelength instead of time-of-flight.

$$ \lambda = \frac{t - t_0(t, L)}{L} $$

while

$$ t_{of} = t - t_0(t, L), $$

assuming $t_0(t, L)$ is an affine function, or slowly varying relative to grid density, linear interpolation in $t$ and $L$ will represent $t_{tof}$ exactly.
But linear interpolation does not represent $\lambda$ exactly even if $t_0$ is slowly varying or affine, so the interpolation error is probably larger when we interpolate $\lambda$ than when we interpolate $t_{of}$.

How to fix this?

If we want to fix this we don't need to change the interfaces or reintroduce tof in the packages.
We just need to pass $t_0$ (or $t_{of}$) instead of wavelength to the interpolation implementation inside WavelengthInterpolator, and then internally in WavelengthInterpolator.__call__ compute wavelength from the interpolated t0.

[jokasimr]

It's very pleasant to review in the new monorepo :) And I like the diff!

[nvaytet]

There's one thing I've been thinking about related to this change, and it is that there might be a bit higher interpolation error when we interpolate wavelength instead of time-of-flight.

λ = t − t 0 ( t , L ) L

while

t o f = t − t 0 ( t , L ) ,

assuming t 0 ( t , L ) is an affine function, or slowly varying relative to grid density, linear interpolation in t and L will represent t t o f exactly. But linear interpolation does not represent λ exactly even if t 0 is slowly varying or affine, so the interpolation error is probably larger when we interpolate λ than when we interpolate t o f .

How to fix this?

If we want to fix this we don't need to change the interfaces or reintroduce tof in the packages. We just need to pass t 0 (or t o f ) instead of wavelength to the interpolation implementation inside WavelengthInterpolator, and then internally in WavelengthInterpolator.__call__ compute wavelength from the interpolated t0.

To get the wavelength, we are using the raw wavelengths given by the tof simulation and averaging and interpolating using them directly.
I'm not sure I understand why they are different, and why you have the extra L factor in the first equation? But maybe talking in-person is easier?

[jokasimr]

To get the wavelength, we are using the raw wavelengths given by the tof simulation and averaging and interpolating using them directly.

Yeah I know. Maybe talking in person is easier.

The point is, assuming $t_0(t, L)$ is piecewise affine (or almost constant/slowly varying):

• $t_{of}(t, L)$ is also piecewise affine. This means we can interpolate it with good accuracy using linear interpolation.
• $\lambda(t, L)$ is not an affine function of $t$ and $L$. This means there will be an interpolation error that is independent of how well $t_0$ can be represented on the interpolation grid.

Here is an illustration explaining the problem:

[nvaytet]

So I tested it., and you are indeed right.

Here is a comparison of the errors we would get inside a square on a very coarse lookup table for 100 random points.
In blue is the relative error on the computed wavelength, while in orange is the relative error on the computed tof.

As you predicted, the error when working with wavelengths is consistently larger than when operating on tofs.

So I will do what you suggest: work internally on tofs, and convert to wavelengths for outputs.

Thanks very much for catching this 🙏 it may have gone unnoticed for years... 😬

[nvaytet]

I attach the notebook I used to make the figure above
interpolation-errors-lut.ipynb

It is a bit of a mess. but at least it's there for the record.

[nvaytet]

How to fix this?

If we want to fix this we don't need to change the interfaces or reintroduce tof in the packages.
We just need to pass t0 (or tof) instead of wavelength to the interpolation implementation inside WavelengthInterpolator, and then internally in WavelengthInterpolator.__call__ compute wavelength from the interpolated t0.

It's not quite that simple I think?
We should have lookup tables that store the tof (instead of wavelength) as a function of (toa, distance), right?
So bascially, we want to do the same as we were doing before these changes, but then convert to wavelength using transform coords internally before we return the data, I think?

Edit: I think I would still rather have wavelength lookup tables stored on disk, and then convert them to tofs internally, because when users will want to mask regions of the table with large uncertainties, it will probably be easier for them to reason in wavelength than in tof (what is the wavelength resolution I am after?)

[jokasimr]

I'd prefer that we keep the WavelengthInterpolator but when calling the InterpolatorImpl we pass lookup.data.values * self._distance_edges and in WavelengthInterpolator.__call__ we return
values=self._interpolator(...) / ltotal (but making that operation in-place is probably better).

Then I don't think we have to make other changes elsewhere or re-introduce tof, right?

[nvaytet]

Yes, I thought about it some more, and basically did just what you suggest here.

[nvaytet]

@YooSunYoung can you take a look to make sure I didn't mess up anything in essnmx?