code: def fold_joint_sfs(s, n1, n2):
      
  
      # check inputs
      s = asarray_ndim(s, 2)
      assert s.shape[0] <= n1 + 1, 'invalid number of chromosomes'
      assert s.shape[1] <= n2 + 1, 'invalid number of chromosomes'
  
      # need to check s has all entries up to m
      if s.shape[0] < n1 + 1:
          sm = np.zeros((n1 + 1, s.shape[1]), dtype=s.dtype)
          sm[:s.shape[0]] = s
          s = sm
  
      # need to check s has all entries up to n
      if s.shape[1] < n2 + 1:
          sn = np.zeros((s.shape[0], n2 + 1), dtype=s.dtype)
          sn[:, :s.shape[1]] = s
          s = sn
  
      # fold
      mf = (n1 + 1) // 2
      nf = (n2 + 1) // 2
      n1 = mf * 2
      n2 = nf * 2
      o = (s[:mf, :nf] +  # top left
           s[mf:n1, :nf][::-1] +  # top right
           s[:mf, nf:n2][:, ::-1] +  # bottom left
           s[mf:n1, nf:n2][::-1, ::-1])  # bottom right
  
      return o
docs: Fold a joint site frequency spectrum.
  
      Parameters
      ----------
      s : array_like, int, shape (m_chromosomes, n_chromosomes)
          Joint site frequency spectrum.
      n1, n2 : int, optional
          The total number of chromosomes called in each population.
  
      Returns
      -------
      joint_sfs_folded : ndarray, int
          Folded joint site frequency spectrum.