code: def time_correlation_by_diagonalization(P, pi, obs1, obs2=None, time=1, rdl=None):
      
      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
docs: 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.