docs: def _opt_soft(eigvectors, rot_matrix, n_clusters):
      """
      Optimizes the PCCA+ rotation matrix such that the memberships are exclusively nonnegative.
  
      Parameters
      ----------
      eigenvectors : ndarray
          A matrix with the sorted eigenvectors in the columns. The stationary eigenvector should
          be first, then the one to the slowest relaxation process, etc.
  
      rot_mat : ndarray (m x m)
          nonoptimized rotation matrix
  
      n_clusters : int
          Number of clusters to group to.
  
      Returns
      -------
      rot_mat : ndarray (m x m)
          Optimized 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).
  
      """
      # only consider first n_clusters eigenvectors
      eigvectors = eigvectors[:, :n_clusters]
  
      # crop first row and first column from rot_matrix
      # rot_crop_matrix = rot_matrix[1:,1:]
      rot_crop_matrix = rot_matrix[1:][:, 1:]
  
      (x, y) = rot_crop_matrix.shape
  
      # reshape rot_crop_matrix into linear vector
      rot_crop_vec = np.reshape(rot_crop_matrix, x * y)
  
      # Susanna Roeblitz' target function for optimization
      def susanna_func(rot_crop_vec, eigvectors):
          # reshape into matrix
          rot_crop_matrix = np.reshape(rot_crop_vec, (x, y))
          # fill matrix
          rot_matrix = _fill_matrix(rot_crop_matrix, eigvectors)
  
          result = 0
          for i in range(0, n_clusters):
              for j in range(0, n_clusters):
                  result += np.power(rot_matrix[j, i], 2) / rot_matrix[0, i]
          return -result
  
      from scipy.optimize import fmin
  
      rot_crop_vec_opt = fmin(susanna_func, rot_crop_vec, args=(eigvectors,), disp=False)
  
      rot_crop_matrix = np.reshape(rot_crop_vec_opt, (x, y))
      rot_matrix = _fill_matrix(rot_crop_matrix, eigvectors)
  
      return rot_matrix