docs: def transition_matrix_reversible_pisym(C, return_statdist=False, **kwargs):
      r"""
      Estimates reversible transition matrix as follows:
  
      ..:math:
          p_{ij} = c_{ij} / c_i where c_i = sum_j c_{ij}
          \pi_j = \sum_j \pi_i p_{ij}
          x_{ij} = \pi_i p_{ij} + \pi_j p_{ji}
          p^{rev}_{ij} = x_{ij} / x_i where x_i = sum_j x_{ij}
  
      In words: takes the nonreversible transition matrix estimate, uses its
      stationary distribution to compute an equilibrium correlation matrix,
      symmetrizes that correlation matrix and then normalizes to the reversible
      transition matrix estimate.
  
      Parameters
      ----------
      C: ndarray, shape (n,n)
          count matrix
  
      Returns
      -------
      T: Estimated transition matrix
  
      """
      # nonreversible estimate
      T_nonrev = transition_matrix_non_reversible(C)
      from msmtools.analysis import stationary_distribution
      pi = stationary_distribution(T_nonrev)
      # correlation matrix
      X = scipy.sparse.diags(pi).dot(T_nonrev)
      X = X.T + X
      # result
      pi_rev = np.array(X.sum(axis=1)).squeeze()
      T_rev = scipy.sparse.diags(1.0/pi_rev).dot(X)
      if return_statdist:
          #np.testing.assert_allclose(pi, stationary_distribution(T_rev))
          #np.testing.assert_allclose(T_rev.T.dot(pi), pi)
          return T_rev, pi
      return T_rev