code: def init_cao(X, n_clusters, dissim):
      
      n_points, n_attrs = X.shape
      centroids = np.empty((n_clusters, n_attrs), dtype='object')
      # Method is based on determining density of points.
      dens = np.zeros(n_points)
      for iattr in range(n_attrs):
          freq = defaultdict(int)
          for val in X[:, iattr]:
              freq[val] += 1
          for ipoint in range(n_points):
              dens[ipoint] += freq[X[ipoint, iattr]] / float(n_points) / float(n_attrs)
  
      # Choose initial centroids based on distance and density.
      centroids[0] = X[np.argmax(dens)]
      if n_clusters > 1:
          # For the remaining centroids, choose maximum dens * dissim to the
          # (already assigned) centroid with the lowest dens * dissim.
          for ik in range(1, n_clusters):
              dd = np.empty((ik, n_points))
              for ikk in range(ik):
                  dd[ikk] = dissim(X, centroids[ikk]) * dens
              centroids[ik] = X[np.argmax(np.min(dd, axis=0))]
  
      return centroids
docs: Initialize centroids according to method by Cao et al. [2009].
  
      Note: O(N * attr * n_clusters**2), so watch out with large n_clusters