def scale_joint_sfs(s):
    """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.

    """

    i = np.arange(s.shape[0])[:, None]
    j = np.arange(s.shape[1])[None, :]
    out = (s * i) * j
    return out