r"""Equilibrium expectation value of a given observable.

    Parameters
    ----------
    T : (M, M) ndarray or scipy.sparse matrix
        Transition matrix
    a : (M,) ndarray
        Observable vector
    mu : (M,) ndarray (optional)
        The stationary distribution of T.  If given, the stationary
        distribution will not be recalculated (saving lots of time)

    Returns
    -------
    val: float
        Equilibrium expectation value fo the given observable

    Notes
    -----
    The equilibrium expectation value of an observable a is defined as follows

    .. math::

        \mathbb{E}_{\mu}[a] = \sum_i \mu_i a_i

    :math:`\mu=(\mu_i)` is the stationary vector of the transition matrix :math:`T`.

    Examples
    --------

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

    >>> T = np.array([[0.9, 0.1, 0.0], [0.5, 0.0, 0.5], [0.0, 0.1, 0.9]])
    >>> a = np.array([1.0, 0.0, 1.0])
    >>> m_a = expectation(T, a)
    >>> m_a # doctest: +ELLIPSIS
    0.909090909...

def expectation(T, a, mu=None):
    r"""Equilibrium expectation value of a given observable.

    Parameters
    ----------
    T : (M, M) ndarray or scipy.sparse matrix
        Transition matrix
    a : (M,) ndarray
        Observable vector
    mu : (M,) ndarray (optional)
        The stationary distribution of T.  If given, the stationary
        distribution will not be recalculated (saving lots of time)

    Returns
    -------
    val: float
        Equilibrium expectation value fo the given observable

    Notes
    -----
    The equilibrium expectation value of an observable a is defined as follows

    .. math::

        \mathbb{E}_{\mu}[a] = \sum_i \mu_i a_i

    :math:`\mu=(\mu_i)` is the stationary vector of the transition matrix :math:`T`.

    Examples
    --------

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

    >>> T = np.array([[0.9, 0.1, 0.0], [0.5, 0.0, 0.5], [0.0, 0.1, 0.9]])
    >>> a = np.array([1.0, 0.0, 1.0])
    >>> m_a = expectation(T, a)
    >>> m_a # doctest: +ELLIPSIS
    0.909090909...

    """
    # check if square matrix and remember size
    T = _types.ensure_ndarray_or_sparse(T, ndim=2, uniform=True, kind='numeric')
    n = T.shape[0]
    a = _types.ensure_ndarray(a, ndim=1, size=n, kind='numeric')
    mu = _types.ensure_ndarray_or_None(mu, ndim=1, size=n, kind='numeric')
    # go
    if not mu:
        mu = stationary_distribution(T)
    return _np.dot(mu, a)