<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def rdl_decomposition_rev(T, norm='reversible', mu=None):
 r"""Decomposition into left and right eigenvectors for reversible transition matrices. Parameters T : (M, M) ndarray Transition matrix norm: {'standard', 'reversible'} standard: (L'R) = Id, L[:,0] is a probability distribution, the stationary distribution mu of T. Right eigenvectors R have a 2-norm of 1. reversible: R and L are related via L=L[:,0]*R. mu : (M,) ndarray, optional Stationary distribution of T Returns ------- R : (M, M) ndarray The normalized (with respect to L) right eigenvectors, such that the column R[:,i] is the right eigenvector corresponding to the eigenvalue w[i], dot(T,R[:,i])=w[i]*R[:,i] D : (M, M) ndarray A diagonal matrix containing the eigenvalues, each repeated according to its multiplicity L : (M, M) ndarray The normalized (with respect to `R`) left eigenvectors, such that the row ``L[i, :]`` is the left eigenvector corresponding to the eigenvalue ``w[i]``, ``dot(L[i, :], T)``=``w[i]*L[i, :]`` Notes ----- The eigenvalues and eigenvectors of the similar symmetric matrix `\sqrt(\mu_i / \mu_j) p_{ij}` will be used to compute the eigenvalues and eigenvectors of T. The stationay distribution will be computed if no precomputed stationary distribution is given. """

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

    """Diagonal matrix of eigenvalues"""
    D = np.diag(val)

    """Right and left eigenvectors"""
    R = eigvec / smu[:, np.newaxis]
    L = eigvec * smu[:, np.newaxis]

    """Ensure that R[:,0] is positive and unity"""
    tmp = R[0, 0]
    R[:, 0] = R[:, 0] / tmp

    """Ensure that L[:, 0] is probability vector"""
    L[:, 0] = L[:, 0] *  tmp

    if norm == 'reversible':
        return R, D, L.T
    elif norm == 'standard':
        """Standard l2-norm of right eigenvectors"""
        w = np.diag(np.dot(R.T, R))
        sw = np.sqrt(w)
        """Don't change normalization of eigenvectors for dominant eigenvalue"""
        sw[0] = 1.0

        R = R / sw[np.newaxis, :]
        L = L * sw[np.newaxis, :]
        return R, D, L.T
    else:
        raise ValueError("Keyword 'norm' has to be either 'standard' or 'reversible'")