r"""Compute eigenvectors of given transition matrix.

    Parameters
    ----------
    T : numpy.ndarray, shape(d,d) or scipy.sparse matrix
        Transition matrix (stochastic matrix)
    k : int (optional)
        Compute the first k eigenvectors
    ncv : int (optional)
        The number of Lanczos vectors generated, `ncv` must be greater than `k`;
        it is recommended that `ncv > 2*k`
    right : bool, optional
        If right=True compute right eigenvectors, left eigenvectors
        otherwise
    reversible : bool, optional
        Indicate that transition matrix is reversible
    mu : (M,) ndarray, optional
        Stationary distribution of T


    Returns
    -------
    eigvec : numpy.ndarray, shape=(d, n)
        The eigenvectors of T ordered with decreasing absolute value of
        the corresponding eigenvalue. If k is None then n=d, if k is
        int then n=k.

    See also
    --------
    rdl_decomposition

    Notes
    -----
    Eigenvectors are computed using the scipy interface
    to the corresponding LAPACK/ARPACK routines.

    If reversible=False, the returned eigenvectors :math:`v_i` are
    normalized such that

    ..  math::

        \langle v_i, v_i \rangle = 1

    This is the case for right eigenvectors :math:`r_i` as well as
    for left eigenvectors :math:`l_i`.

    If you desire orthonormal left and right eigenvectors please use the
    rdl_decomposition method.

    If reversible=True the the eigenvectors of the similar symmetric
    matrix `\sqrt(\mu_i / \mu_j) p_{ij}` will be used to compute the
    eigenvectors of T.

    The precomputed stationary distribution will only be used if
    reversible=True.

    Examples
    --------

    >>> import numpy as np
    >>> from msmtools.analysis import eigenvectors

    >>> T = np.array([[0.9, 0.1, 0.0], [0.5, 0.0, 0.5], [0.0, 0.1, 0.9]])
    >>> R = eigenvectors(T)

    Matrix with right eigenvectors as columns

    >>> R # doctest: +ELLIPSIS
    array([[  5.77350269e-01,   7.07106781e-01,   9.90147543e-02], ...

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

    Parameters
    ----------
    T : numpy.ndarray, shape(d,d) or scipy.sparse matrix
        Transition matrix (stochastic matrix)
    k : int (optional)
        Compute the first k eigenvectors
    ncv : int (optional)
        The number of Lanczos vectors generated, `ncv` must be greater than `k`;
        it is recommended that `ncv > 2*k`
    right : bool, optional
        If right=True compute right eigenvectors, left eigenvectors
        otherwise
    reversible : bool, optional
        Indicate that transition matrix is reversible
    mu : (M,) ndarray, optional
        Stationary distribution of T


    Returns
    -------
    eigvec : numpy.ndarray, shape=(d, n)
        The eigenvectors of T ordered with decreasing absolute value of
        the corresponding eigenvalue. If k is None then n=d, if k is
        int then n=k.

    See also
    --------
    rdl_decomposition

    Notes
    -----
    Eigenvectors are computed using the scipy interface
    to the corresponding LAPACK/ARPACK routines.

    If reversible=False, the returned eigenvectors :math:`v_i` are
    normalized such that

    ..  math::

        \langle v_i, v_i \rangle = 1

    This is the case for right eigenvectors :math:`r_i` as well as
    for left eigenvectors :math:`l_i`.

    If you desire orthonormal left and right eigenvectors please use the
    rdl_decomposition method.

    If reversible=True the the eigenvectors of the similar symmetric
    matrix `\sqrt(\mu_i / \mu_j) p_{ij}` will be used to compute the
    eigenvectors of T.

    The precomputed stationary distribution will only be used if
    reversible=True.

    Examples
    --------

    >>> import numpy as np
    >>> from msmtools.analysis import eigenvectors

    >>> T = np.array([[0.9, 0.1, 0.0], [0.5, 0.0, 0.5], [0.0, 0.1, 0.9]])
    >>> R = eigenvectors(T)

    Matrix with right eigenvectors as columns

    >>> R # doctest: +ELLIPSIS
    array([[  5.77350269e-01,   7.07106781e-01,   9.90147543e-02], ...

    """
    T = _types.ensure_ndarray_or_sparse(T, ndim=2, uniform=True, kind='numeric')
    if _issparse(T):
        if right:
            return sparse.decomposition.eigenvectors(T, k=k, right=right, ncv=ncv)
        else:
            return sparse.decomposition.eigenvectors(T, k=k, right=right, ncv=ncv).T

    else:
        if right:
            return dense.decomposition.eigenvectors(T, k=k, right=right)
        else:
            return dense.decomposition.eigenvectors(T, k=k, right=right).T