docs: def time_correlations_direct(P, pi, obs1, obs2=None, times=[1]):
      r"""Compute time-correlations of obs1, or time-cross-correlation with obs2.
  
      The time-correlation at time=k is computed by the matrix-vector expression:
      cor(k) = obs1' diag(pi) P^k obs2
  
  
      Parameters
      ----------
      P : ndarray, shape=(n, n) or scipy.sparse matrix
          Transition matrix
      obs1 : ndarray, shape=(n)
          Vector representing observable 1 on discrete states
      obs2 : ndarray, shape=(n)
          Vector representing observable 2 on discrete states. If not given,
          the autocorrelation of obs1 will be computed
      pi : ndarray, shape=(n)
          stationary distribution vector. Will be computed if not given
      times : array-like, shape(n_t)
          Vector of time points at which the (auto)correlation will be evaluated
  
      Returns
      -------
  
      """
      n_t = len(times)
      times = np.sort(times)  # sort it to use caching of previously computed correlations
      f = np.zeros(n_t)
  
      # maximum time > number of rows?
      if times[-1] > P.shape[0]:
          use_diagonalization = True
          R, D, L = rdl_decomposition(P)
          # discard imaginary part, if all elements i=0
          if not np.any(np.iscomplex(R)):
              R = np.real(R)
          if not np.any(np.iscomplex(D)):
              D = np.real(D)
          if not np.any(np.iscomplex(L)):
              L = np.real(L)
          rdl = (R, D, L)
  
      if use_diagonalization:
          for i in range(n_t):
              f[i] = time_correlation_by_diagonalization(P, pi, obs1, obs2, times[i], rdl)
      else:
          start_values = None
          for i in range(n_t):
              f[i], start_values = \
                  time_correlation_direct_by_mtx_vec_prod(P, pi, obs1, obs2,
                                                          times[i], start_values, True)
      return f