docs: def _init_unit_centers(X, n_clusters, random_state, init):
      """Initializes unit norm centers.
  
      Parameters
      ----------
      X : array-like or sparse matrix, shape=(n_samples, n_features)
  
      n_clusters : int, optional, default: 8
          The number of clusters to form as well as the number of
          centroids to generate.
  
      random_state : integer or numpy.RandomState, optional
          The generator used to initialize the centers. If an integer is
          given, it fixes the seed. Defaults to the global numpy random
          number generator.
  
      init:  (string) one of
          k-means++ : uses sklearn k-means++ initialization algorithm
          spherical-k-means : use centroids from one pass of spherical k-means
          random : random unit norm vectors
          random-orthonormal : random orthonormal vectors
          If an ndarray is passed, it should be of shape (n_clusters, n_features)
          and gives the initial centers.
      """
      n_examples, n_features = np.shape(X)
      if isinstance(init, np.ndarray):
          n_init_clusters, n_init_features = init.shape
          assert n_init_clusters == n_clusters
          assert n_init_features == n_features
  
          # ensure unit normed centers
          centers = init
          for cc in range(n_clusters):
              centers[cc, :] = centers[cc, :] / np.linalg.norm(centers[cc, :])
  
          return centers
  
      elif init == "spherical-k-means":
          labels, inertia, centers, iters = spherical_kmeans._spherical_kmeans_single_lloyd(
              X, n_clusters, x_squared_norms=np.ones((n_examples,)), init="k-means++"
          )
  
          return centers
  
      elif init == "random":
          centers = np.random.randn(n_clusters, n_features)
          for cc in range(n_clusters):
              centers[cc, :] = centers[cc, :] / np.linalg.norm(centers[cc, :])
  
          return centers
  
      elif init == "k-means++":
          centers = _init_centroids(
              X,
              n_clusters,
              "k-means++",
              random_state=random_state,
              x_squared_norms=np.ones((n_examples,)),
          )
  
          for cc in range(n_clusters):
              centers[cc, :] = centers[cc, :] / np.linalg.norm(centers[cc, :])
  
          return centers
  
      elif init == "random-orthonormal":
          centers = np.random.randn(n_clusters, n_features)
          q, r = np.linalg.qr(centers.T, mode="reduced")
  
          return q.T
  
      elif init == "random-class":
          centers = np.zeros((n_clusters, n_features))
          for cc in range(n_clusters):
              while np.linalg.norm(centers[cc, :]) == 0:
                  labels = np.random.randint(0, n_clusters, n_examples)
                  centers[cc, :] = X[labels == cc, :].sum(axis=0)
  
          for cc in range(n_clusters):
              centers[cc, :] = centers[cc, :] / np.linalg.norm(centers[cc, :])
  
          return centers