r"""Compute eigenvectors of reversible transition matrix.

    Parameters
    ----------
    T : (d, d) ndarray
        Transition matrix (stochastic matrix)
    right : bool, optional
        If right=True compute right eigenvectors, left eigenvectors
        otherwise
    mu : (d,) ndarray, optional
        Stationary distribution of T

    Returns
    -------
    eigvec : (d, d) ndarray
        The eigenvectors of T ordered with decreasing absolute value
        of the corresponding eigenvalue

def eigenvectors_rev(T, right=True, mu=None):
    r"""Compute eigenvectors of reversible transition matrix.

    Parameters
    ----------
    T : (d, d) ndarray
        Transition matrix (stochastic matrix)
    right : bool, optional
        If right=True compute right eigenvectors, left eigenvectors
        otherwise
    mu : (d,) ndarray, optional
        Stationary distribution of T

    Returns
    -------
    eigvec : (d, d) ndarray
        The eigenvectors of T ordered with decreasing absolute value
        of the corresponding eigenvalue

    """
    if mu is None:
        mu = stationary_distribution(T)
    """ symmetrize T """
    smu = np.sqrt(mu)
    S = smu[:,None] * T / smu
    val, eigvec = eigh(S)
    """Sort eigenvectors"""
    perm = np.argsort(np.abs(val))[::-1]
    eigvec = eigvec[:, perm]
    if right:
        return eigvec / smu[:, np.newaxis]
    else:
        return eigvec * smu[:, np.newaxis]