code: def scale_joint_sfs(s):
      
  
      i = np.arange(s.shape[0])[:, None]
      j = np.arange(s.shape[1])[None, :]
      out = (s * i) * j
      return out
docs: Scale a joint site frequency spectrum.
  
      Parameters
      ----------
      s : array_like, int, shape (n1, n2)
          Joint site frequency spectrum.
  
      Returns
      -------
      joint_sfs_scaled : ndarray, int, shape (n1, n2)
          Scaled joint site frequency spectrum.