Replace maths library with Python libraries

[timothy-nunn]

Removing some of the maths routines that can be replaced with simple routines from libraries like numpy and scipy.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 80.88235% with 13 lines in your changes missing coverage. Please review.

Project coverage is 31.37%. Comparing base (67f90e6) to head (972dd21).

Files with missing lines	Patch %	Lines
process/pfcoil.py	14.28%	12 Missing ⚠️
process/availability.py	50.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3418      +/-   ##
==========================================
+ Coverage   31.21%   31.37%   +0.15%     
==========================================
  Files          81       81              
  Lines       19432    19477      +45     
==========================================
+ Hits         6065     6110      +45     
  Misses      13367    13367              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[timothy-nunn]

The changes to PF coil require a bit of explaining. The single value decomposition from scipy returns different results to the maths library implementation. This causes the integration tests to fail, but none of the regression tests.

In some cases, like test_solv the values of ccls changed, but their sum did not.
In test_efc, the ccls values changed (one of the PF coils groups ended up with 0 current!?) but the ssq value did not, weirdly, ssq is calculated from the ccls.

Overall, I suspect this is nothing to worry about, but its worth careful consideration before approving.

[jonmaddock]

I see, this is tricky. There are some quite substantial (to me) differences in the test results; we just need to be happy that this is due to differences in the SVD implementations and not any mistakes in the refactoring. It's particularly strange about the sum and ssq values, as you say.

[jonmaddock]

Great, just the question of the SVD differences. I'm not sure how serious these differences are; it looks correct, but maybe we need to convince ourselves fully.

[jonmaddock]

I'm sure this is intentional, but can you explain why curpfb and curfps values have been swapped?

[timothy-nunn]

Im not sure this is happening anymore

[jonmaddock]

These are quite different, as you say.

[jonmaddock]

I see, this is tricky. There are some quite substantial (to me) differences in the test results; we just need to be happy that this is due to differences in the SVD implementations and not any mistakes in the refactoring. It's particularly strange about the sum and ssq values, as you say.

[timothy-nunn]

The SVD differences can be explained as follows. The original SVD implementation does not return a correctly factorised matrix. The following code demonstrates that when attempting to recreate the original matrix (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.svd.html) for the test_solv, the outputs do not create a 3x3 matrix of all 2's when running with OLD=True, but do when running with OLD=False (the new integration test values).

import numpy as np

OLD = True

if OLD:
    U = np.array([
        [-0.81649658, -0.57735027, 0.0],
        [0.40824829, -0.57735027, -0.70710678],
        [0.40824829, -0.57735027, 0.70710678],
    ])
    V = np.array([
        [-0.81649658, -0.57735027, 0.0],
        [0.40824829, -0.57735027, -0.70710678],
        [0.40824829, -0.57735027, 0.70710678],
    ])
    s = np.array([5.1279005e-16, 6.0000000e00, 0.0000000e00])
else:
    U = np.array([
        [-0.57735027, -0.57735027, -0.57735027],
        [-0.57735027, -0.21132487, 0.78867513],
        [-0.57735027, 0.78867513, -0.21132487],
    ])
    V = np.array([
        [-0.57735027, -0.57735027, -0.57735027],
        [0.0, -0.70710678, 0.70710678],
        [0.81649658, -0.40824829, -0.40824829],
    ])
    s = np.array([6.0, 0.0, 0.0])


# reconstruct
sigma = np.zeros((3, 3))
for i in range(3):
    sigma[i, i] = s[i]

a1 = np.dot(U, np.dot(sigma, V))

print(f"Checking {'old' if OLD else 'new'} SVD results")
print("Original matrix was 3x3 all 2's")
print(a1)