Repository: martinfleis/clustergram

Issue:
Optionally measure silhouette and other metrics
Add an option to measure additional metrics to assess the results of clustering, like a silhouette score or Calinski-Harabasz.

Discussion:

Patch:
diff --git a/README.md b/README.md
index 37fa78e..2e70303 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@ Clustergram is a diagram proposed by Matthias Schonlau in his paper *[The cluste
 
 The clustergram was later implemented in R by [Tal Galili](https://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/), who also gives a thorough explanation of the concept.
 
-This is a Python translation of Tal's script written for `scikit-learn` and RAPIDS `cuML` implementations of K-Means and Gaussian Mixture Model (scikit-learn only) clustering.
+This is a Python translation of Tal's script written for `scikit-learn` and RAPIDS `cuML` implementations of K-Means, Mini Batch K-Means and Gaussian Mixture Model (scikit-learn only) clustering.
 
 ## Getting started
 
@@ -24,7 +24,7 @@ conda install clustergram -c conda-forge
 pip install clustergram
 ```
 
-In any case, you still need to install your selected backend 
+In any case, you still need to install your selected backend
 (`scikit-learn` or `cuML`).
 
 The example of clustergram on Palmer penguins dataset:
@@ -76,16 +76,16 @@ cgram.plot(
 On the `y` axis, a clustergram can use mean values as in the original paper by Matthias Schonlau or PCA weighted mean values as in the implementation by Tal Galili.
 
 ```python
-cgram = Clustergram(range(1, 8), pca_weighted=True)
+cgram = Clustergram(range(1, 8))
 cgram.fit(data)
-cgram.plot(figsize=(12, 8))
+cgram.plot(figsize=(12, 8), pca_weighted=True)
 ```
 ![Default clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/pca_true.png)
 
 ```python
-cgram = Clustergram(range(1, 8), pca_weighted=False)
+cgram = Clustergram(range(1, 8))
 cgram.fit(data)
-cgram.plot(figsize=(12, 8))
+cgram.plot(figsize=(12, 8), pca_weighted=False)
 ```
 ![Default clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/pca_false.png)
 
@@ -93,7 +93,7 @@ cgram.plot(figsize=(12, 8))
 
 Clustergram offers two backends for the computation - `scikit-learn` which uses CPU and RAPIDS.AI `cuML`, which uses GPU. Note that both are optional dependencies, but you will need at least one of them to generate clustergram.
 
-Using scikit-learn (default):
+Using `scikit-learn` (default):
 
 ```python
 cgram = Clustergram(range(1, 8), backend='sklearn')
@@ -101,7 +101,7 @@ cgram.fit(data)
 cgram.plot()
 ```
 
-Using cuML:
+Using `cuML`:
 
 ```python
 cgram = Clustergram(range(1, 8), backend='cuML')
@@ -113,7 +113,7 @@ cgram.plot()
 
 ## Supported methods
 
-Clustergram currently supports K-Means and Gaussian Mixture Model clustering methods. Note tha GMM and Mini Batch K-Means are supported only for `scikit-learn` backend.
+Clustergram currently supports K-Means, Mini Batch K-Means and Gaussian Mixture Model clustering methods. Note tha GMM and Mini Batch K-Means are supported only for `scikit-learn` backend.
 
 Using K-Means (default):
 
@@ -126,7 +126,7 @@ cgram.plot()
 Using Mini Batch K-Means, which can provide significant speedup over K-Means:
 
 ```python
-cgram = Clustergram(range(1, 8), method='minibatchkmeans')
+cgram = Clustergram(range(1, 8), method='minibatchkmeans', batch_size=100)
 cgram.fit(data)
 cgram.plot()
 ```
@@ -155,6 +155,79 @@ cgram.plot(k_range=range(3, 10), figsize=(12, 8))
 ```
 ![Limited clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/limited_plot.png)
 
+## Additional clustering performance evaluation
+
+Clustergam includes handy wrappers around a selection of clustering performance metrics offered by
+`scikit-learn`. Data which were originally computed on GPU are converted to numpy on the fly.
+
+### Silhouette score
+
+Compute the mean Silhouette Coefficient of all samples. See [`scikit-learn` documentation](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html#sklearn.metrics.silhouette_score) for details.
+
+```python
+>>> cgram.silhouette_score()
+2    0.531540
+3    0.447219
+4    0.400154
+5    0.377720
+6    0.372128
+7    0.331575
+Name: silhouette_score, dtype: float64
+```
+Once computed, resulting Series is available as `cgram.silhouette`. Calling the original method will recompute the score.
+
+### Calinski and Harabasz score
+
+Compute the Calinski and Harabasz score, also known as the Variance Ratio Criterion. See [`scikit-learn` documentation](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.calinski_harabasz_score.html#sklearn.metrics.calinski_harabasz_score) for details.
+
+```python
+>>> cgram.calinski_harabasz_score()
+2    482.191469
+3    441.677075
+4    400.392131
+5    411.175066
+6    382.731416
+7    352.447569
+Name: calinski_harabasz_score, dtype: float64
+```
+Once computed, resulting Series is available as `cgram.calinski_harabasz`. Calling the original method will recompute the score.
+
+### Davies-Bouldin score
+
+Compute the Davies-Bouldin score. See [`scikit-learn` documentation](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.davies_bouldin_score.html#sklearn.metrics.davies_bouldin_score) for details.
+
+```python
+>>> cgram.davies_bouldin_score()
+2    0.714064
+3    0.943553
+4    0.943320
+5    0.973248
+6    0.950910
+7    1.074937
+Name: davies_bouldin_score, dtype: float64
+```
+Once computed, resulting Series is available as `cgram.davies_bouldin`. Calling the original method will recompute the score.
+
+## Acessing labels
+
+`Clustergram` stores resulting labels for each of the tested options, which can be accessed as:
+
+```python
+>>> cgram.labels
+     1  2  3  4  5  6  7
+0    0  0  2  2  3  2  1
+1    0  0  2  2  3  2  1
+2    0  0  2  2  3  2  1
+3    0  0  2  2  3  2  1
+4    0  0  2  2  0  0  3
+..  .. .. .. .. .. .. ..
+337  0  1  1  3  2  5  0
+338  0  1  1  3  2  5  0
+339  0  1  1  1  1  1  4
+340  0  1  1  3  2  5  5
+341  0  1  1  1  1  1  5
+```
+
 ## Saving clustergram
 
 You can save both plot and `clustergram.Clustergram` to a disk.
diff --git a/clustergram/clustergram.py b/clustergram/clustergram.py
index 0baa9ae..8f6f362 100644
--- a/clustergram/clustergram.py
+++ b/clustergram/clustergram.py
@@ -1,7 +1,7 @@
 """
 Clustergram - visualization and diagnostics for cluster analysis in Python
 
-Copyright (C) 2020  Martin Fleischmann
+Copyright (C) 2020-2021  Martin Fleischmann
 
 Clustergram is a Python implementation of R function written by Tal Galili.
 https://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/
@@ -11,6 +11,8 @@ Original idea is by Matthias Schonlau - http://www.schonlau.net/clustergram.html
 """
 
 from time import time
+import pandas as pd
+import numpy as np
 
 
 class Clustergram:
@@ -31,19 +33,22 @@ class Clustergram:
     Parameters
     ----------
     k_range : iterable
-        iterable of integer values to be tested as k.
+        iterable of integer values to be tested as ``k`` (number of cluster or
+        components).
     backend : {'sklearn', 'cuML'} (default 'sklearn')
-        Whether to use ``sklearn``'s implementation of KMeans and PCA or ``cuML`` version.
-        ``sklearn`` does computation on CPU, ``cuml`` on GPU.
+        Whether to use ``sklearn``'s implementation of KMeans and PCA or ``cuML``
+        version. ``sklearn`` does computation on CPU, ``cuml`` on GPU.
     method : {'kmeans', 'gmm', 'minibatchkmeans'} (default 'kmeans')
-        Clustering method. ``kmeans`` uses K-Means clustering, ``gmm`` Gaussian Mixture Model,
-        ``minibatchkmeans`` uses Mini Batch K-Means.
-        ``gmm`` and ``minibatchkmeans`` are currently supported only with ``sklearn`` backend.
-    pca_weighted : bool (default True)
-        Whether use PCA weighted mean of clusters or standard mean of clusters.
-    pca_kwargs : dict (default {})
-        Additional arguments passed to the PCA object,
-        e.g. ``svd_solver``. Applies only if ``pca_weighted=True``.
+        Clustering method.
+
+        * ``kmeans`` uses K-Means clustering, either as ``sklearn.cluster.KMeans``
+          or ``cuml.KMeans``.
+        * ``gmm`` uses Gaussian Mixture Model as ``sklearn.mixture.GaussianMixture``
+        * ``minibatchkmeans`` uses Mini Batch K-Means as
+          ``sklearn.cluster.MiniBatchKMeans``
+
+        Note that ``gmm`` and ``minibatchkmeans`` are currently supported only
+        with ``sklearn`` backend.
     verbose : bool (default True)
         Print progress and time of individual steps.
     **kwargs
@@ -54,8 +59,11 @@ class Clustergram:
     Attributes
     ----------
 
-    means : DataFrame
-        DataFrame with (weighted) means of clusters.
+    labels : DataFrame
+        DataFrame with cluster labels for each iteration.
+    cluster_centers : dict
+        Dictionary with cluster centers for each iteration.
+
 
     Examples
     --------
@@ -66,7 +74,7 @@ class Clustergram:
     Specifying parameters:
 
     >>> c_gram2 = clustergram.Clustergram(
-    ...     range(1, 9), backend="cuML", pca_weighted=False, random_state=0
+    ...     range(1, 9), backend="cuML", random_state=0
     ... )
     >>> c_gram2.fit(cudf_data)
     >>> c_gram2.plot(figsize=(12, 12))
@@ -83,17 +91,14 @@ class Clustergram:
     """
 
     def __init__(
-        self,
-        k_range,
-        backend="sklearn",
-        method="kmeans",
-        pca_weighted=True,
-        pca_kwargs={},
-        verbose=True,
-        **kwargs,
+        self, k_range, backend="sklearn", method="kmeans", verbose=True, **kwargs,
     ):
         self.k_range = k_range
 
+        # cleanup after API change
+        kwargs.pop("pca_weighted", None)
+        kwargs.pop("pca_kwargs", None)
+
         if backend not in ["sklearn", "cuML"]:
             raise ValueError(
                 f'"{backend}" is not a supported backend. Use "sklearn" or "cuML".'
@@ -104,25 +109,33 @@ class Clustergram:
         supported = ["kmeans", "gmm", "minibatchkmeans"]
         if method not in supported:
             raise ValueError(
-                f'"{method}" is not a supported method. Only {supported} are supported now.'
+                f"'{method}' is not a supported method. "
+                f"Only {supported} are supported now."
             )
         else:
             self.method = method
 
-        self.pca_weighted = pca_weighted
         self.engine_kwargs = kwargs
-        self.pca_kwargs = pca_kwargs
         self.verbose = verbose
 
+        if self.backend == "sklearn":
+            self.plot_data_pca = pd.DataFrame()
+            self.plot_data = pd.DataFrame()
+        else:
+            import cudf
+
+            self.plot_data_pca = cudf.DataFrame()
+            self.plot_data = cudf.DataFrame()
+
     def __repr__(self):
         return (
             f"Clustergram(k_range={self.k_range}, backend='{self.backend}', "
-            f"method='{self.method}', pca_weighted={self.pca_weighted}, kwargs={self.engine_kwargs})"
+            f"method='{self.method}', kwargs={self.engine_kwargs})"
         )
 
     def fit(self, data, **kwargs):
         """
-        Compute (weighted) means of clusters.
+        Compute clustering for each k within set range.
 
         Parameters
         ----------
@@ -138,12 +151,18 @@ class Clustergram:
         self
             Fitted clustergram.
 
+        Examples
+        --------
+        >>> c_gram = clustergram.Clustergram(range(1, 9))
+        >>> c_gram.fit(data)
+
         """
+        self.data = data
         if self.backend == "sklearn":
             if self.method == "kmeans":
-                self.means = self._kmeans_sklearn(data, minibatch=False, **kwargs)
+                self._kmeans_sklearn(data, minibatch=False, **kwargs)
             elif self.method == "minibatchkmeans":
-                self.means = self._kmeans_sklearn(data, minibatch=True, **kwargs)
+                self._kmeans_sklearn(data, minibatch=True, **kwargs)
             elif self.method == "gmm":
                 self.means = self._gmm_sklearn(data, **kwargs)
         if self.backend == "cuML":
@@ -153,20 +172,11 @@ class Clustergram:
         """Use scikit-learn KMeans"""
         try:
             from sklearn.cluster import KMeans, MiniBatchKMeans
-            from sklearn.decomposition import PCA
-            from pandas import DataFrame
-            import numpy as np
         except ImportError:
-            raise ImportError(
-                "scikit-learn, pandas and numpy are required to use `sklearn` backend."
-            )
+            raise ImportError("scikit-learn is required to use `sklearn` backend.")
 
-        df = DataFrame()
-        if self.pca_weighted:
-            s = time()
-            self.pca_kwargs.pop("n_components", 1)
-            pca = PCA(n_components=1, **self.pca_kwargs).fit(data)
-            print(f"PCA computed in {time() - s} seconds.") if self.verbose else None
+        self.labels = pd.DataFrame()
+        self.cluster_centers = {}
 
         for n in self.k_range:
             s = time()
@@ -176,85 +186,56 @@ class Clustergram:
                 )
             else:
                 results = KMeans(n_clusters=n, **self.engine_kwargs).fit(data, **kwargs)
-            cluster = results.labels_
-            if self.pca_weighted:
-                means = results.cluster_centers_.dot(pca.components_[0])
-            else:
-                means = np.mean(results.cluster_centers_, axis=1)
-            df[n] = np.take(means, cluster)
+
+            self.labels[n] = results.labels_
+            self.cluster_centers[n] = results.cluster_centers_
+
             print(f"K={n} fitted in {time() - s} seconds.") if self.verbose else None
-        return df
 
     def _kmeans_cuml(self, data, **kwargs):
         """Use cuML KMeans"""
         try:
-            from cuml import KMeans, PCA
-            from cudf import DataFrame
-            import cupy as cp
-            import numpy as np
+            from cuml import KMeans
+            import cudf
         except ImportError:
             raise ImportError(
                 "cuML, cuDF and cupy packages are required to use `cuML` backend."
             )
 
-        df = DataFrame()
-        if self.pca_weighted:
-            s = time()
-            self.pca_kwargs.pop("n_components", 1)
-            pca = PCA(n_components=1, **self.pca_kwargs).fit(data)
-            print(f"PCA computed in {time() - s} seconds.") if self.verbose else None
+        self.labels = cudf.DataFrame()
+        self.cluster_centers = {}
 
         for n in self.k_range:
             s = time()
             results = KMeans(n_clusters=n, **self.engine_kwargs).fit(data, **kwargs)
-            cluster = results.labels_
-            if self.pca_weighted:
-                if isinstance(results.cluster_centers_, DataFrame):
-                    means = results.cluster_centers_.values.dot(
-                        pca.components_.values[0]
-                    )
-                else:
-                    means = results.cluster_centers_.dot(pca.components_[0])
-                df[n] = cp.take(means, cluster)
-            else:
-                means = results.cluster_centers_.mean(axis=1)
-                if isinstance(means, (cp.core.core.ndarray, np.ndarray)):
-                    df[n] = means.take(cluster)
-                else:
-                    df[n] = means.take(cluster).to_array()
+            self.labels[n] = results.labels_
+            self.cluster_centers[n] = results.cluster_centers_
+
             print(f"K={n} fitted in {time() - s} seconds.") if self.verbose else None
-        return df
 
     def _gmm_sklearn(self, data, **kwargs):
         """Use sklearn.mixture.GaussianMixture"""
         try:
             from sklearn.mixture import GaussianMixture
-            from sklearn.decomposition import PCA
-            from pandas import DataFrame
             import numpy as np
             from scipy.stats import multivariate_normal
         except ImportError:
             raise ImportError(
-                "scikit-learn, pandas and numpy are required to use `sklearn` backend."
+                "scikit-learn and scipy are required to use `sklearn` "
+                "backend and `gmm`."
             )
 
-        if isinstance(data, DataFrame):
+        if isinstance(data, pd.DataFrame):
             data = data.values
 
-        df = DataFrame()
-
-        if self.pca_weighted:
-            s = time()
-            self.pca_kwargs.pop("n_components", 1)
-            pca = PCA(n_components=1, **self.pca_kwargs).fit(data)
-            print(f"PCA computed in {time() - s} seconds.") if self.verbose else None
+        self.labels = pd.DataFrame()
+        self.cluster_centers = {}
 
         for n in self.k_range:
             s = time()
             results = GaussianMixture(n_components=n, **self.engine_kwargs).fit(
                 data, **kwargs
             )
-            cluster = results.predict(data)
             centers = np.empty(shape=(results.n_components, data.shape[1]))
             for i in range(results.n_components):
                 density = multivariate_normal(
@@ -263,13 +244,266 @@ class Clustergram:
                     allow_singular=True,
                 ).logpdf(data)
                 centers[i, :] = data[np.argmax(density)]
-            if self.pca_weighted:
-                means = centers.dot(pca.components_[0])
-            else:
-                means = np.mean(centers, axis=1)
-            df[n] = np.take(means, cluster)
+
+            self.labels[n] = results.predict(data)
+            self.cluster_centers[n] = centers
+
             print(f"K={n} fitted in {time() - s} seconds.") if self.verbose else None
-        return df
+
+    def silhouette_score(self, **kwargs):
+        """
+        Compute the mean Silhouette Coefficient of all samples.
+
+        See the documentation of ``sklearn.metrics.silhouette_score`` for details.
+
+        Once computed, resulting Series is available as ``Clustergram.silhouette``.
+        Calling the original method will compute the score from the beginning.
+
+        Parameters
+        ----------
+
+        **kwargs
+            Additional arguments passed to the silhouette_score function,
+            e.g. ``sample_size``.
+
+        Returns
+        -------
+        self.silhouette : pd.Series
+
+        Notes
+        -----
+        The algortihm uses ``sklearn``.
+        With ``cuML`` backend, data are converted on the fly.
+
+        Examples
+        --------
+        >>> c_gram = clustergram.Clustergram(range(1, 9))
+        >>> c_gram.fit(data)
+        >>> c_gram.silhouette_score()
+        2    0.702450
+        3    0.644272
+        4    0.767728
+        5    0.948991
+        6    0.769985
+        7    0.575644
+        Name: silhouette_score, dtype: float64
+
+        Once computed:
+
+        >>> c_gram.silhouette
+        2    0.702450
+        3    0.644272
+        4    0.767728
+        5    0.948991
+        6    0.769985
+        7    0.575644
+        Name: silhouette_score, dtype: float64
+
+        """
+        from sklearn import metrics
+
+        self.silhouette = pd.Series(name="silhouette_score", dtype="float64")
+
+        if self.backend == "sklearn":
+            for k in self.k_range:
+                if k > 1:
+                    self.silhouette.loc[k] = metrics.silhouette_score(
+                        self.data, self.labels[k], **kwargs
+                    )
+        else:
+            if hasattr(self.data, "to_pandas"):
+                data = self.data.to_pandas()
+            else:
+                data = self.data.get()
+
+            for k in self.k_range:
+                if k > 1:
+                    self.silhouette.loc[k] = metrics.silhouette_score(
+                        data, self.labels[k].to_pandas(), **kwargs
+                    )
+
+        return self.silhouette
+
+    def calinski_harabasz_score(self):
+        """
+        Compute the Calinski and Harabasz score.
+
+        See the documentation of ``sklearn.metrics.calinski_harabasz_score``
+        for details.
+
+        Once computed, resulting Series is available as
+        ``Clustergram.calinski_harabasz``. Calling the original method will
+        compute the score from the beginning.
+
+        Returns
+        -------
+        self.calinski_harabasz : pd.Series
+
+        Notes
+        -----
+        The algortihm uses ``sklearn``.
+        With ``cuML`` backend, data are converted on the fly.
+
+        Examples
+        --------
+        >>> c_gram = clustergram.Clustergram(range(1, 9))
+        >>> c_gram.fit(data)
+        >>> c_gram.calinski_harabasz_score()
+        2      23.176629
+        3      30.643018
+        4      55.223336
+        5    3116.435184
+        6    3899.068689
+        7    4439.306049
+        Name: calinski_harabasz_score, dtype: float64
+
+        Once computed:
+
+        >>> c_gram.calinski_harabasz
+        2      23.176629
+        3      30.643018
+        4      55.223336
+        5    3116.435184
+        6    3899.068689
+        7    4439.306049
+        Name: calinski_harabasz_score, dtype: float64
+
+        """
+        from sklearn import metrics
+
+        self.calinski_harabasz = pd.Series(
+            name="calinski_harabasz_score", dtype="float64"
+        )
+
+        if self.backend == "sklearn":
+            for k in self.k_range:
+                if k > 1:
+                    self.calinski_harabasz.loc[k] = metrics.calinski_harabasz_score(
+                        self.data, self.labels[k]
+                    )
+        else:
+            if hasattr(self.data, "to_pandas"):
+                data = self.data.to_pandas()
+            else:
+                data = self.data.get()
+
+            for k in self.k_range:
+                if k > 1:
+                    self.calinski_harabasz.loc[k] = metrics.calinski_harabasz_score(
+                        data, self.labels[k].to_pandas()
+                    )
+        return self.calinski_harabasz
+
+    def davies_bouldin_score(self):
+        """
+        Compute the Davies-Bouldin score.
+
+        See the documentation of ``sklearn.metrics.davies_bouldin_score`` for details.
+
+        Once computed, resulting Series is available as ``Clustergram.davies_bouldin``.
+        Calling the original method will recompute the score.
+
+        Returns
+        -------
+        self.davies_bouldin : pd.Series
+
+        Notes
+        -----
+        The algortihm uses ``sklearn``.
+        With ``cuML`` backend, data are converted on the fly.
+
+        Examples
+        --------
+        >>> c_gram = clustergram.Clustergram(range(1, 9))
+        >>> c_gram.fit(data)
+        >>> c_gram.davies_bouldin_score()
+        2    0.249366
+        3    0.351812
+        4    0.347580
+        5    0.055679
+        6    0.030516
+        7    0.025207
+        Name: davies_bouldin_score, dtype: float64
+
+        Once computed:
+
+        >>> c_gram.davies_bouldin
+        2    0.249366
+        3    0.351812
+        4    0.347580
+        5    0.055679
+        6    0.030516
+        7    0.025207
+        Name: davies_bouldin_score, dtype: float64
+
+        """
+        from sklearn import metrics
+
+        self.davies_bouldin = pd.Series(name="davies_bouldin_score", dtype="float64")
+
+        if self.backend == "sklearn":
+            for k in self.k_range:
+                if k > 1:
+                    self.davies_bouldin.loc[k] = metrics.davies_bouldin_score(
+                        self.data, self.labels[k]
+                    )
+        else:
+            if hasattr(self.data, "to_pandas"):
+                data = self.data.to_pandas()
+            else:
+                data = self.data.get()
+
+            for k in self.k_range:
+                if k > 1:
+                    self.davies_bouldin.loc[k] = metrics.davies_bouldin_score(
+                        data, self.labels[k].to_pandas()
+                    )
+
+        return self.davies_bouldin
+
+    def _compute_pca_means_sklearn(self, **pca_kwargs):
+        """Compute PCA weighted cluster mean values using sklearn backend"""
+        from sklearn.decomposition import PCA
+
+        self.pca = PCA(n_components=1, **pca_kwargs).fit(self.data).components_[0]
+
+        for n in self.k_range:
+            means = self.cluster_centers[n].dot(self.pca)
+            self.plot_data_pca[n] = np.take(means, self.labels[n].values)
+
+    def _compute_means_sklearn(self):
+        """Compute cluster mean values using sklearn backend"""
+        for n in self.k_range:
+            means = np.mean(self.cluster_centers[n], axis=1)
+            self.plot_data[n] = np.take(means, self.labels[n].values)
+
+    def _compute_pca_means_cuml(self, **pca_kwargs):
+        """Compute PCA weighted cluster mean values using cuML backend"""
+        from cuml import PCA
+        import cudf
+        import cupy as cp
+
+        self.pca = PCA(n_components=1, **pca_kwargs).fit(self.data)
+
+        for n in self.k_range:
+            if isinstance(self.data, cudf.DataFrame):
+                means = self.cluster_centers[n].values.dot(
+                    self.pca.components_.values[0]
+                )
+            else:
+                means = self.cluster_centers[n].dot(self.pca.components_[0])
+            self.plot_data_pca[n] = cp.take(means, self.labels[n].values)
+
+    def _compute_means_cuml(self):
+        """Compute cluster mean values using cuML backend"""
+        import cupy as cp
+
+        for n in self.k_range:
+            means = self.cluster_centers[n].mean(axis=1)
+            if isinstance(means, (cp.core.core.ndarray, np.ndarray)):
+                self.plot_data[n] = means.take(self.labels[n].values)
+            else:
+                self.plot_data[n] = means.take(self.labels[n].values).to_array()
 
     def plot(
         self,
@@ -280,6 +514,8 @@ class Clustergram:
         line_style=None,
         figsize=None,
         k_range=None,
+        pca_weighted=True,
+        pca_kwargs={},
     ):
         """
         Generate clustergram plot based on cluster centre mean values.
@@ -306,11 +542,43 @@ class Clustergram:
         k_range : iterable (default None)
             iterable of integer values to be plotted. In none, ``Clustergram.k_range``
             will be used. Has to be a subset of ``Clustergram.k_range``.
+        pca_weighted : bool (default True)
+            Whether use PCA weighted mean of clusters or standard mean of clusters on
+            y-axis.
+        pca_kwargs : dict (default {})
+            Additional arguments passed to the PCA object,
+            e.g. ``svd_solver``. Applies only if ``pca_weighted=True``.
 
         Returns
         -------
         ax : matplotlib axis instance
+
+        Examples
+        --------
+        >>> c_gram = clustergram.Clustergram(range(1, 9))
+        >>> c_gram.fit(data)
+        >>> c_gram.plot()
+
+        Notes
+        -----
+        Before plotting, ``Clustergram`` needs to compute the summary values.
+        Those are computed on the first call of each option (pca_weighted=True/False).
+
         """
+        if pca_weighted:
+            if self.plot_data_pca.empty:
+                pca_kwargs.pop("n_components", None)
+
+                if self.backend == "sklearn":
+                    self._compute_pca_means_sklearn(**pca_kwargs)
+                else:
+                    self._compute_pca_means_cuml(**pca_kwargs)
+        else:
+            if self.plot_data.empty:
+                if self.backend == "sklearn":
+                    self._compute_means_sklearn()
+                else:
+                    self._compute_means_cuml()
 
         if ax is None:
             import matplotlib.pyplot as plt
@@ -333,14 +601,22 @@ class Clustergram:
         if k_range is None:
             k_range = self.k_range
 
+        if pca_weighted:
+            means = self.plot_data_pca
+            ax.set_ylabel("PCA weighted mean of the clusters")
+        else:
+            means = self.plot_data
+            ax.set_ylabel("Mean of the clusters")
+        ax.set_xlabel("Number of clusters (k)")
+
         for i in k_range:
-            cl = self.means[i].value_counts()
+            cl = means[i].value_counts()
 
             if self.backend == "sklearn":
                 ax.scatter(
                     [i] * i,
                     [cl.index],
-                    cl * ((500 / len(self.means)) * size),
+                    cl * ((500 / len(means)) * size),
                     zorder=cl_zorder,
                     color=cl_c,
                     edgecolor=cl_ec,
@@ -351,7 +627,7 @@ class Clustergram:
                 ax.scatter(
                     [i] * i,
                     cl.index.to_array(),
-                    (cl * ((500 / len(self.means)) * size)).to_array(),
+                    (cl * ((500 / len(means)) * size)).to_array(),
                     zorder=cl_zorder,
                     color=cl_c,
                     edgecolor=cl_ec,
@@ -361,16 +637,14 @@ class Clustergram:
 
             try:
                 if self.backend == "sklearn":
-                    sub = self.means.groupby([i, i + 1]).count().reset_index()
+                    sub = means.groupby([i, i + 1]).count().reset_index()
                 else:
-                    sub = (
-                        self.means.groupby([i, i + 1]).count().reset_index().to_pandas()
-                    )
+                    sub = means.groupby([i, i + 1]).count().reset_index().to_pandas()
                 for r in sub.itertuples():
                     ax.plot(
                         [i, i + 1],
                         [r[1], r[2]],
-                        linewidth=r[3] * ((50 / len(self.means)) * linewidth),
+                        linewidth=r[3] * ((50 / len(means)) * linewidth),
                         color=l_c,
                         zorder=l_zorder,
                         solid_capstyle=solid_capstyle,
@@ -378,9 +652,4 @@ class Clustergram:
                     )
             except (KeyError, ValueError):
                 pass
-        if self.pca_weighted:
-            ax.set_ylabel("PCA weighted mean of the clusters")
-        else:
-            ax.set_ylabel("Mean of the clusters")
-        ax.set_xlabel("Number of clusters (k)")
         return ax
diff --git a/doc/api.rst b/doc/api.rst
index 02fca57..63d446b 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -10,4 +10,4 @@ clustergram
 .. currentmodule:: clustergram
 
 .. autoclass:: Clustergram
-   :members: fit, plot
\ No newline at end of file
+   :members: fit, plot, silhouette_score, calinski_harabasz_score, davies_bouldin_score
\ No newline at end of file
diff --git a/doc/index.md b/doc/index.md
index 718e97d..7d6fd98 100644
--- a/doc/index.md
+++ b/doc/index.md
@@ -8,7 +8,7 @@ Clustergram is a diagram proposed by Matthias Schonlau in his paper *[The cluste
 
 The clustergram was later implemented in R by [Tal Galili](https://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/), who also gives a thorough explanation of the concept.
 
-This is a Python translation of Tal's script written for `scikit-learn` and RAPIDS `cuML` implementations of K-Means and Gaussian Mixture Model (scikit-learn only) clustering.
+This is a Python translation of Tal's script written for `scikit-learn` and RAPIDS `cuML` implementations of K-Means, Mini Batch K-Means and Gaussian Mixture Model (scikit-learn only) clustering.
 
 ## Getting started
 
@@ -22,7 +22,7 @@ conda install clustergram -c conda-forge
 pip install clustergram
 ```
 
-In any case, you still need to install your selected backend 
+In any case, you still need to install your selected backend
 (`scikit-learn` or `cuML`).
 
 The example of clustergram on Palmer penguins dataset:
@@ -74,16 +74,16 @@ cgram.plot(
 On the `y` axis, a clustergram can use mean values as in the original paper by Matthias Schonlau or PCA weighted mean values as in the implementation by Tal Galili.
 
 ```python
-cgram = Clustergram(range(1, 8), pca_weighted=True)
+cgram = Clustergram(range(1, 8))
 cgram.fit(data)
-cgram.plot(figsize=(12, 8))
+cgram.plot(figsize=(12, 8), pca_weighted=True)
 ```
 ![Default clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/pca_true.png)
 
 ```python
-cgram = Clustergram(range(1, 8), pca_weighted=False)
+cgram = Clustergram(range(1, 8))
 cgram.fit(data)
-cgram.plot(figsize=(12, 8))
+cgram.plot(figsize=(12, 8), pca_weighted=False)
 ```
 ![Default clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/pca_false.png)
 
@@ -91,7 +91,7 @@ cgram.plot(figsize=(12, 8))
 
 Clustergram offers two backends for the computation - `scikit-learn` which uses CPU and RAPIDS.AI `cuML`, which uses GPU. Note that both are optional dependencies, but you will need at least one of them to generate clustergram.
 
-Using scikit-learn (default):
+Using `scikit-learn` (default):
 
 ```python
 cgram = Clustergram(range(1, 8), backend='sklearn')
@@ -99,7 +99,7 @@ cgram.fit(data)
 cgram.plot()
 ```
 
-Using cuML:
+Using `cuML`:
 
 ```python
 cgram = Clustergram(range(1, 8), backend='cuML')
@@ -111,7 +111,7 @@ cgram.plot()
 
 ## Supported methods
 
-Clustergram currently supports K-Means and Gaussian Mixture Model clustering methods. Note tha GMM and Mini Batch K-Means are supported only for `scikit-learn` backend.
+Clustergram currently supports K-Means, Mini Batch K-Means and Gaussian Mixture Model clustering methods. Note tha GMM and Mini Batch K-Means are supported only for `scikit-learn` backend.
 
 Using K-Means (default):
 
@@ -124,7 +124,7 @@ cgram.plot()
 Using Mini Batch K-Means, which can provide significant speedup over K-Means:
 
 ```python
-cgram = Clustergram(range(1, 8), method='minibatchkmeans')
+cgram = Clustergram(range(1, 8), method='minibatchkmeans', batch_size=100)
 cgram.fit(data)
 cgram.plot()
 ```
@@ -153,6 +153,79 @@ cgram.plot(k_range=range(3, 10), figsize=(12, 8))
 ```
 ![Limited clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/limited_plot.png)
 
+## Additional clustering performance evaluation
+
+Clustergam includes handy wrappers around a selection of clustering performance metrics offered by
+`scikit-learn`. Data which were originally computed on GPU are converted to numpy on the fly.
+
+### Silhouette score
+
+Compute the mean Silhouette Coefficient of all samples. See [`scikit-learn` documentation](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html#sklearn.metrics.silhouette_score) for details.
+
+```python
+>>> cgram.silhouette_score()
+2    0.531540
+3    0.447219
+4    0.400154
+5    0.377720
+6    0.372128
+7    0.331575
+Name: silhouette_score, dtype: float64
+```
+Once computed, resulting Series is available as `cgram.silhouette`. Calling the original method will recompute the score.
+
+### Calinski and Harabasz score
+
+Compute the Calinski and Harabasz score, also known as the Variance Ratio Criterion. See [`scikit-learn` documentation](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.calinski_harabasz_score.html#sklearn.metrics.calinski_harabasz_score) for details.
+
+```python
+>>> cgram.calinski_harabasz_score()
+2    482.191469
+3    441.677075
+4    400.392131
+5    411.175066
+6    382.731416
+7    352.447569
+Name: calinski_harabasz_score, dtype: float64
+```
+Once computed, resulting Series is available as `cgram.calinski_harabasz`. Calling the original method will recompute the score.
+
+### Davies-Bouldin score
+
+Compute the Davies-Bouldin score. See [`scikit-learn` documentation](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.davies_bouldin_score.html#sklearn.metrics.davies_bouldin_score) for details.
+
+```python
+>>> cgram.davies_bouldin_score()
+2    0.714064
+3    0.943553
+4    0.943320
+5    0.973248
+6    0.950910
+7    1.074937
+Name: davies_bouldin_score, dtype: float64
+```
+Once computed, resulting Series is available as `cgram.davies_bouldin`. Calling the original method will recompute the score.
+
+## Acessing labels
+
+`Clustergram` stores resulting labels for each of the tested options, which can be accessed as:
+
+```python
+>>> cgram.labels
+     1  2  3  4  5  6  7
+0    0  0  2  2  3  2  1
+1    0  0  2  2  3  2  1
+2    0  0  2  2  3  2  1
+3    0  0  2  2  3  2  1
+4    0  0  2  2  0  0  3
+..  .. .. .. .. .. .. ..
+337  0  1  1  3  2  5  0
+338  0  1  1  3  2  5  0
+339  0  1  1  1  1  1  4
+340  0  1  1  3  2  5  5
+341  0  1  1  1  1  1  5
+```
+
 ## Saving clustergram
 
 You can save both plot and `clustergram.Clustergram` to a disk.


Test patch:
diff --git a/clustergram/test_clustergram.py b/clustergram/test_clustergram.py
index 3ce2ee5..7ca807b 100644
--- a/clustergram/test_clustergram.py
+++ b/clustergram/test_clustergram.py
@@ -1,10 +1,12 @@
 from sklearn.datasets import make_blobs
 import pandas as pd
+import numpy as np
 import pytest
 
 try:
     import cudf
     import cuml
+    import cupy as cp
 
     RAPIDS = True
 except (ImportError, ModuleNotFoundError):
@@ -12,7 +14,7 @@ except (ImportError, ModuleNotFoundError):
 
 from clustergram import Clustergram
 
-n_samples = 1000
+n_samples = 10
 n_features = 2
 
 n_clusters = 5
@@ -30,171 +32,368 @@ data = pd.DataFrame(device_data)
 
 
 def test_sklearn_kmeans():
-    clustergram = Clustergram(range(1, 10), backend="sklearn")
+    clustergram = Clustergram(range(1, 8), backend="sklearn", random_state=random_state)
     clustergram.fit(data)
 
-    pd.testing.assert_series_equal(
-        clustergram.means.mean(),
-        pd.Series(
-            [
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30046431,
-                -0.3003699,
-                -0.30086877,
-                -0.30016977,
-            ],
-            index=list(range(1, 10)),
-        ),
-        rtol=6,
+    for i in range(1, 8):
+        assert clustergram.labels[i].nunique() == i
+    assert clustergram.labels.shape == (10, 7)
+    assert clustergram.labels.notna().all().all()
+
+    expected = [
+        2.3448529748438847,
+        2.793993337861089,
+        2.8270525934022026,
+        2.513999519667852,
+        2.344852974843884,
+        2.059790669470653,
+        2.271614946625764,
+    ]
+    assert expected == [
+        pytest.approx(np.mean(clustergram.cluster_centers[x]), rel=1e-12)
+        for x in range(1, 8)
+    ]
+
+    assert clustergram.plot_data_pca.empty
+    ax = clustergram.plot(pca_kwargs=dict(random_state=random_state))
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data.empty
+    ax = clustergram.plot(pca_weighted=False)
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data_pca.mean().mean() == pytest.approx(
+        -0.21726383145643727, rel=1e-15
     )
+    assert clustergram.plot_data.mean().mean() == pytest.approx(
+        2.3448529748438847, rel=1e-15
+    )
+
 
-    ax = clustergram.plot()
+def test_sklearn_minibatchkmeans():
+    clustergram = Clustergram(
+        range(1, 8),
+        backend="sklearn",
+        method="minibatchkmeans",
+        random_state=random_state,
+    )
+    clustergram.fit(data)
+
+    for i in range(1, 8):
+        assert clustergram.labels[i].nunique() == i
+    assert clustergram.labels.shape == (10, 7)
+    assert clustergram.labels.notna().all().all()
+
+    expected = [
+        2.35104516103195,
+        2.7950777425566464,
+        2.8277733533088534,
+        2.522309579003138,
+        2.3452078554650075,
+        2.0602593594306526,
+        2.2717280475246584,
+    ]
+    assert expected == [
+        pytest.approx(np.mean(clustergram.cluster_centers[x]), rel=1e-12)
+        for x in range(1, 8)
+    ]
+
+    assert clustergram.plot_data_pca.empty
+    ax = clustergram.plot(pca_kwargs=dict(random_state=random_state))
     ax.get_geometry() == (1, 1, 1)
 
-    clustergram = Clustergram(range(1, 10), backend="sklearn", pca_weighted=False)
-    clustergram.fit(device_data)
+    assert clustergram.plot_data.empty
+    ax = clustergram.plot(pca_weighted=False)
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data_pca.mean().mean() == pytest.approx(
+        -0.19729006113986378, rel=1e-15
+    )
+    assert clustergram.plot_data.mean().mean() == pytest.approx(
+        2.3478208185941405, rel=1e-15
+    )
+
+
+def test_sklearn_gmm():
+    clustergram = Clustergram(
+        range(1, 8), backend="sklearn", method="gmm", random_state=random_state
+    )
+    clustergram.fit(data)
+
+    for i in range(1, 8):
+        assert clustergram.labels[i].nunique() == i
+    assert clustergram.labels.shape == (10, 7)
+    assert clustergram.labels.notna().all().all()
+
+    expected = [
+        2.6841643400794335,
+        3.1705471246923214,
+        2.52976769800815,
+        2.74946247582996,
+        2.333257525098582,
+        2.0444680314214723,
+        2.2697082687156915,
+    ]
+    assert expected == [
+        pytest.approx(np.mean(clustergram.cluster_centers[x]), rel=1e-12)
+        for x in range(1, 8)
+    ]
+
+    assert clustergram.plot_data_pca.empty
+    ax = clustergram.plot(pca_kwargs=dict(random_state=random_state))
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data.empty
+    ax = clustergram.plot(pca_weighted=False)
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data_pca.mean().mean() == pytest.approx(
+        -0.5099395641282745, rel=1e-15
+    )
+    assert clustergram.plot_data.mean().mean() == pytest.approx(
+        2.4439850629293924, rel=1e-15
+    )
+
+
+@pytest.mark.skipif(
+    not RAPIDS, reason="RAPIDS not available.",
+)
+def test_cuml_kmeans():
+    n_samples = 10
+    n_features = 2
+
+    n_clusters = 5
+    random_state = 0
+
+    device_data, device_labels = cuml.make_blobs(
+        n_samples=n_samples,
+        n_features=n_features,
+        centers=n_clusters,
+        random_state=random_state,
+        cluster_std=0.1,
+    )
+
+    data = cudf.DataFrame(device_data)
+
+    # cudf.DataFrame
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
+    clustergram.fit(data)
+
+    for i in range(1, 8):
+        assert clustergram.labels[i].nunique() == i
+    assert clustergram.labels.shape == (10, 7)
+    assert clustergram.labels.notna().all().all()
+
+    expected = [
+        3.7674055099487305,
+        2.7064273357391357,
+        3.451129913330078,
+        4.223802089691162,
+        4.125243663787842,
+        2.953890800476074,
+        3.4818685054779053,
+    ]
+    assert expected == [
+        pytest.approx(float(clustergram.cluster_centers[x].mean().mean()), rel=1e-6)
+        for x in range(1, 8)
+    ]
+
+    assert clustergram.plot_data_pca.empty
+    ax = clustergram.plot(pca_kwargs=dict(random_state=random_state))
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data.empty
+    ax = clustergram.plot(pca_weighted=False)
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data_pca.mean().mean() == pytest.approx(
+        1.1016593081610544, rel=1e-15
+    )
+    assert clustergram.plot_data.mean().mean() == pytest.approx(
+        3.7674053737095425, rel=1e-15
+    )
+
+    # cupy array
+    data = device_data
+
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
+    clustergram.fit(data)
+
+    for i in range(1, 8):
+        assert clustergram.labels[i].nunique() == i
+    assert clustergram.labels.shape == (10, 7)
+    assert clustergram.labels.notna().all().all()
+
+    expected = [
+        3.7674055099487305,
+        2.7064273357391357,
+        3.451129913330078,
+        4.223802089691162,
+        4.125243663787842,
+        2.953890800476074,
+        3.4818685054779053,
+    ]
+    assert expected == [
+        pytest.approx(float(cp.mean(clustergram.cluster_centers[x])), rel=1e-6)
+        for x in range(1, 8)
+    ]
+
+    assert clustergram.plot_data_pca.empty
+    ax = clustergram.plot(pca_kwargs=dict(random_state=random_state))
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data.empty
+    ax = clustergram.plot(pca_weighted=False)
+    ax.get_geometry() == (1, 1, 1)
+
+    assert clustergram.plot_data_pca.mean().mean() == pytest.approx(
+        1.1016593081610544, rel=1e-15
+    )
+    assert clustergram.plot_data.mean().mean() == pytest.approx(
+        3.7674053737095425, rel=1e-15
+    )
+
+
+def test_errors():
+    with pytest.raises(ValueError):
+        Clustergram(range(1, 3), backend="nonsense")
+    with pytest.raises(ValueError):
+        Clustergram(range(1, 3), method="nonsense")
+
+
+def test_repr_():
+    expected = (
+        "Clustergram(k_range=range(1, 30), backend='sklearn', "
+        "method='kmeans', kwargs={'n_init': 10})"
+    )
+    clustergram = Clustergram(range(1, 30), n_init=10)
+    assert expected == clustergram.__repr__()
+
+
+def test_silhouette_score():
+    clustergram = Clustergram(range(1, 8), backend="sklearn", random_state=random_state)
+    clustergram.fit(data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean(),
+        clustergram.silhouette_score(),
         pd.Series(
-            [
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05522513,
-                2.05582383,
-                2.05585776,
-                2.05591782,
-            ],
-            index=list(range(1, 10)),
+            [0.70244987, 0.64427202, 0.76772759, 0.94899084, 0.76998519, 0.57564372],
+            index=list(range(2, 8)),
+            name="silhouette_score",
         ),
-        rtol=6,
     )
 
-    ax = clustergram.plot()
-    ax.get_geometry() == (1, 1, 1)
+    pd.testing.assert_series_equal(
+        clustergram.silhouette,
+        pd.Series(
+            [0.70244987, 0.64427202, 0.76772759, 0.94899084, 0.76998519, 0.57564372],
+            index=list(range(2, 8)),
+            name="silhouette_score",
+        ),
+    )
 
 
-def test_sklearn_minibatchkmeans():
-    clustergram = Clustergram(range(1, 10), backend="sklearn", method="minibatchkmeans")
+@pytest.mark.skipif(
+    not RAPIDS, reason="RAPIDS not available.",
+)
+def test_silhouette_score_cuml():
+    n_samples = 10
+    n_features = 2
+
+    n_clusters = 5
+    random_state = 0
+
+    device_data, device_labels = cuml.make_blobs(
+        n_samples=n_samples,
+        n_features=n_features,
+        centers=n_clusters,
+        random_state=random_state,
+        cluster_std=0.1,
+    )
+
+    data = cudf.DataFrame(device_data)
+
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
     clustergram.fit(data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean(),
+        clustergram.silhouette_score(),
         pd.Series(
             [
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30046431,
-                -0.3003699,
-                -0.30086877,
-                -0.30016977,
+                0.7494349479675293,
+                0.9806153178215027,
+                0.6721830368041992,
+                0.39418715238571167,
+                0.44574037194252014,
+                0.08033210784196854,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="silhouette_score",
         ),
-        rtol=6,
     )
 
-    ax = clustergram.plot()
-    ax.get_geometry() == (1, 1, 1)
-
-    clustergram = Clustergram(
-        range(1, 10), backend="sklearn", method="minibatchkmeans", pca_weighted=False
-    )
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
     clustergram.fit(device_data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean(),
+        clustergram.silhouette_score(),
         pd.Series(
             [
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05522513,
-                2.05582383,
-                2.05585776,
-                2.05591782,
+                0.7494349479675293,
+                0.9806153178215027,
+                0.6721830368041992,
+                0.39418715238571167,
+                0.44574037194252014,
+                0.08033210784196854,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="silhouette_score",
         ),
-        rtol=6,
     )
 
-    ax = clustergram.plot()
-    ax.get_geometry() == (1, 1, 1)
 
-
-def test_sklearn_gmm():
-    clustergram = Clustergram(range(1, 10), backend="sklearn", method="gmm")
+def test_calinski_harabasz_score():
+    clustergram = Clustergram(range(1, 8), backend="sklearn", random_state=random_state)
     clustergram.fit(data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean(),
+        clustergram.calinski_harabasz_score(),
         pd.Series(
             [
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30042205,
-                -0.30046431,
-                -0.3003699,
-                -0.30086877,
-                -0.30016977,
+                23.17662874,
+                30.64301789,
+                55.22333618,
+                3116.43518408,
+                3899.06868932,
+                4439.30604863,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="calinski_harabasz_score",
         ),
-        rtol=6,
     )
 
-    ax = clustergram.plot()
-    ax.get_geometry() == (1, 1, 1)
-
-    clustergram = Clustergram(
-        range(1, 10), backend="sklearn", method="gmm", pca_weighted=False
-    )
-    clustergram.fit(device_data)
-
     pd.testing.assert_series_equal(
-        clustergram.means.mean(),
+        clustergram.calinski_harabasz,
         pd.Series(
             [
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05533664,
-                2.05522513,
-                2.05582383,
-                2.05585776,
-                2.05591782,
+                23.17662874,
+                30.64301789,
+                55.22333618,
+                3116.43518408,
+                3899.06868932,
+                4439.30604863,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="calinski_harabasz_score",
         ),
-        rtol=6,
     )
 
-    ax = clustergram.plot()
-    ax.get_geometry() == (1, 1, 1)
-
 
 @pytest.mark.skipif(
-    not RAPIDS,
-    reason="RAPIDS not available.",
+    not RAPIDS, reason="RAPIDS not available.",
 )
-def test_cuml_kmeans():
-    n_samples = 1000
+def test_calinski_harabasz_score_cuml():
+    n_samples = 10
     n_features = 2
 
     n_clusters = 5
@@ -210,53 +409,79 @@ def test_cuml_kmeans():
 
     data = cudf.DataFrame(device_data)
 
-    clustergram = Clustergram(range(1, 10), backend="cuML")
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
     clustergram.fit(data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean().to_pandas(),
+        clustergram.calinski_harabasz_score(),
         pd.Series(
             [
-                0.30042205,
-                0.30042205,
-                0.30042205,
-                0.30042205,
-                0.30042205,
-                0.30046431,
-                0.3003699,
-                0.30086877,
-                0.30016977,
+                25.619150510634366,
+                15374.042816067375,
+                10813.16845006968,
+                8818.1163716754,
+                8070.657293970755,
+                7259.89764652579,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="calinski_harabasz_score",
         ),
-        rtol=6,
     )
 
-    data = device_data
-
-    clustergram = Clustergram(range(1, 10), backend="cuML")
-    clustergram.fit(data)
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
+    clustergram.fit(device_data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean().to_pandas(),
+        clustergram.calinski_harabasz_score(),
         pd.Series(
             [
-                0.30042205,
-                0.30042205,
-                0.30042205,
-                0.30042205,
-                0.30042205,
-                0.30046431,
-                0.3003699,
-                0.30086877,
-                0.30016977,
+                25.619150510634366,
+                15374.042816067375,
+                10813.16845006968,
+                8818.1163716754,
+                8070.657293970755,
+                7259.89764652579,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="calinski_harabasz_score",
+        ),
+    )
+
+
+def test_davies_bouldin_score():
+    clustergram = Clustergram(range(1, 8), backend="sklearn", random_state=random_state)
+    clustergram.fit(data)
+
+    pd.testing.assert_series_equal(
+        clustergram.davies_bouldin_score(),
+        pd.Series(
+            [0.2493657, 0.35181197, 0.34758021, 0.05567944, 0.03051626, 0.02520726],
+            index=list(range(2, 8)),
+            name="davies_bouldin_score",
+        ),
+    )
+
+    pd.testing.assert_series_equal(
+        clustergram.davies_bouldin,
+        pd.Series(
+            [0.2493657, 0.35181197, 0.34758021, 0.05567944, 0.03051626, 0.02520726],
+            index=list(range(2, 8)),
+            name="davies_bouldin_score",
         ),
-        rtol=6,
     )
 
-    device_data, device_labels = make_blobs(
+
+@pytest.mark.skipif(
+    not RAPIDS, reason="RAPIDS not available.",
+)
+def test_davies_bouldin_score_cuml():
+    n_samples = 10
+    n_features = 2
+
+    n_clusters = 5
+    random_state = 0
+
+    device_data, device_labels = cuml.make_blobs(
         n_samples=n_samples,
         n_features=n_features,
         centers=n_clusters,
@@ -264,65 +489,42 @@ def test_cuml_kmeans():
         cluster_std=0.1,
     )
 
-    clustergram = Clustergram(range(1, 10), backend="cuML", pca_weighted=False)
+    data = cudf.DataFrame(device_data)
+
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
     clustergram.fit(data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean().to_pandas(),
+        clustergram.davies_bouldin_score(),
         pd.Series(
             [
-                2.31314663,
-                2.31314663,
-                2.31314663,
-                2.31314663,
-                2.31314663,
-                2.31353356,
-                2.31314663,
-                2.31326447,
-                2.31284071,
+                0.3107512701086121,
+                0.02263161666570639,
+                0.2261582258142144,
+                0.3839688146565784,
+                0.13388392354928222,
+                0.279734367840293,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="davies_bouldin_score",
         ),
-        rtol=6,
     )
 
-    clustergram = Clustergram(range(1, 10), backend="cuML", pca_weighted=False)
-    clustergram.fit(cudf.DataFrame(device_data))
+    clustergram = Clustergram(range(1, 8), backend="cuML", random_state=random_state)
+    clustergram.fit(device_data)
 
     pd.testing.assert_series_equal(
-        clustergram.means.mean().to_pandas(),
+        clustergram.davies_bouldin_score(),
         pd.Series(
             [
-                2.31314663,
-                2.31314663,
-                2.31314663,
-                2.31314663,
-                2.31314663,
-                2.31353356,
-                2.31314663,
-                2.31326447,
-                2.31284071,
+                0.3107512701086121,
+                0.02263161666570639,
+                0.2261582258142144,
+                0.3839688146565784,
+                0.13388392354928222,
+                0.279734367840293,
             ],
-            index=list(range(1, 10)),
+            index=list(range(2, 8)),
+            name="davies_bouldin_score",
         ),
-        rtol=6,
-    )
-
-    ax = clustergram.plot()
-    ax.get_geometry() == (1, 1, 1)
-
-
-def test_errors():
-    with pytest.raises(ValueError):
-        Clustergram(range(1, 3), backend="nonsense")
-    with pytest.raises(ValueError):
-        Clustergram(range(1, 3), method="nonsense")
-
-
-def test_repr_():
-    expected = (
-        "Clustergram(k_range=range(1, 30), backend='sklearn', "
-        "method='kmeans', pca_weighted=True, kwargs={'n_init': 10})"
     )
-    clustergram = Clustergram(range(1, 30), n_init=10)
-    assert expected == clustergram.__repr__()