Document how to use a component reading to create a new source

[nvaytet]

We add section to the docs where we show how to use a component reading to create a more efficient source.

We histogram the neutrons at the detector in wavelength and birth_time, smooth the histogram, and use it as a new distribution to sample from.

Can speed up calculations.
We should consider making the final function into a classmethod of the Source: Source.from_reading(...).

[jokasimr]

Interesting! Could give a big speedup in some cases.

However, I'd like to suggest a modification. Currently the probability of sampling neutrons from point $(\lambda, t, t^o)$ where ($t$ is the time in the detector, and $t^o$ is the birth time) with the new source will change from what it was using the original source. It might be close to the same, and maybe it does not matter in practice. But we can actually make them be the same, with a little bit more work.

If we have a source distribution $q_0$ and we have identified a region $(\lambda, t^o) \in M$ such that outside $M$ the probability the neutron passes the chopper cascade is essentially zero. Then we can define a new source probability density

$$q(\lambda, t^o) = \mathcal{1}_M \frac{q_0(\lambda, t^o)}{\int_M q_0(\lambda, t^o) d\lambda \ d t^o}$$

such that the probability density of $\lambda, t^o$ given that the neutron passed the chopper cascade is the same for the source $q_0$ as it is for $q$.

Intuitively, with the new source we only increase the probability of sampling in the region where we actually need to sample, but in that region, the ratio of the probabilities of sampling two points does not change.

After applying the gaussian kernel to $p$ I'd say we have identified $M$. It's the region where $p>0$ (or perhaps where $p>\varepsilon$ for some small $\varepsilon$ because the gaussian kernel might introduce some nonzero weight everywhere, I'm not sure).

I'm imagining it would look something like this:

def source_from_reading(\n",
    reading,
    original_source,
    neutrons: int,
    time_smooth_kernel: sc.Variable | None = None,
    wavelength_smooth_kernel: sc.Variable | None = None,
):
    from scipp.scipy.ndimage import gaussian_filter

    p = reading.data.hist(wavelength=original_source.coords['wavelength'], birth_time=original_source.coords['birth_time'])
    if time_smooth_kernel is None:
        time_smooth_kernel = 2 * time_binning
    if wavelength_smooth_kernel is None:
        wavelength_smooth_kernel = 2 * wavelength_binning
    p = gaussian_filter(
        p,
        sigma={
            'birth_time': time_smooth_kernel,
            'wavelength': wavelength_smooth_kernel,
        },
    )
    q = original_source * (p > 1e-8) / (original_source * (p > 1e-8)).sum()
    return tof.Source.from_distribution(neutrons=neutrons, p=q)"

[nvaytet]

Thanks for looking into this 👍

So you are saying that you'd rather apply a mask on the original source. I like the idea.
However, as it's written now

q = original_source * (p > 1e-8) / (original_source * (p > 1e-8)).sum()

would require p to have the same shape as original_source.

Basically it would be really nice if we could come up with a procedure where the mask we put on the source does not depend on either the initial binning for the histogram, nor the width of the gaussian kernel.
My intuition is that we can probably achieve the former but not the latter (i.e. the gaussian smoothing will always have an effect).

Currently, the ess source is sampled from a data array with a defined shape, so we could use that original shape to define our bins for histogramming?

[jokasimr]

So you are saying that you'd rather apply a mask on the original source.

Yes that's it. Because then the distribution in the detector will be unchanged, assuming we sampled a reasonable number of points to figure out what region of the source makes it to the detector.

As written now, would require p to have the same shape as original_source.
It would be really nice if we could come up with a procedure where the mask we put on the source does not depend on either the initial binning for the histogram, nor the width of the gaussian kernel.

I don't understand why we want the shape of the new source to be independent of the original source.
Why does that matter?

If we want the new source to be finer than the original source, then we can first up-sample the original source (for example by using scipy.ndimage.zoom, or scipp.rebin) to the desired shape, and then do the same thing.

Making the new source coarser than the original source does not seem useful, but I might be missing something.

[nvaytet]

I meant that to do an operation such as q = original_source * p, Scipp requires the DataArrays p and original_source to have the same shape.

In my current version, p is defined by some binning I randomly chose when looking at the detector events.

I had missed that in the code you proposed, you suggested to use for histogramming the bin edges from the original source, which ensures consistency.