docs: def pcca(P, m):
      """
      PCCA+ spectral clustering method with optimized memberships [1]_
  
      Clusters the first m eigenvectors of a transition matrix in order to cluster the states.
      This function does not assume that the transition matrix is fully connected. Disconnected sets
      will automatically define the first metastable states, with perfect membership assignments.
  
      Parameters
      ----------
      P : ndarray (n,n)
          Transition matrix.
  
      m : 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.
  
      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).
      [2] F. Noe, multiset PCCA and HMMs, in preparation.
  
      """
      # imports
      from msmtools.estimation import connected_sets
      from msmtools.analysis import eigenvalues, is_transition_matrix, hitting_probability
  
      # validate input
      n = np.shape(P)[0]
      if (m > n):
          raise ValueError("Number of metastable states m = " + str(m)+
                           " exceeds number of states of transition matrix n = " + str(n))
      if not is_transition_matrix(P):
          raise ValueError("Input matrix is not a transition matrix.")
  
      # prepare output
      chi = np.zeros((n, m))
  
      # test connectivity
      components = connected_sets(P)
      # print "all labels ",labels
      n_components = len(components)  # (n_components, labels) = connected_components(P, connection='strong')
      # print 'n_components'
  
      # store components as closed (with positive equilibrium distribution)
      # or as transition states (with vanishing equilibrium distribution)
      closed_components = []
      transition_states = []
      for i in range(n_components):
          component = components[i]  # np.argwhere(labels==i).flatten()
          rest = list(set(range(n)) - set(component))
          # is component closed?
          if (np.sum(P[component, :][:, rest]) == 0):
              closed_components.append(component)
          else:
              transition_states.append(component)
      n_closed_components = len(closed_components)
      closed_states = np.concatenate(closed_components)
      if len(transition_states) == 0:
          transition_states = np.array([], dtype=int)
      else:
          transition_states = np.concatenate(transition_states)
  
      # check if we have enough clusters to support the disconnected sets
      if (m < len(closed_components)):
          raise ValueError("Number of metastable states m = " + str(m) + " is too small. Transition matrix has " +
                           str(len(closed_components)) + " disconnected components")
  
      # We collect eigenvalues in order to decide which
      closed_components_Psub = []
      closed_components_ev = []
      closed_components_enum = []
      for i in range(n_closed_components):
          component = closed_components[i]
          # print "component ",i," ",component
          # compute eigenvalues in submatrix
          Psub = P[component, :][:, component]
          closed_components_Psub.append(Psub)
          closed_components_ev.append(eigenvalues(Psub))
          closed_components_enum.append(i * np.ones((component.size), dtype=int))
  
      # flatten
      closed_components_ev_flat = np.hstack(closed_components_ev)
      closed_components_enum_flat = np.hstack(closed_components_enum)
      # which components should be clustered?
      component_indexes = closed_components_enum_flat[np.argsort(closed_components_ev_flat)][0:m]
      # cluster each component
      ipcca = 0
      for i in range(n_closed_components):
          component = closed_components[i]
          # how many PCCA states in this component?
          m_by_component = np.shape(np.argwhere(component_indexes == i))[0]
  
          # if 1, then the result is trivial
          if (m_by_component == 1):
              chi[component, ipcca] = 1.0
              ipcca += 1
          elif (m_by_component > 1):
              #print "submatrix: ",closed_components_Psub[i]
              chi[component, ipcca:ipcca + m_by_component] = _pcca_connected(closed_components_Psub[i], m_by_component)
              ipcca += m_by_component
          else:
              raise RuntimeError("Component " + str(i) + " spuriously has " + str(m_by_component) + " pcca sets")
  
      # finally assign all transition states
      # print "chi\n", chi
      # print "transition states: ",transition_states
      # print "closed states: ", closed_states
      if (transition_states.size > 0):
          # make all closed states absorbing, so we can see which closed state we hit first
          Pabs = P.copy()
          Pabs[closed_states, :] = 0.0
          Pabs[closed_states, closed_states] = 1.0
          for i in range(closed_states.size):
              # hitting probability to each closed state
              h = hitting_probability(Pabs, closed_states[i])
              for j in range(transition_states.size):
                  # transition states belong to closed states with the hitting probability, and inherit their chi
                  chi[transition_states[j]] += h[transition_states[j]] * chi[closed_states[i]]
  
      # check if we have m metastable sets. If less than m, we must raise
      nmeta = np.count_nonzero(chi.sum(axis=0))
      assert m <= nmeta, str(m) + " metastable states requested, but transition matrix only has " + str(nmeta) \
                         + ". Consider using a prior or request less metastable states. "
  
      # print "chi\n", chi
      return chi