Singular matrix in `max_vol` for EIM w/ low-energy Coulomb scattering

[beykyle]

To reproduce (in branch feature-kd-global):

import rose
import numpy as np
from matplotlib import pyplot as plt

# set up kinematics
from rose.koning_delaroche import KDGlobal, Projectile

# for 27-Al
A = 27
Z = 13

# lab bombarding energy
Elab = 0.1 # MeV

# get kinematics and default KD params
omp = rose.koning_delaroche.KDGlobal(Projectile.proton)
(mu, Ecom, k, eta, R_C), parameters = omp.get_params(A, Z, Elab)`

interactions = rose.InteractionEIMSpace(
    rose.koning_delaroche.KD_simple,
    len(parameters),
    mu,
    Ecom,
    is_complex=True,
    spin_orbit_potential=rose.koning_delaroche.KD_simple_so,
    training_info=np.array(
        [parameters - np.fabs(0.5 * parameters), parameters + np.fabs(0.5 * parameters)]
    ).T,
    Z_1=1,
    Z_2=13,
    R_C=R_C,
    l_max=1,
)

Error:

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[3], line 1
----> 1 interactions = rose.InteractionEIMSpace(
      2     rose.koning_delaroche.KD_simple,
      3     len(parameters),
      4     mu,
      5     Ecom,
      6     is_complex=True,
      7     spin_orbit_potential=rose.koning_delaroche.KD_simple_so,
      8     training_info=np.array(
      9         [parameters - np.fabs(0.5 * parameters), parameters + np.fabs(0.5 * parameters)]
     10     ).T,
     11     Z_1=1,
     12     Z_2=13,
     13     R_C=R_C,
     14     l_max=1,
     15 )

File ~/umich/rose/src/rose/interaction_eim.py:261, in InteractionEIMSpace.__init__(self, coordinate_space_potential, n_theta, mu, energy, training_info, l_max, Z_1, Z_2, R_C, is_complex, spin_orbit_potential, n_basis, explicit_training, n_train, rho_mesh, match_points)
    258 else:
    259     for l in range(l_max+1):
    260         self.interactions.append(
--> 261             [InteractionEIM(coordinate_space_potential, n_theta, mu,
    262                 energy, l, training_info, Z_1=Z_1, Z_2=Z_2, R_C=R_C,
    263                 is_complex=is_complex,
    264                 spin_orbit_term=SpinOrbitTerm(spin_orbit_potential,
    265                 lds), n_basis=n_basis,
    266                 explicit_training=explicit_training, n_train=n_train,
    267                 rho_mesh=rho_mesh, match_points=match_points) for lds in
    268                 couplings(l)]
    269         )

File ~/umich/rose/src/rose/interaction_eim.py:261, in <listcomp>(.0)
    258 else:
    259     for l in range(l_max+1):
    260         self.interactions.append(
--> 261             [InteractionEIM(coordinate_space_potential, n_theta, mu,
    262                 energy, l, training_info, Z_1=Z_1, Z_2=Z_2, R_C=R_C,
    263                 is_complex=is_complex,
    264                 spin_orbit_term=SpinOrbitTerm(spin_orbit_potential,
    265                 lds), n_basis=n_basis,
    266                 explicit_training=explicit_training, n_train=n_train,
    267                 rho_mesh=rho_mesh, match_points=match_points) for lds in
    268                 couplings(l)]
    269         )

File ~/umich/rose/src/rose/interaction_eim.py:141, in InteractionEIM.__init__(self, coordinate_space_potential, n_theta, mu, energy, ell, training_info, Z_1, Z_2, R_C, is_complex, spin_orbit_term, n_basis, explicit_training, n_train, rho_mesh, match_points)
    139 di = i_max // (n_basis - 1)
    140 i_init = np.arange(0, i_max + 1, di)
--> 141 self.match_indices = max_vol(self.snapshots, i_init)
    142     # np.random.randint(0, self.snapshots.shape[0], size=self.snapshots.shape[1]))
    143 self.match_points = rho_mesh[self.match_indices]

File ~/umich/rose/src/rose/interaction_eim.py:36, in max_vol(basis, indxGuess)
     33     indexing[[ ij,indxGuess[ij] ]] = indexing[[ indxGuess[ij],ij ]]
     35 for iIn in range(1, 100):
---> 36     B = np.dot(interpBasis, np.linalg.inv(interpBasis[:nbases]))
     37     b = np.max(B)
     38     if b > 1:

File <__array_function__ internals>:180, in inv(*args, **kwargs)

File ~/mambaforge/envs/om/lib/python3.10/site-packages/numpy/linalg/linalg.py:552, in inv(a)
    550 signature = 'D->D' if isComplexType(t) else 'd->d'
    551 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 552 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    553 return wrap(ainv.astype(result_t, copy=False))

File ~/mambaforge/envs/om/lib/python3.10/site-packages/numpy/linalg/linalg.py:89, in _raise_linalgerror_singular(err, flag)
     88 def _raise_linalgerror_singular(err, flag):
---> 89     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix