code: def peak_reproducibility(
      data,
      vardata,
      domain,
      ngroups,
      sigma,
      method="crfx",
      swap=False,
      verbose=0,
      **kwargs,
  ):
      """Return a measure of cluster-level reproducibility
      of activation patterns
      (i.e. how far clusters are from each other)
  
      Parameters
      ----------
      data: array of shape (nvox,nsubj)
            the input data from which everything is computed
      vardata: array of shape (nvox,nsubj)
               the variance of the data that is also available
      domain: referential- and domain-defining image
      ngroups (int),
               Number of subbgroups to be drawn
      sigma: float, parameter that encodes how far far is
      threshold: float, binarization threshold
      method: string to be chosen among 'crfx', 'cmfx' or 'cffx',
             inference method under study
      swap = False: if True, a random sign swap of the data is performed
           This is used to simulate a null hypothesis on the data.
      verbose=0 : verbosity mode
  
      Returns
      -------
      score (float): the desired  cluster-level reproducibility index
      """
      tiny = 1.0e-15
      nsubj = data.shape[1]
      samples = draw_samples(nsubj, ngroups)
      all_pos = []
  
      # compute the positions in the different subgroups
      for i in range(ngroups):
          x = data[:, samples[i]]
  
          if swap:
              # apply a random sign swap to x
              x *= 2 * (np.random.rand(len(samples[i])) > 0.5) - 1
  
          if method != "crfx":
              vx = vardata[:, samples[i]]
          if method != "bsa":
              threshold = kwargs["threshold"]
  
              if method == "crfx":
                  stat_map = ttest(x)
              elif method == "cmfx":
                  stat_map = mfx_ttest(x, vx)
              elif method == "cffx":
                  stat_map = fttest(x, vx)
              elif method == "cjt":
                  if "k" in kwargs:
                      k = kwargs["k"]
                  else:
                      k = nsubj / 2
                  stat_map = conjunction(x, vx, k)
  
              pos = get_peak_position_from_thresholded_map(
                  stat_map, domain, threshold
              )
              all_pos.append(pos)
          else:
              # method='bsa' is a special case
              tx = x / (tiny + np.sqrt(vx))
              afname = kwargs["afname"]
              theta = kwargs["theta"]
              dmax = kwargs["dmax"]
              ths = kwargs["ths"]
              thq = kwargs["thq"]
              smin = kwargs["smin"]
              niter = kwargs["niter"]
              afname = afname + "_%02d_%04d.pic" % (niter, i)
              pos = coord_bsa(
                  domain, tx, theta, dmax, ths, thq, smin, afname
              )
              all_pos.append(pos)
  
      # derive a kernel-based goodness measure from the pairwise comparison
      # of sets of positions
      score = 0
      for i in range(ngroups):
          for j in range(i):
              score += statistics_from_position(
                  all_pos[i], all_pos[j], sigma
              )
              score += statistics_from_position(
                  all_pos[j], all_pos[i], sigma
              )
      score /= ngroups * (ngroups - 1)
      return score
  
imports: import numpy
  from nipy.io.nibcompat import get_affine
  from nipy.labs.spatial_models.discrete_domain import grid_domain_from_binary_array
label_desc: Use image node_id 16 for calling a global function with example usage: peak_reproducibility(data, vardata, domain, ngroups, sigma, method, swap, verbose) and returns: score