code: def rumba_deconv(
      data, kernel, n_iter=600, recon_type="smf", n_coils=1
  ):
      r"""
      Fit fODF and GM/CSF volume fractions for a voxel using RUMBA-SD [1]_.
  
      Deconvolves the kernel from the diffusion-weighted signal by computing a
      maximum likelihood estimation of the fODF. Minimizes the negative
      log-likelihood of the data under Rician or Noncentral Chi noise
      distributions by adapting the iterative technique developed in
      Richardson-Lucy deconvolution.
  
      Parameters
      ----------
      data : 1d ndarray (N,)
          Signal values for a single voxel.
      kernel : 2d ndarray (N, M)
          Deconvolution kernel mapping volume fractions of the M compartments to
          N-length signal. Last two columns should be for GM and CSF.
      n_iter : int, optional
          Number of iterations for fODF estimation. Must be a positive int.
          Default: 600
      recon_type : {'smf', 'sos'}, optional
          MRI reconstruction method: spatial matched filter (SMF) or
          sum-of-squares (SoS). SMF reconstruction generates Rician noise while
          SoS reconstruction generates Noncentral Chi noise. Default: 'smf'
      n_coils : int, optional
          Number of coils in MRI scanner -- only relevant in SoS reconstruction.
          Must be a positive int. Default: 1
  
      Returns
      -------
      fit_vec : 1d ndarray (M,)
          Vector containing fODF and GM/CSF volume fractions. First M-2
          components are fODF while last two are GM and CSF respectively.
  
      Notes
      -----
      The diffusion MRI signal measured at a given voxel is a sum of
      contributions from each intra-voxel compartment, including parallel white
      matter (WM) fiber populations in a given orientation as well as effects
      from GM and CSF. The equation governing these  contributions is:
  
      $S_i = S_0\left(\sum_{j=1}^{M}f_j\exp(-b_i\textbf{v}_i^T\textbf{D}_j
      \textbf{v}_i) + f_{GM}\exp(-b_iD_{GM})+f_{CSF}\exp(-b_iD_{CSF})\right)$
  
      Where $S_i$ is the resultant signal along the diffusion-sensitizing
      gradient unit vector $\textbf{v_i}; i = 1, ..., N$ with a b-value of $b_i$.
      $f_j; j = 1, ..., M$ is the volume fraction of the $j^{th}$ fiber
      population with an anisotropic diffusion tensor $\textbf{D_j}$.
  
      $f_{GM}$ and $f_{CSF}$ are the volume fractions and $D_{GM}$ and $D_{CSF}$
      are the mean diffusivities of GM and CSF respectively.
  
      This equation is linear in $f_j, f_{GM}, f_{CSF}$ and can be simplified to
      a single matrix multiplication:
  
      $\textbf{S} = \textbf{Hf}$
  
      Where $\textbf{S}$ is the signal vector at a certain voxel, $\textbf{H}$ is
      the deconvolution kernel, and $\textbf{f}$ is the vector of volume
      fractions for each compartment.
  
      Modern MRI scanners produce noise following a Rician or Noncentral Chi
      distribution, depending on their signal reconstruction technique [2]_.
      Using this linear model, it can be shown that the likelihood of a signal
      under a Noncentral Chi noise model is:
  
      $P(\textbf{S}|\textbf{H}, \textbf{f}, \sigma^2, n) = \prod_{i=1}^{N}\left(
      \frac{S_i}{\bar{S_i}}\right)^n\exp\left\{-\frac{1}{2\sigma^2}\left[
      S_i^2 + \bar{S}_i^2\right]\right\}I_{n-1}\left(\frac{S_i\bar{S}_i}
      {\sigma^2}\right)u(S_i)$
  
      Where $S_i$ and $\bar{S}_i = \textbf{Hf}$ are the measured and expected
      signals respectively, and $n$ is the number of coils in the scanner, and
      $I_{n-1}$ is the modified Bessel function of first kind of order $n-1$.
      This gives the likelihood under a Rician distribution when $n$ is set to 1.
  
      By taking the negative log of this with respect to $\textbf{f}$ and setting
      the derivative to 0, the $\textbf{f}$ maximizing likelihood is found to be:
  
      $\textbf{f} = \textbf{f} \circ \frac{\textbf{H}^T\left[\textbf{S}\circ
      \frac{I_n(\textbf{S}\circ \textbf{Hf}/\sigma^2)} {I_{n-1}(\textbf{S}
      \circ\textbf{Hf}\sigma^2)} \right ]} {\textbf{H}^T\textbf{Hf}}$
  
      The solution can be found using an iterative scheme, just as in the
      Richardson-Lucy algorithm:
  
      $\textbf{f}^{k+1} = \textbf{f}^k \circ \frac{\textbf{H}^T\left[\textbf{S}
      \circ\frac{I_n(\textbf{S}\circ\textbf{Hf}^k/\sigma^2)} {I_{n-1}(\textbf{S}
      \circ\textbf{Hf}^k/\sigma^2)} \right ]} {\textbf{H}^T\textbf{Hf}^k}$
  
      In order to apply this, a reasonable estimate of $\sigma^2$ is required.
      To find this, a separate iterative scheme is found using the derivative
      of the negative log with respect to $\sigma^2$, and is run in parallel.
      This is shown here:
  
      $\alpha^{k+1} = \frac{1}{nN}\left\{ \frac{\textbf{S}^T\textbf{S} +
      \textbf{f}^T\textbf{H}^T\textbf{Hf}}{2} - \textbf{1}^T_N\left[(\textbf{S}
      \circ\textbf{Hf})\circ\frac{I_n(\textbf{S}\circ\textbf{Hf}/\alpha^k)}
      {I_{n-1}(\textbf{S}\circ\textbf{Hf}/\alpha^k)} \right ]\right \}$
  
      For more details, see [1]_.
  
      References
      ----------
      .. [1] Canales-Rodríguez, E. J., Daducci, A., Sotiropoulos, S. N., Caruyer,
             E., Aja-Fernández, S., Radua, J., Mendizabal, J. M. Y.,
             Iturria-Medina, Y., Melie-García, L., Alemán-Gómez, Y., Thiran,
             J.-P.,Sarró, S., Pomarol-Clotet, E., & Salvador, R. (2015).
             Spherical Deconvolution of Multichannel Diffusion MRI Data with
             Non-Gaussian Noise Models and Spatial Regularization. PLOS ONE,
             10(10), e0138910. https://doi.org/10.1371/journal.pone.0138910
  
      .. [2] Constantinides, C. D., Atalar, E., & McVeigh, E. R. (1997).
             Signal-to-Noise Measurements in Magnitude Images from NMR Phased
             Arrays. Magnetic Resonance in Medicine: Official Journal of the
             Society of Magnetic Resonance in Medicine / Society of Magnetic
             Resonance in Medicine, 38(5), 852–857.
      """
  
      n_comp = kernel.shape[1]  # number of compartments
      n_grad = len(data)  # gradient directions
  
      fodf = np.ones((n_comp, 1))  # initial guess is iso-probable
      fodf = fodf / np.sum(fodf, axis=0)  # normalize initial guess
  
      if recon_type == "smf":
          n_order = (
              1  # Rician noise (same as Noncentral Chi with order 1)
          )
      elif recon_type == "sos":
          n_order = n_coils  # Noncentral Chi noise (order = # of coils)
      else:
          raise ValueError(
              "Invalid recon_type. Should be 'smf' or 'sos', "
              + f"received {recon_type}"
          )
  
      data = data.reshape(-1, 1)
      reblurred = np.matmul(kernel, fodf)
  
      # For use later
      kernel_t = kernel.T
      f_zero = 0  # minimum value allowed in fODF
  
      # Initialize variance
      sigma0 = 1 / 15
      sigma2 = sigma0**2 * np.ones(data.shape)  # Expand into vector
  
      reblurred_s = data * reblurred / sigma2
  
      for _ in range(n_iter):
          fodf_i = fodf
          ratio = mbessel_ratio(n_order, reblurred_s)
          rl_factor = np.matmul(kernel_t, data * ratio) / (
              np.matmul(kernel_t, reblurred) + _EPS
          )
  
          fodf = fodf_i * rl_factor  # result of iteration
          fodf = np.maximum(f_zero, fodf)  # apply positivity constraint
  
          # Update other variables
          reblurred = np.matmul(kernel, fodf)
          reblurred_s = data * reblurred / sigma2
  
          # Iterate variance
          sigma2_i = (1 / (n_grad * n_order)) * np.sum(
              (data**2 + reblurred**2) / 2
              - (sigma2 * reblurred_s) * ratio,
              axis=0,
          )
          sigma2_i = np.minimum(
              (1 / 8) ** 2, np.maximum(sigma2_i, (1 / 80) ** 2)
          )
          sigma2 = sigma2_i * np.ones(data.shape)  # Expand into vector
  
      # Normalize final result
      fit_vec = np.squeeze(fodf / (np.sum(fodf, axis=0) + _EPS))
  
      return fit_vec
  
imports: import logging
  import warnings
  import numpy
  from dipy.core.geometry import vec2vec_rotmat
  from dipy.core.gradients import gradient_table, unique_bvals_tolerance, get_bval_indices
  from dipy.core.onetime import auto_attr
  from dipy.core.sphere import Sphere
  from dipy.data import get_sphere
  from dipy.reconst.shm import lazy_index, normalize_data
  from dipy.reconst.odf import OdfModel, OdfFit
  from dipy.reconst.csdeconv import AxSymShResponse
  from dipy.segment.mask import bounding_box, crop
  from dipy.sims.voxel import single_tensor, all_tensor_evecs
label_desc: Use image node_id 1 for calling a global function with example usage: rumba_deconv(data, kernel, n_iter, recon_type, n_coils) and returns: fit_vec