Efficient calculation in a (unpolarised) reflectivity analysis

Efficient calculation

The Abeles matrix transfer method is an efficient way of calculating reflectivity. In a calculation there are a few aspects that dominate calculation time:

1. sqrt of complex numbers to calculate wavevectors
2. exp of complex numbers to implement Nevot Croce smearing, and in calculation of characteristic matrices
3. matrix multiplication of the 2x2 characteristic matrices.

There are typically many datapoints measured (NQ), and potentially hundreds of layers (NL). Calculation over all these would typically require a nested loop in C (outer loop over NQ, then inner loop over NL), so thought should be given to how much calculation can be lifted out of those loops, and how to make the innermost loops efficient as possible. In Python (or Cython) those nested loops can be completely removed by using numpy vectorisation. However, in my experience this still isn't as fast as straight C. In Julia it's better to follow the C route of nested loops. When compiling C code there may be some compiler flags that are ok to use, that give minor speedups (e.g. -funsafe-math-optimizations, -ffinite-math-only). This will be compiler and platform specific, and make assumptions about input and output. It is not worth trying to write functions to perform complex arithmetic, those provided by the maths (C/C++) libraries are typically faster, and also lookout correct mathematical performance for things like overflow and the correct branch cuts to use when calculating square roots.

There are some language specific ways of performing quick matrix multiplication, it's typically faster to hardcode the matrix multiplication for each of the entries (but less clear as to what's being done).

As a relative guide for a two layer model with 1000 datapoints it should be possible to calculate a reflectivity curve in C in ~250 microseconds, or about 433 microseconds in Python (using numpy). The reflectivity calculation should be O(N) for both NQ and NL.

Lastly, the use of the fastest kernel calculation speed (C?) may not be desirable in all circumstances. For example many optimisation analysis packages (using Bayesian MCMC or scalar minimisation) may be able to use gradient or jacobian information in their operation. Historically these gradients have been obtained via numerical differentiation, but Automatic Differentiation (AD) packages such as JAX/Autograd (Python) or ForwardDiff.jl (Julia) can calculate gradients directly from the function of interest. These AD packages typically require that function are coded in a certain manner, and can't autodiff C code. The gains from AD may outweigh the cost of a slightly lower reflectivity calculation speed.