code: def pairwise_dxy(pos, gac, start=None, stop=None, is_accessible=None):
      
  
      if not isinstance(pos, SortedIndex):
          pos = SortedIndex(pos, copy=False)
      gac = asarray_ndim(gac, 3)
      # compute this once here, to avoid repeated evaluation within the loop
      gan = np.sum(gac, axis=2)
      m = gac.shape[1]
      dist = list()
      for i, j in itertools.combinations(range(m), 2):
          ac1 = gac[:, i, ...]
          an1 = gan[:, i]
          ac2 = gac[:, j, ...]
          an2 = gan[:, j]
          d = sequence_divergence(pos, ac1, ac2, an1=an1, an2=an2,
                                  start=start, stop=stop,
                                  is_accessible=is_accessible)
          dist.append(d)
      return np.array(dist)
docs: Convenience function to calculate a pairwise distance matrix using
      nucleotide divergence (a.k.a. Dxy) as the distance metric.
  
      Parameters
      ----------
      pos : array_like, int, shape (n_variants,)
          Variant positions.
      gac : array_like, int, shape (n_variants, n_samples, n_alleles)
          Per-genotype allele counts.
      start : int, optional
          Start position of region to use.
      stop : int, optional
          Stop position of region to use.
      is_accessible : array_like, bool, shape (len(contig),), optional
          Boolean array indicating accessibility status for all positions in the
          chromosome/contig.
  
      Returns
      -------
      dist : ndarray
          Distance matrix in condensed form.
  
      See Also
      --------
      allel.model.ndarray.GenotypeArray.to_allele_counts