calculates time correlation. Raises P to power 'times' by diagonalization.
    If rdl tuple (R, D, L) is given, it will be used for
    further calculation.

def time_correlation_by_diagonalization(P, pi, obs1, obs2=None, time=1, rdl=None):
    """
    calculates time correlation. Raises P to power 'times' by diagonalization.
    If rdl tuple (R, D, L) is given, it will be used for
    further calculation.
    """
    if rdl is None:
        raise ValueError("no rdl decomposition")
    R, D, L = rdl

    d_times = np.diag(D) ** time
    diag_inds = np.diag_indices_from(D)
    D_time = np.zeros(D.shape, dtype=d_times.dtype)
    D_time[diag_inds] = d_times
    P_time = np.dot(np.dot(R, D_time), L)

    # multiply element-wise obs1 and pi. this is obs1' diag(pi)
    l = np.multiply(obs1, pi)
    m = np.dot(P_time, obs2)
    result = np.dot(l, m)
    return result