docs: def bootstrap(nside, rand, nbar, *data):
      """ This function will bootstrap data based on the sky coverage of rand.
          It is different from bootstrap in the traditional sense, but for correlation
          functions it gives the correct answer with less computation.
  
          nbar : number density of rand, used to estimate the effective area of a pixel
  
          nside : number of healpix pixels per side to use
  
          *data : a list of data -- will be binned on the same regions.
  
          small regions (incomplete pixels) are combined such that the total
          area is about the same (a healpix pixel) in each returned boot strap sample
  
          Yields: area, random, *data
  
          rand and *data are in (RA, DEC)
  
          Example:
  
          >>> for area, ran, data1, data2 in bootstrap(4, ran, 100., data1, data2):
          >>>    # Do stuff
          >>>    pass
      """
  
      def split(data, indices, axis):
          """ This function splits array. It fixes the bug
              in numpy that zero length array are improperly handled.
  
              In the future this will be fixed.
          """
          s = []
          s.append(slice(0, indices[0]))
          for i in range(len(indices) - 1):
              s.append(slice(indices[i], indices[i+1]))
          s.append(slice(indices[-1], None))
  
          rt = []
          for ss in s:
              ind = [slice(None, None, None) for i in range(len(data.shape))]
              ind[axis] = ss
              ind = tuple(ind)
              rt.append(data[ind])
          return rt
  
      def hpsplit(nside, data):
          # data is (RA, DEC)
          RA, DEC = data
          pix = radec2pix(nside, RA, DEC)
          n = numpy.bincount(pix)
          a = numpy.argsort(pix)
          data = numpy.array(data)[:, a]
          rt = split(data, n.cumsum(), axis=-1)
          return rt
  
      # mean area of sky.
  
      Abar =  41252.96 / nside2npix(nside)
      rand = hpsplit(nside, rand)
      if len(data) > 0:
          data = [list(i) for i in zip(*[hpsplit(nside, d1) for d1 in data])]
      else:
          data = [[] for i in range(len(rand))]
  
      heap = []
      j = 0
      for r, d in zip(rand, data):
          if len(r[0]) == 0: continue
          a = 1.0 * len(r[0]) / nbar
          j = j + 1
          if len(heap) == 0:
              heapq.heappush(heap, (a, j, r, d))
          else:
              a0, j0, r0, d0 = heapq.heappop(heap)
              if a0 + a < Abar:
                  a0 += a
                  d0 = [
                       numpy.concatenate((d0[i], d[i]), axis=-1)
                       for i in range(len(d))
                      ]
                  r0 = numpy.concatenate((r0, r), axis=-1)
              else:
                  heapq.heappush(heap, (a, j, r, d))
              heapq.heappush(heap, (a0, j0, r0, d0))
  
      for i in range(len(heap)):
          area, j, r, d = heapq.heappop(heap)
          rt = [area, r] + d
          yield rt