code: def _nnls(A, b, maxiter=None, tol=None):
      """
      This is a single RHS algorithm from ref [2] above. For multiple RHS
      support, the algorithm is given in  :doi:`10.1002/cem.889`
      """
      m, n = A.shape
  
      AtA = A.T @ A
      Atb = b @ A  # Result is 1D - let NumPy figure it out
  
      if not maxiter:
          maxiter = 3 * n
      if tol is None:
          tol = 10 * max(m, n) * np.spacing(1.0)
  
      # Initialize vars
      x = np.zeros(n, dtype=np.float64)
  
      # Inactive constraint switches
      P = np.zeros(n, dtype=bool)
  
      # Projected residual
      resid = Atb.copy().astype(np.float64)  # x=0. Skip (-AtA @ x) term
  
      # Overall iteration counter
      # Outer loop is not counted, inner iter is counted across outer spins
      iter = 0
  
      while (not P.all()) and (resid[~P] > tol).any():  # B
          # Get the "most" active coeff index and move to inactive set
          resid[P] = -np.inf
          k = np.argmax(resid)  # B.2
          P[k] = True  # B.3
  
          # Iteration solution
          s = np.zeros(n, dtype=np.float64)
          P_ind = P.nonzero()[0]
          s[P] = solve(
              AtA[P_ind[:, None], P_ind[None, :]],
              Atb[P],
              assume_a="sym",
              check_finite=False,
          )  # B.4
  
          # Inner loop
          while (iter < maxiter) and (s[P].min() <= tol):  # C.1
              alpha_ind = ((s < tol) & P).nonzero()
              alpha = (
                  x[alpha_ind] / (x[alpha_ind] - s[alpha_ind])
              ).min()  # C.2
              x *= 1 - alpha
              x += alpha * s
              P[x < tol] = False
              s[P] = solve(
                  AtA[np.ix_(P, P)],
                  Atb[P],
                  assume_a="sym",
                  check_finite=False,
              )
              s[~P] = 0  # C.6
              iter += 1
  
          x[:] = s[:]
          resid = Atb - AtA @ x
  
          if iter == maxiter:
              # Typically following line should return
              # return x, np.linalg.norm(A@x - b), -1
              # however at the top level, -1 raises an exception wasting norm
              # Instead return dummy number 0.
              return x, 0.0, -1
  
      return x, np.linalg.norm(A @ x - b), 1
  
label_desc: Use image node_id 2 for calling a global function with example usage: _nnls(A, b, maxiter, tol) and returns: x,  np,  int, x,  int,  unknown