docs: def _pcca_connected(P, n, return_rot=False):
      """
      PCCA+ spectral clustering method with optimized memberships [1]_
  
      Clusters the first n_cluster eigenvectors of a transition matrix in order to cluster the states.
      This function assumes that the transition matrix is fully connected.
  
      Parameters
      ----------
      P : ndarray (n,n)
          Transition matrix.
  
      n : int
          Number of clusters to group to.
  
      Returns
      -------
      chi by default, or (chi,rot) if return_rot = True
  
      chi : ndarray (n x m)
          A matrix containing the probability or membership of each state to be assigned to each cluster.
          The rows sum to 1.
  
      rot_mat : ndarray (m x m)
          A rotation matrix that rotates the dominant eigenvectors to yield the PCCA memberships, i.e.:
          chi = np.dot(evec, rot_matrix
  
      References
      ----------
      [1] S. Roeblitz and M. Weber, Fuzzy spectral clustering by PCCA+:
          application to Markov state models and data classification.
          Adv Data Anal Classif 7, 147-179 (2013).
  
      """
  
      # test connectivity
      from msmtools.estimation import connected_sets
  
      labels = connected_sets(P)
      n_components = len(labels)  # (n_components, labels) = connected_components(P, connection='strong')
      if (n_components > 1):
          raise ValueError("Transition matrix is disconnected. Cannot use pcca_connected.")
  
      from msmtools.analysis import stationary_distribution
  
      pi = stationary_distribution(P)
      # print "statdist = ",pi
  
      from msmtools.analysis import is_reversible
  
      if not is_reversible(P, mu=pi):
          raise ValueError("Transition matrix does not fulfill detailed balance. "
                           "Make sure to call pcca with a reversible transition matrix estimate")
      # TODO: Susanna mentioned that she has a potential fix for nonreversible matrices by replacing each complex conjugate
      #      pair by the real and imaginary components of one of the two vectors. We could use this but would then need to
      #      orthonormalize all eigenvectors e.g. using Gram-Schmidt orthonormalization. Currently there is no theoretical
      #      foundation for this, so I'll skip it for now.
  
      # right eigenvectors, ordered
      from msmtools.analysis import eigenvectors
  
      evecs = eigenvectors(P, n)
  
      # orthonormalize
      for i in range(n):
          evecs[:, i] /= math.sqrt(np.dot(evecs[:, i] * pi, evecs[:, i]))
      # make first eigenvector positive
      evecs[:, 0] = np.abs(evecs[:, 0])
  
      # Is there a significant complex component?
      if not np.alltrue(np.isreal(evecs)):
          warnings.warn(
              "The given transition matrix has complex eigenvectors, so it doesn't exactly fulfill detailed balance "
              + "forcing eigenvectors to be real and continuing. Be aware that this is not theoretically solid.")
      evecs = np.real(evecs)
  
      # create initial solution using PCCA+. This could have negative memberships
      (chi, rot_matrix) = _pcca_connected_isa(evecs, n)
  
      #print "initial chi = \n",chi
  
      # optimize the rotation matrix with PCCA++.
      rot_matrix = _opt_soft(evecs, rot_matrix, n)
  
      # These memberships should be nonnegative
      memberships = np.dot(evecs[:, :], rot_matrix)
  
      # We might still have numerical errors. Force memberships to be in [0,1]
      # print "memberships unnormalized: ",memberships
      memberships = np.maximum(0.0, memberships)
      memberships = np.minimum(1.0, memberships)
      # print "memberships unnormalized: ",memberships
      for i in range(0, np.shape(memberships)[0]):
          memberships[i] /= np.sum(memberships[i])
  
      # print "final chi = \n",chi
  
      return memberships