code: def joint_sfs(dac1, dac2, n1=None, n2=None):
      
  
      # check inputs
      dac1, n1 = _check_dac_n(dac1, n1)
      dac2, n2 = _check_dac_n(dac2, n2)
  
      # compute site frequency spectrum
      x = n1 + 1
      y = n2 + 1
      # need platform integer for bincount
      tmp = (dac1 * y + dac2).astype(int, copy=False)
      s = np.bincount(tmp)
      s.resize(x, y)
      return s
docs: Compute the joint site frequency spectrum between two populations.
  
      Parameters
      ----------
      dac1 : array_like, int, shape (n_variants,)
          Derived allele counts for the first population.
      dac2 : array_like, int, shape (n_variants,)
          Derived allele counts for the second population.
      n1, n2 : int, optional
          The total number of chromosomes called in each population.
  
      Returns
      -------
      joint_sfs : ndarray, int, shape (m_chromosomes, n_chromosomes)
          Array where the (i, j)th element is the number of variant sites with i
          derived alleles in the first population and j derived alleles in the
          second population.