Clean up voltageToEfieldConverter and use exact solution if possible

[sjoerd-bouma]

Removes some unused code from the voltageToEfieldConverter, and uses the exact solution for the matrix inversion in the case of only two antennas. This was previously discussed in #527. It is useful e.g. for LOFAR where unfolding is always done per antenna pair. The output is unchanged, but it makes things a bit quicker (the inversion is sped up by ~a factor 2)

MWE below to demonstrate that the output is unchanged:

import numpy as np
from NuRadioReco.utilities import fft

def stacked_lstsq(L, b, rcond=1e-10): # old version of the function
    """
    Solve L x = b, via SVD least squares cutting of small singular values
    L is an array of shape (..., M, N) and b of shape (..., M).
    Returns x of shape (..., N)
    """
    u, s, v = np.linalg.svd(L, full_matrices=False)
    s_max = s.max(axis=-1, keepdims=True)
    s_min = rcond * s_max
    inv_s = np.zeros_like(s)
    inv_s[s >= s_min] = 1 / s[s >= s_min]
    x = np.einsum('...ji,...j->...i', v,
                  inv_s * np.einsum('...ji,...j->...i', u, b.conj()))
    return np.conj(x, x)

def solve_2x2(L, b):
    denom = (L[:,0,0] * L[:,1,1] - L[:,0,1] * L[:,1,0])
    e_theta = (b[:,0] * L[:,1,1] - b[:,1] * L[:,0,1]) / denom
    e_phi = (b[:,1] - L[:,1,0] * e_theta) / L[:,1,1]
    return np.stack((e_theta, e_phi), axis=-1)


np.random.seed(1234)
L = np.random.normal(size=(1025, 2, 2)) + 1.j * np.random.normal(size=(1025, 2, 2))
b = fft.time2freq(np.random.normal(size=(2, 2048)), sampling_rate=2).T

sol_svd = stacked_lstsq(L, b)
sol_ana = np.sum(np.linalg.inv(L)* b[:, None], axis=-1)
sol_ana_2x2 = solve_2x2(L, b)

print(np.allclose(sol_svd, sol_ana, rtol=1e-8, atol=1e-8))
print(np.allclose(sol_svd, sol_ana_2x2, rtol=1e-8, atol=1e-8))
>>> True
>>> True

[cg-laser]

nice cleanup, thanks!

I always liked to have the explicit equations in the code, but your generic analytic matrix inversion code is neater as it generalized to N dimensions. However, it might be slightly slower than the explicit implementation of the inverse for two dimensions? Then we could also catch the special case with two dimensions and use the explicit equation? Up to you.

[sjoerd-bouma]

nice cleanup, thanks!

I always liked to have the explicit equations in the code, but your generic analytic matrix inversion code is neater as it generalized to N dimensions. However, it might be slightly slower than the explicit implementation of the inverse for two dimensions? Then we could also catch the special case with two dimensions and use the explicit equation? Up to you.

Thanks, good point. I've added it back in (and updated the MWE above). The performance does improve by a lot. Unfortunately the bottleneck is not in this part of the voltageToEfieldConverter so we don't gain much, but I suppose it doesn't hurt either.

[cg-laser]

great, thanks