ADR: Convolution

[henrikjacobsenfys]

Introduction to Convolution

This documents the design decisions for the Convolution class.

Convolution consists of four classes and a data structure: Convolution, AnalyticalConvolution, NumericalConvolution, ConvolutionBase and energy_grid. Users are only expected to interact indirectly with Convolution via the Analysis class. It takes the following input:
energy: The energy where the convolution is calculated.
sample_components: ComponentCollection or ModelComponent
resolution_components: ComponentCollection or ModelComponent
energy_offset: a fitable Parameter or number that shifts the energy axis slightly, to account for e.g. instrument misalignment.
various settings for numerical convolution.

Convolution builds a convolution_plan, which separates the sample_components into three groups: delta functions, analytical components (Gauss, Lorentzian, Voigt) and the rest.

The primary method, convolution calculates the convolution of these three groups. The delta functions are handled in the Convolution class, as they are simple. The analytical components are passed to an AnalyticalConvolution object, and the numerical components are passed to a NumericalConvolution object along with any extra args for control the accuracy and speed. If temperature is included then all but the delta functions are handled numerically.

Convolution inherits from NumericalConvolutionBase

ConvolutionBase

ConvolutionBase contains basic functionality shared among all convolutions: Input validation, properties, setters and getters. To make it easy to handle energy_offset it has a energy_with_offset property/

AnalyticalConvolution

AnalyticalConvolution inherits from ConvolutionBase. It contains lists of all possible analytical pairs and methods to calculate their convolution. For reference: all areas are multiplied when calculating convolutions.
The convolution of a Gaussian with a Gaussian is another Gaussian with width (standard deviation) $\sigma = \sqrt{\sigma_1^2 + \sigma_2^2}
The convolution of a Lorentzian with a Lorentzian is another Lorentzian with width (hald width at half max) $\Gamma = \Gamma_1 + \Gamma_2
The convolution of a Gaussian with a Lorentzian is a Voigt, which can be calculated efficiently using voigt_profile from scipy.special.
The convolution of a Gaussian with a Voigt adds the Gaussian widths in quadrature and leaves the Lorentzian width unchanged.
The convolution of a Lorentzian with a Voigt adds the Lorentzian widths and leaves the Gaussian width unchanged.
The convolution of a Voigt with a Voigt adds the Gaussian widths in quadrature and adds the Lorentzian widths.

EnergyGrid

To improve numerical accuracy, the numerical convolution is calculated on an upsampled energy grid as described below. All the relevant values are saved in an EnergyGrid, which contains:
energy_dense: the dense energy_grid
energy_dense_centered: centered energy grid which is used to evaluate the resolution. (for example, if data goes from -5 to +15 meV, the convolution would give incorrect results if the resolution was evaluated here.)
energy_dense_step: the step size of the dense grid for normalization
energy_span_dense: the total span of the dense energy grid for warnings
energy_even_length_offset: a small offset that needs to be applied when the length of the energy is even

NumericalConvolutionBase

Inherits from ConvolutionBase. Also takes the following input:
upsample_factor: defaults to 5, which means the convolution is evaluated on 5 times as many points as the input for accuracy
extension_factor: defaults to 0.2, which means the energy grid is extended by 10% of the span in both directions. I.e. if the energy goes from -1 to +1 meV, the calculations will be from -1.1 to +1.1 meV. This is to avoid edge effects
'temperature': for detailed balancing
'normalize_detailed_balance': whether or not to normalize detailed balancing
It also takes units for temperature and energy, but these should rarely be required.

It builds the extended energy grid and updates it when e.g. energy changes. The grid is not updated when energy_offset is changed. The reason is that energy_offset should always be small and is a fit parameter. Recreating the grid at every step in a fit would be inefficient.

When the length of the energy is even, the fft-based convolution gives a slightly wrong result. The convolution of two arrays of length N is of length 2N-1. When using 'same' mode, only the central N points are kept, so the output has the same length as the input. However, if N is even, the center falls between two points, leading to a half-bin offset.
For example, if N=4, the convolution has length 7, and when we select the 4 central points we either get indices [2,3,4,5] or [1,2,3,4], both of which are offset by 0.5*dx from the true center at index 3.5. The effect is small but non-zero.

'NumericalConvolutionBase' also has methods to check if the width of any functions is too narrow or too large for numerical stability. In such cases, warnings are given.The thresholds are not guaranteed to be correct in all cases, but should serve as a guide for the user. They are showcased in tests/performance_tests/convolution

NumericalConvolution

Inherits from NumericalConvolutionBase and only adds the convolution method. At the end it interpolates back to the original energy values from the grid.

Possible extensions in the future

Another class to handle convolutions with measured resolutions is planned. This class may require methods to subtract background, smooth the data and perhaps upsample the data using interpolation.