Rescale matrices and vectors in invert_regularised_nnls and allow kwargs.

CHANGELOG.md

@@ -15,6 +15,7 @@ New:
 * Improved parsing of metadata from the ADAS ADF15 'bnd' files for H-like ions. Raises a runtime error if the metadata cannot be parsed. (#424)
 * **Beam dispersion calculation has changed from sigma(z) = sigma + z * tan(alpha) to sigma(z) = sqrt(sigma^2 + (z * tan(alpha))^2) for consistancy with the Gaussian beam model. Attention!!! The results of BES and CX spectroscopy are affected by this change. (#414)**
 * Improved beam direction calculation to allow for natural broadening of the BES line shape due to beam divergence. (#414)
+* Add kwargs to invert_regularised_nnls to pass them to scipy.optimize.nnls. (#438)
 
 Bug fixes:
 * Fix deprecated transforms being cached in LaserMaterial after laser.transform update (#420)

cherab/tools/inversions/nnls.py

@@ -21,21 +21,25 @@
 import scipy
 
 
-def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None):
-    """
+def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None, **kwargs):
+    r"""
     Solves :math:`\mathbf{b} = \mathbf{W} \mathbf{x}` for the vector :math:`\mathbf{x}`,
     using Tikhonov regulariastion.
 
     This is a thin wrapper around scipy.optimize.nnls, which modifies
     the arguments to include the supplied Tikhonov regularisation matrix.
 
+    The values of w_matrix, b_vector and alpha * tikhonov_matrix are notmalised
+    by max(b_vector) before passing them to scipy.optimize.nnls().
+
     :param np.ndarray w_matrix: The sensitivity matrix describing the coupling between the
       detectors and the voxels. Must be an array with shape :math:`(N_d, N_s)`.
     :param np.ndarray b_vector: The measured power/radiance vector with shape :math:`(N_d)`.
     :param float alpha: The regularisation hyperparameter :math:`\alpha` which determines
       the regularisation strength of the tikhonov matrix.
     :param np.ndarray tikhonov_matrix: The tikhonov regularisation matrix operator, an array
       with shape :math:`(N_s, N_s)`. If None, the identity matrix is used.
+    :param **kwargs: Keyword arguments passed to scipy.optimize.nnls.
     :return: (x, norm), the solution vector and the residual norm.
 
     .. code-block:: pycon
@@ -60,6 +64,9 @@ def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None
     d_vector = np.zeros(m+n)
     d_vector[0:m] = b_vector[:]
 
-    x_vector, rnorm = scipy.optimize.nnls(c_matrix, d_vector)
+    # Normalise c_matrix and d_vector to avoid possible issues with the nnls termination criteria.
+    vmax = d_vector.max()
+
+    x_vector, rnorm = scipy.optimize.nnls(c_matrix / vmax, d_vector / vmax, **kwargs)
 
-    return x_vector, rnorm
+    return x_vector, rnorm * vmax