r"""Time-correlation for equilibrium experiment.

    Parameters
    ----------
    T : (M, M) ndarray or scipy.sparse matrix
        Transition matrix
    obs1 : (M,) ndarray
        Observable, represented as vector on state space
    obs2 : (M,) ndarray (optional)
        Second observable, for cross-correlations
    times : array-like of int (optional), default=(1)
        List of times (in tau) at which to compute correlation
    maxtime : int, optional, default=None
        Maximum time step to use. Equivalent to . Alternative to times.
    k : int (optional)
        Number of eigenvalues and eigenvectors to use for computation
    ncv : int (optional)
        The number of Lanczos vectors generated, `ncv` must be greater than k;
        it is recommended that ncv > 2*k

    Returns
    -------
    correlations : ndarray
        Correlation values at given times
    times : ndarray, optional
        time points at which the correlation was computed (if return_times=True)

    References
    ----------
    .. [1] Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D
        Chodera and J Smith. 2010. Dynamical fingerprints for probing
        individual relaxation processes in biomolecular dynamics with
        simulations and kinetic experiments. PNAS 108 (12): 4822-4827.

    Notes
    -----

    **Auto-correlation**

    The auto-correlation of an observable :math:`a(x)` for a system in
    equilibrium is

    .. math:: \mathbb{E}_{\mu}[a(x,0)a(x,t)]=\sum_x \mu(x) a(x, 0) a(x, t)

    :math:`a(x,0)=a(x)` is the observable at time :math:`t=0`.  It can
    be propagated forward in time using the t-step transition matrix
    :math:`p^{t}(x, y)`.

    The propagated observable at time :math:`t` is :math:`a(x,
    t)=\sum_y p^t(x, y)a(y, 0)`.

    Using the eigenvlaues and eigenvectors of the transition matrix
    the autocorrelation can be written as

    .. math:: \mathbb{E}_{\mu}[a(x,0)a(x,t)]=\sum_i \lambda_i^t \langle a, r_i\rangle_{\mu} \langle l_i, a \rangle.

    **Cross-correlation**

    The cross-correlation of two observables :math:`a(x)`,
    :math:`b(x)` is similarly given

    .. math:: \mathbb{E}_{\mu}[a(x,0)b(x,t)]=\sum_x \mu(x) a(x, 0) b(x, t)

    Examples
    --------

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

    >>> 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, 0.0])
    >>> times = np.array([1, 5, 10, 20])

    >>> corr = correlation(T, a, times=times)
    >>> corr
    array([ 0.40909091,  0.34081364,  0.28585667,  0.23424263])

def correlation(T, obs1, obs2=None, times=(1), maxtime=None, k=None, ncv=None, return_times=False):
    r"""Time-correlation for equilibrium experiment.

    Parameters
    ----------
    T : (M, M) ndarray or scipy.sparse matrix
        Transition matrix
    obs1 : (M,) ndarray
        Observable, represented as vector on state space
    obs2 : (M,) ndarray (optional)
        Second observable, for cross-correlations
    times : array-like of int (optional), default=(1)
        List of times (in tau) at which to compute correlation
    maxtime : int, optional, default=None
        Maximum time step to use. Equivalent to . Alternative to times.
    k : int (optional)
        Number of eigenvalues and eigenvectors to use for computation
    ncv : int (optional)
        The number of Lanczos vectors generated, `ncv` must be greater than k;
        it is recommended that ncv > 2*k

    Returns
    -------
    correlations : ndarray
        Correlation values at given times
    times : ndarray, optional
        time points at which the correlation was computed (if return_times=True)

    References
    ----------
    .. [1] Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D
        Chodera and J Smith. 2010. Dynamical fingerprints for probing
        individual relaxation processes in biomolecular dynamics with
        simulations and kinetic experiments. PNAS 108 (12): 4822-4827.

    Notes
    -----

    **Auto-correlation**

    The auto-correlation of an observable :math:`a(x)` for a system in
    equilibrium is

    .. math:: \mathbb{E}_{\mu}[a(x,0)a(x,t)]=\sum_x \mu(x) a(x, 0) a(x, t)

    :math:`a(x,0)=a(x)` is the observable at time :math:`t=0`.  It can
    be propagated forward in time using the t-step transition matrix
    :math:`p^{t}(x, y)`.

    The propagated observable at time :math:`t` is :math:`a(x,
    t)=\sum_y p^t(x, y)a(y, 0)`.

    Using the eigenvlaues and eigenvectors of the transition matrix
    the autocorrelation can be written as

    .. math:: \mathbb{E}_{\mu}[a(x,0)a(x,t)]=\sum_i \lambda_i^t \langle a, r_i\rangle_{\mu} \langle l_i, a \rangle.

    **Cross-correlation**

    The cross-correlation of two observables :math:`a(x)`,
    :math:`b(x)` is similarly given

    .. math:: \mathbb{E}_{\mu}[a(x,0)b(x,t)]=\sum_x \mu(x) a(x, 0) b(x, t)

    Examples
    --------

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

    >>> 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, 0.0])
    >>> times = np.array([1, 5, 10, 20])

    >>> corr = correlation(T, a, times=times)
    >>> corr
    array([ 0.40909091,  0.34081364,  0.28585667,  0.23424263])

    """
    # check if square matrix and remember size
    T = _types.ensure_ndarray_or_sparse(T, ndim=2, uniform=True, kind='numeric')
    n = T.shape[0]
    obs1 = _types.ensure_ndarray(obs1, ndim=1, size=n, kind='numeric')
    obs2 = _types.ensure_ndarray_or_None(obs2, ndim=1, size=n, kind='numeric')
    times = _types.ensure_int_vector(times, require_order=True)

    # check input
    # go
    if _issparse(T):
        return sparse.fingerprints.correlation(T, obs1, obs2=obs2, times=times, k=k, ncv=ncv)
    else:
        return dense.fingerprints.correlation(T, obs1, obs2=obs2, times=times, k=k)