Enable time-of-flight indexing and Laue/ToF refinement

[toastisme]

Apologies for this being a large PR - it's difficult to separate out the different components. This adds indexing and refinement for time-of-flight data and Laue data more broadly. The main contributions are Laue specific refinement classes that also minimise degrees of freedom with respect to the calculated and observed wavelength of each reflection.

Usage:
dials.index imported.expt strong.refl reflections.weighting_strategy.override=<constant, statistical> reflections.weighting_strategy.wavelength_weight=<weight>

Background

The diffraction condition can be written

$$\mathbf{p^*_0}\cdot \mathbf{p^*_0} + 2\mathbf{p^*_0}\cdot\mathbf{S_0} = 0,$$

where $\mathbf{p^*_0}$ is the reciprocal lattice vector obtained from the candidate UB matrix for a given Miller index. For a rotation experiment we satisfy the condition by identifying the angle $\mathbf{p^*_0}$ must be rotated by to cross the Ewald sphere. For a Laue experiment, we can instead find the required wavelength by rearranging the diffraction condition to give

$$\lambda = -2\frac{\mathbf{\hat{S_0}} \cdot \mathbf{p^*_0}}{\mathbf{p^*_0}\cdot \mathbf{p^*_0}}.$$

The target function to minimise during refinement is a weighted sum of squared residuals over all $n$ observed reflections

$$L = \frac{1}{2}\sum^n_{i=1}w_{i,X}(X-X_{obs})^2 + w_{i,Y}(Y-Y_{obs})^2 + w_{i,\lambda}(\lambda - \lambda_{obs})^2,$$

where $X, Y, \lambda$ correspond to calculated reflection centroid positions in the x, y panel positions, and a wavelength calculated as above, respectively. Likewise, $X_{obs}, Y_{obs}, \lambda_{obs}$ refer to the centroid positions obtained from spot finding. $w_{i,X}, w_{i,Y}, w_{i,\lambda}$ are weight coefficients based on the strategy selected by the user. $L$ has first derivatives

$$\frac{\partial L}{\partial p} = \sum^n_{i=1}w_{i,X}\frac{\partial X}{\partial p}(X-X_{obs}) + w_{i,Y}\frac{\partial Y}{\partial p}(Y-Y_{obs}) + w_{i,\lambda}\frac{\partial \lambda}{\partial p}(\lambda - \lambda_{obs}),$$

where $p$ refers to a degree of freedom of the model obtained from indexing (e.g. unit cell angles, lengths, panel positions etc.). To refine centroids with respect to their calculated wavelengths we need to calculate $\frac{\partial \lambda}{\partial p}$. To do this we follow the same approach as for rotational data[1] and use the implicit function theorem. This states that if a given function $G(\lambda, p) = c$, $\lambda(p)$ has derivatives

$$\frac{\partial}{\partial p}\lambda(p) = -\frac{\frac{\partial G}{\partial p}}{\frac{\partial G}{\partial \lambda}}.$$

In this case $G$ is the diffraction condition, and hence

$$\frac{\partial \lambda}{\partial p} = -\frac{\frac{\partial}{\partial p}(\mathbf{p^*_0}\cdot \mathbf{p^*_0} + 2\mathbf{p^*_0}\cdot\mathbf{S_0})}{\frac{\partial}{\partial\lambda}(\mathbf{p^*_0}\cdot \mathbf{p^*_0} + 2\mathbf{p^*_0}\cdot\mathbf{S_0})}.$$

Terms in the numerator have already been worked out in DIALS[1], so only the denominator needs calculating.

For the first term in the denominator,

$$\frac{\partial}{\partial \lambda}\mathbf{p^*_0}\cdot\mathbf{p^*_0} = \mathbf{p^*_0}\cdot\frac{\partial \mathbf{p^*_0}}{\partial \lambda} + \frac{\partial \mathbf{p^*_0}}{\partial \lambda}\cdot\mathbf{p^*_0}$$ $$\begin{split} \frac{\partial \mathbf{p^*_0}}{\partial \lambda} &= \frac{\partial}{\partial\lambda}(\mathbf{S'}-\mathbf{S_0})\\\ &= \frac{\partial}{\partial\lambda}(\mathbf{S'}-\mathbf{\hat{S_0}}\lambda^{-1})\\\ &= \frac{1}{\lambda}\mathbf{S_0}, \end{split}$$

and hence

$$\frac{\partial}{\partial \lambda}\mathbf{p^*_0}\cdot\mathbf{p^*_0} = \frac{2}{\lambda}\mathbf{p^*_0}\cdot\mathbf{S_0}.$$

For the second term,

$$\frac{\partial}{\partial\lambda}(\mathbf{p^*_0}\cdot\mathbf{S_0}) = \mathbf{p^*_0}\cdot\frac{\partial \mathbf{S_0}}{\partial \lambda} + \mathbf{S_0}\cdot\frac{\partial \mathbf{p^*_0}}{\partial \lambda}$$ $$\frac{\partial\mathbf{S_0}}{\partial\lambda} = -\frac{\mathbf{S_0}}{\lambda},$$

and so

$$\begin{split} \frac{\partial}{\partial\lambda}(\mathbf{p^*_0}\cdot\mathbf{S_0}) &= \frac{1}{\lambda}(\mathbf{S_0}\cdot\mathbf{S_0}-\mathbf{p^*_0}\cdot\mathbf{S_0}). \end{split}$$

Putting it all together

$$\frac{\partial}{\partial\lambda}(\mathbf{p^*_0}\cdot \mathbf{p^*_0} + 2\mathbf{p^*_0}\cdot\mathbf{S_0}) = \frac{2}{\lambda}\mathbf{S_0}\cdot\mathbf{S_0}$$ $$\begin{split} \frac{\partial \lambda}{\partial p} &= -\frac{\frac{\partial}{\partial p}(\mathbf{p^*_0}\cdot \mathbf{p^*_0} + 2\mathbf{p^*_0}\cdot\mathbf{S_0})}{\frac{2}{\lambda}\mathbf{S_0}\cdot\mathbf{S_0}}\\\ &= -\frac{\mathbf{p^*_0}\cdot\frac{\partial\mathbf{p^*_0}}{\partial p} + \mathbf{S_0}\cdot\frac{\partial\mathbf{p^*_0}}{\partial p} + \mathbf{p^*_0}\frac{\partial\mathbf{S_0}}{\partial p}}{\frac{1}{\lambda}\mathbf{S_0}\cdot\mathbf{S_0}}\\\ &= -\frac{\lambda(\frac{\partial\mathbf{p^*_0}}{\partial p}\cdot\mathbf{S'}+\mathbf{p^*_0}\cdot\frac{\partial\mathbf{S_0}}{\partial p})}{\mathbf{S_0}\cdot\mathbf{S_0}}, \end{split}$$

[1] D. G. Waterman, G. Winter, R. J. Gildea, J. M. Parkhurst, A. S. Brewster, N. K. Sauter, and G. Evans.
Diffraction-geometry refinement in the DIALS framework. Acta Crystallographica Section D, 72(4):558–
575, Apr 2016.

Weighting strategies

I looked at constant and statistical weighting strategies for $w_{\lambda}$. At least for time-of-flight experiments, the variance along the wavelength direction for each spot is less obviously useful, as the spot profile typically has a long tail. Instead I used the variance in the X and Y multiplied by a constant. Results for different values are given below. exp (blue dashed line) refers to expected results for a given unit cell parameter. Looking at the average mean squared error across the datasets, statistical weighting seems generally better, with the best weighting being 1E4.

[biochem-fan]

This adds indexing and refinement for time-of-flight data and Laue data more broadly. The main contributions are Laue specific refinement classes that also minimise degrees of freedom with respect to the calculated and observed wavelength of each reflection.

Just to make sure, is this PR for Laue data collected with ToF detectors? In other words, we cannot use this on X-ray Laue or pink beam datasets without ToF information (i.e. unknown per-spot $\lambda_{obs}$), can we?

[dagewa]

I'm not sure this is the right error message. This class gets called in quite a number of places, not just by the indexer where I assume the max cell issue is coming from. So we might end up here in some other situation. Maybe it should be "Cannot create an ExperimentsPredictor for ToF and non-ToF experiments at the same time"

[toastisme]

This adds indexing and refinement for time-of-flight data and Laue data more broadly. The main contributions are Laue specific refinement classes that also minimise degrees of freedom with respect to the calculated and observed wavelength of each reflection.

Just to make sure, is this PR for Laue data collected with ToF detectors? In other words, we cannot use this on X-ray Laue or pink beam datasets without ToF information (i.e. unknown per-spot λobs), can we?

I've written it so you should be able to use this for any data where you have assigned wavelengths to observed reflections, and you have enough confidence in the assignment to minimise your model against them. All your data should need is a wavelength column in your reflection table. The only thing missing is a consistent way for DIALS to recognise you have a Laue experiment in the absence of a reflection table, which is the motivation for #733.

If you have any Laue data I could use a basis of a test that would be useful!

[biochem-fan]

@toastisme
Most of public datasets are from pink beam serial crystallography (e.g. https://zenodo.org/records/8354296, https://www.cxidb.org/id-180.html)

The only Laue dataset I found is https://zenodo.org/records/10199220.

[dagewa]

There are quite a lot of new classes added here, but no added tests.

• LauePredictionParameterisation
• LaueReflectionManager
• LaueReflectionPredictor
• LaueExperimentsPredictor
• TOFExperimentsPredictor
• LaueReflectionManager
• TOFReflectionManager
• TOFLeastSquaresResidualWithRmsdCutoff
• LaueLeastSquaresResidualWithRmsdCutoff
• LaueStatisticalWeightingStrategy
• LaueMixedWeightingStrategy

I don't think it is necessary to write individual tests for all of these, as not all of the rotation equivalents of these classes are tested directly in a unit test like manner either. However, some of them are. For example, derivatives are tested versus finite difference approximations in places like test_beam_parameters.py, test_stills_spherical_relp_derivatives.py, etc., and for the full prediction equation parameterisation in test_prediction_parameters.py.

A similar test vs finite-differences might be useful for LauePredictionParameterisation.

There's also a test of the gradients of the target function LeastSquaresPositionalResidualWithRmsdCutoff in test_finite_diffs.py (bad module name). Perhaps something similar for TOFLeastSquaresResidualWithRmsdCutoff / LaueLeastSquaresResidualWithRmsdCutoff?

The whole machinery of rotation method refinement is also tested against generated reflection positions using ideal geometry in test_orientation_refinement.py. Might something like that be useful here?

[toastisme]

There are quite a lot of new classes added here, but no added tests.

* `LauePredictionParameterisation`

* `LaueReflectionManager`

* `LaueReflectionPredictor`

* `LaueExperimentsPredictor`

* `TOFExperimentsPredictor`

* `LaueReflectionManager`

* `TOFReflectionManager`

* `TOFLeastSquaresResidualWithRmsdCutoff`

* `LaueLeastSquaresResidualWithRmsdCutoff`

* `LaueStatisticalWeightingStrategy`

* `LaueMixedWeightingStrategy`

I don't think it is necessary to write individual tests for all of these, as not all of the rotation equivalents of these classes are tested directly in a unit test like manner either. However, some of them are. For example, derivatives are tested versus finite difference approximations in places like test_beam_parameters.py, test_stills_spherical_relp_derivatives.py, etc., and for the full prediction equation parameterisation in test_prediction_parameters.py.

A similar test vs finite-differences might be useful for LauePredictionParameterisation.

There's also a test of the gradients of the target function LeastSquaresPositionalResidualWithRmsdCutoff in test_finite_diffs.py (bad module name). Perhaps something similar for TOFLeastSquaresResidualWithRmsdCutoff / LaueLeastSquaresResidualWithRmsdCutoff?

The whole machinery of rotation method refinement is also tested against generated reflection positions using ideal geometry in test_orientation_refinement.py. Might something like that be useful here?

Yes completely agree. Thanks for highlighting the key areas to test. I had some dials data / dials data files prs to add some neutron data to write these kind of tests that are now merged. I was hoping to have some X-ray Laue data too.