code: def recursive_response(
      gtab,
      data,
      mask=None,
      sh_order=8,
      peak_thr=0.01,
      init_fa=0.08,
      init_trace=0.0021,
      iter=8,
      convergence=0.001,
      parallel=False,
      num_processes=None,
      sphere=default_sphere,
  ):
      """Recursive calibration of response function using peak threshold
  
      Parameters
      ----------
      gtab : GradientTable
      data : ndarray
          diffusion data
      mask : ndarray, optional
          mask for recursive calibration, for example a white matter mask. It has
          shape `data.shape[0:3]` and dtype=bool. Default: use the entire data
          array.
      sh_order : int, optional
          maximal spherical harmonics order. Default: 8
      peak_thr : float, optional
          peak threshold, how large the second peak can be relative to the first
          peak in order to call it a single fiber population [1]. Default: 0.01
      init_fa : float, optional
          FA of the initial 'fat' response function (tensor). Default: 0.08
      init_trace : float, optional
          trace of the initial 'fat' response function (tensor). Default: 0.0021
      iter : int, optional
          maximum number of iterations for calibration. Default: 8.
      convergence : float, optional
          convergence criterion, maximum relative change of SH
          coefficients. Default: 0.001.
      parallel : bool, optional
          Whether to use parallelization in peak-finding during the calibration
          procedure. Default: True
      num_processes : int, optional
          If `parallel` is True, the number of subprocesses to use
          (default multiprocessing.cpu_count()). If < 0 the maximal number of
          cores minus ``num_processes + 1`` is used (enter -1 to use as many
          cores as possible). 0 raises an error.
      sphere : Sphere, optional.
          The sphere used for peak finding. Default: default_sphere.
  
      Returns
      -------
      response : ndarray
          response function in SH coefficients
  
      Notes
      -----
      In CSD there is an important pre-processing step: the estimation of the
      fiber response function. Using an FA threshold is not a very robust method.
      It is dependent on the dataset (non-informed used subjectivity), and still
      depends on the diffusion tensor (FA and first eigenvector),
      which has low accuracy at high b-value. This function recursively
      calibrates the response function, for more information see [1].
  
      References
      ----------
      .. [1] Tax, C.M.W., et al. NeuroImage 2014. Recursive calibration of
             the fiber response function for spherical deconvolution of
             diffusion MRI data.
      """
      S0 = 1.0
      evals = fa_trace_to_lambdas(init_fa, init_trace)
      res_obj = (evals, S0)
  
      if mask is None:
          data = data.reshape(-1, data.shape[-1])
      else:
          data = data[mask]
  
      n = np.arange(0, sh_order + 1, 2)
      where_dwi = lazy_index(~gtab.b0s_mask)
      response_p = np.ones(len(n))
  
      for _ in range(iter):
          r_sh_all = np.zeros(len(n))
          csd_model = ConstrainedSphericalDeconvModel(
              gtab, res_obj, sh_order=sh_order
          )
  
          csd_peaks = peaks_from_model(
              model=csd_model,
              data=data,
              sphere=sphere,
              relative_peak_threshold=peak_thr,
              min_separation_angle=25,
              parallel=parallel,
              num_processes=num_processes,
          )
  
          dirs = csd_peaks.peak_dirs
          vals = csd_peaks.peak_values
          single_peak_mask = (vals[:, 1] / vals[:, 0]) < peak_thr
          data = data[single_peak_mask]
          dirs = dirs[single_peak_mask]
  
          for num_vox in range(data.shape[0]):
              rotmat = vec2vec_rotmat(
                  dirs[num_vox, 0], np.array([0, 0, 1])
              )
  
              rot_gradients = np.dot(rotmat, gtab.gradients.T).T
  
              x, y, z = rot_gradients[where_dwi].T
              r, theta, phi = cart2sphere(x, y, z)
              # for the gradient sphere
              B_dwi = real_sh_descoteaux_from_index(
                  0, n, theta[:, None], phi[:, None]
              )
              r_sh_all += np.linalg.lstsq(
                  B_dwi, data[num_vox, where_dwi], rcond=-1
              )[0]
  
          response = r_sh_all / data.shape[0]
          res_obj = AxSymShResponse(
              data[:, gtab.b0s_mask].mean(), response
          )
  
          change = abs((response_p - response) / response_p)
          if all(change < convergence):
              break
  
          response_p = response
  
      return res_obj
  
imports: import warnings
  import numpy
  import numbers
  from scipy.integrate import quad
  from scipy.special import lpn, gamma
  import scipy.linalg
  import scipy.linalg.lapack
  from dipy.data import small_sphere, get_sphere, default_sphere
  from dipy.core.geometry import cart2sphere
  from dipy.core.ndindex import ndindex
  from dipy.sims.voxel import single_tensor
  from dipy.reconst.multi_voxel import multi_voxel_fit
  from dipy.reconst.dti import TensorModel, fractional_anisotropy
  from dipy.reconst.shm import sph_harm_ind_list, real_sh_descoteaux_from_index, sph_harm_lookup, lazy_index, SphHarmFit, real_sh_descoteaux, sh_to_rh, forward_sdeconv_mat, SphHarmModel
  from dipy.reconst.utils import _roi_in_volume, _mask_from_roi
  from dipy.direction.peaks import peaks_from_model
  from dipy.core.geometry import vec2vec_rotmat
  from dipy.utils.deprecator import deprecate_with_version, deprecated_params
label_desc: Use image node_id 13 for calling a global function with example usage: recursive_response(gtab, data, mask, sh_order, peak_thr, init_fa, init_trace, iter, convergence, parallel, num_processes, sphere) and returns: res_obj