docs: def get_hypergeometric_stats(N, indices):
      """Calculates hypergeom. p-values and fold enrichments for all cutoffs.
  
      Parameters
      ----------
      N: int
          The length of the list
      indices:  `numpy.ndarray` with ``dtype=np.uint16``
          The (sorted) indices of the "1's" in the list.
      """
      assert isinstance(N, (int, np.integer))
      assert isinstance(indices, np.ndarray) and \
              np.issubdtype(indices.dtype, np.uint16)
  
      K = indices.size
  
      pvals = np.empty(N+1, dtype=np.float64)
      folds = np.empty(N+1, dtype=np.float64)
      pvals[0] = 1.0
      folds[0] = 1.0
  
      n = 0
      k = 0
      p = 1.0
      while n < N:
          if k < K and indices[k] == n:
              # "add one"
              # calculate f(k+1; N,K,n+1) from f(k; N,K,n)
              p *= (float((n+1) * (K-k)) / \
                    float((N-n) * (k+1)))
              k += 1
          else:
              # "add zero"
              # calculate f(k; N,K,n+1) from f(k; N,K,n)
              p *= (float((n+1) * (N-K-n+k)) /
                    float((N-n) * (n-k+1)))
          n += 1
          # calculate hypergeometric p-value
          pvals[n] = get_hgp(p, k, N, K, n)
          # calculate fold enrichment
          folds[n] = k / (K*(n/float(N)))
  
      return pvals, folds