Repository: martinfleis/clustergram

Issue:
Support hierarchical clustering
Supporting hierarchical clustering as in the original Schonlau's paper would be nice. Not sure if using scipy or sklearn, will have to explore.

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=h1) Report
> Merging [#11](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=desc) (126edde) into [main](https://codecov.io/gh/martinfleis/clustergram/commit/5bf2cdcbd53b35cac8d43bbe6d15c4157311a0b0?el=desc) (5bf2cdc) will **increase** coverage by `1.70%`.
> The diff coverage is `89.74%`.

[![Impacted file tree graph](https://codecov.io/gh/martinfleis/clustergram/pull/11/graphs/tree.svg?width=650&height=150&src=pr&token=4IL16GQR4P)](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##             main      #11      +/-   ##
==========================================
+ Coverage   55.30%   57.00%   +1.70%     
==========================================
  Files           3        3              
  Lines         481      507      +26     
==========================================
+ Hits          266      289      +23     
- Misses        215      218       +3     
```


| [Impacted Files](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [clustergram/clustergram.py](https://codecov.io/gh/martinfleis/clustergram/pull/11/diff?src=pr&el=tree#diff-Y2x1c3RlcmdyYW0vY2x1c3RlcmdyYW0ucHk=) | `75.22% <89.74%> (+1.75%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=footer). Last update [5bf2cdc...126edde](https://codecov.io/gh/martinfleis/clustergram/pull/11?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/README.md b/README.md
index 2e70303..de31fbe 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, Mini Batch 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, plus hierarchical/agglomerative clustering using `SciPy`.
 
 ## Getting started
 
@@ -25,7 +25,7 @@ pip install clustergram
 ```
 
 In any case, you still need to install your selected backend
-(`scikit-learn` or `cuML`).
+(`scikit-learn` and `scipy` or `cuML`).
 
 The example of clustergram on Palmer penguins dataset:
 
@@ -89,9 +89,9 @@ cgram.plot(figsize=(12, 8), pca_weighted=False)
 ```
 ![Default clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/pca_false.png)
 
-## Scikit-learn and RAPIDS cuML backends
+## Scikit-learn, SciPy and RAPIDS cuML backends
 
-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.
+Clustergram offers three backends for the computation - `scikit-learn` and `scipy` which use CPU and RAPIDS.AI `cuML`, which uses GPU. Note that all are optional dependencies but you will need at least one of them to generate clustergram.
 
 Using `scikit-learn` (default):
 
@@ -113,7 +113,7 @@ cgram.plot()
 
 ## Supported methods
 
-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.
+Clustergram currently supports K-Means, Mini Batch K-Means, Gaussian Mixture Model and SciPy's hierarchical clustering methods. Note tha GMM and Mini Batch K-Means are supported only for `scikit-learn` backend and hierarchical methods are supported only for `scipy` backend.
 
 Using K-Means (default):
 
@@ -139,6 +139,14 @@ cgram.fit(data)
 cgram.plot()
 ```
 
+Using Ward's hierarchical clustering:
+
+```python
+cgram = Clustergram(range(1, 8), method='hierarchical', linkage='ward')
+cgram.fit(data)
+cgram.plot()
+```
+
 ## Partial plot
 
 `Clustergram.plot()` can also plot only a part of the diagram, if you want to focus on a limited range of `k`.
diff --git a/clustergram/clustergram.py b/clustergram/clustergram.py
index 8f6f362..d849574 100644
--- a/clustergram/clustergram.py
+++ b/clustergram/clustergram.py
@@ -25,20 +25,26 @@ class Clustergram:
     as k-means and for hierarchical cluster algorithms when the number of
     observations is large enough to make dendrograms impractical.
 
-    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
+    Clustergram offers three backends for the computation - ``scikit-learn`` and
+    ``scipy`` which use CPU and RAPIDS.AI ``cuML``, which uses GPU. Note that all
+    are optional dependencies but you will need at least one of them to
     generate clustergram.
 
     Parameters
     ----------
-    k_range : iterable
+    k_range : iterable (default None)
         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.
-    method : {'kmeans', 'gmm', 'minibatchkmeans'} (default 'kmeans')
+        components). Not required for hierarchical clustering but will be applied
+        if given. It is recommended to always use limited range for hierarchical
+        methods as unlimited clustergram can take a while to compute and for large
+        number of observations is not legible.
+    backend : {'sklearn', 'cuML', 'scipy'} (default None)
+        Specify computational backend. Defaults to ``sklearn`` for ``'kmeans'``,
+        ``'gmm'``, and ``'minibatchkmeans'`` methods and to ``'scipy'`` for any of
+        hierarchical clustering methods. ``'scipy'`` uses ``sklearn`` for PCA
+        computation if that is required.
+        ``sklearn`` does computation on CPU, ``cuml`` on GPU.
+    method : {'kmeans', 'gmm', 'minibatchkmeans', 'hierarchical'} (default 'kmeans')
         Clustering method.
 
         * ``kmeans`` uses K-Means clustering, either as ``sklearn.cluster.KMeans``
@@ -46,6 +52,8 @@ class Clustergram:
         * ``gmm`` uses Gaussian Mixture Model as ``sklearn.mixture.GaussianMixture``
         * ``minibatchkmeans`` uses Mini Batch K-Means as
           ``sklearn.cluster.MiniBatchKMeans``
+        * ``hierarchical`` uses hierarchical/agglomerative clustering as
+          ``scipy.cluster.hierarchy.linkage``. See
 
         Note that ``gmm`` and ``minibatchkmeans`` are currently supported only
         with ``sklearn`` backend.
@@ -53,7 +61,10 @@ class Clustergram:
         Print progress and time of individual steps.
     **kwargs
         Additional arguments passed to the model (e.g. ``KMeans``),
-        e.g. ``random_state``.
+        e.g. ``random_state``. Pass ``linkage`` to specify linkage method in
+        case of hierarchical clustering (e.g. ``linkage='ward'``). See the
+        documentation of scipy for details:
+        https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html#scipy.cluster.hierarchy.linkage
 
 
     Attributes
@@ -63,6 +74,8 @@ class Clustergram:
         DataFrame with cluster labels for each iteration.
     cluster_centers : dict
         Dictionary with cluster centers for each iteration.
+    linkage : scipy.cluster.hierarchy.linkage
+        Linkage object for hierarchical methods.
 
 
     Examples
@@ -91,7 +104,12 @@ class Clustergram:
     """
 
     def __init__(
-        self, k_range, backend="sklearn", method="kmeans", verbose=True, **kwargs,
+        self,
+        k_range=None,
+        backend=None,
+        method="kmeans",
+        verbose=True,
+        **kwargs,
     ):
         self.k_range = k_range
 
@@ -99,14 +117,19 @@ class Clustergram:
         kwargs.pop("pca_weighted", None)
         kwargs.pop("pca_kwargs", None)
 
-        if backend not in ["sklearn", "cuML"]:
+        if backend is None:
+            backend = "scipy" if method == "hierarchical" else "sklearn"
+
+        allowed_backends = ["sklearn", "cuML", "scipy"]
+        if backend not in allowed_backends:
             raise ValueError(
-                f'"{backend}" is not a supported backend. Use "sklearn" or "cuML".'
+                f'"{backend}" is not a supported backend. '
+                f"Use one of {allowed_backends}."
             )
         else:
             self.backend = backend
 
-        supported = ["kmeans", "gmm", "minibatchkmeans"]
+        supported = ["kmeans", "gmm", "minibatchkmeans", "hierarchical"]
         if method not in supported:
             raise ValueError(
                 f"'{method}' is not a supported method. "
@@ -115,10 +138,23 @@ class Clustergram:
         else:
             self.method = method
 
+        if self.k_range is None and self.method != "hierarchical":
+            raise ValueError(f"'k_range' is mandatory for '{self.method}' method.")
+
+        if (
+            (self.backend == "cuML" and self.method != "kmeans")
+            or (self.backend == "scipy" and self.method != "hierarchical")
+            or (self.backend == "sklearn" and self.method == "hierarchical")
+        ):
+            raise ValueError(
+                f"'{self.method}' method is not implemented "
+                f"for '{self.backend}' backend. Use supported combination."
+            )
+
         self.engine_kwargs = kwargs
         self.verbose = verbose
 
-        if self.backend == "sklearn":
+        if self.backend in ["sklearn", "scipy"]:
             self.plot_data_pca = pd.DataFrame()
             self.plot_data = pd.DataFrame()
         else:
@@ -164,9 +200,11 @@ class Clustergram:
             elif self.method == "minibatchkmeans":
                 self._kmeans_sklearn(data, minibatch=True, **kwargs)
             elif self.method == "gmm":
-                self.means = self._gmm_sklearn(data, **kwargs)
+                self._gmm_sklearn(data, **kwargs)
         if self.backend == "cuML":
-            self.means = self._kmeans_cuml(data, **kwargs)
+            self._kmeans_cuml(data, **kwargs)
+        if self.backend == "scipy":
+            self._scipy_hierarchical(data, **kwargs)
 
     def _kmeans_sklearn(self, data, minibatch, **kwargs):
         """Use scikit-learn KMeans"""
@@ -217,7 +255,6 @@ class Clustergram:
         """Use sklearn.mixture.GaussianMixture"""
         try:
             from sklearn.mixture import GaussianMixture
-            import numpy as np
             from scipy.stats import multivariate_normal
         except ImportError:
             raise ImportError(
@@ -250,6 +287,33 @@ class Clustergram:
 
             print(f"K={n} fitted in {time() - s} seconds.") if self.verbose else None
 
+    def _scipy_hierarchical(self, data, **kwargs):
+        """Use scipy.cluster.hierarchy.linkage"""
+        try:
+            from scipy.cluster import hierarchy
+        except ImportError:
+            raise ImportError("scipy is required to use `scipy` backend.")
+
+        method = self.engine_kwargs.pop("linkage", "single")
+        self.linkage = hierarchy.linkage(data, method=method, **self.engine_kwargs)
+        rootnode, nodelist = hierarchy.to_tree(self.linkage, rd=True)
+        distances = [node.dist for node in nodelist if node.dist > 0][::-1]
+
+        self.labels = pd.DataFrame()
+        self.cluster_centers = {}
+
+        if self.k_range is None:
+            self.k_range = range(1, len(distances) + 1)
+
+        if not isinstance(data, pd.DataFrame):
+            data = pd.DataFrame(data)
+
+        for i in self.k_range:
+            d = distances[i - 1]
+            lab = hierarchy.fcluster(self.linkage, d, criterion="distance")
+            self.labels[i] = lab - 1
+            self.cluster_centers[i] = data.groupby(lab).mean().values
+
     def silhouette_score(self, **kwargs):
         """
         Compute the mean Silhouette Coefficient of all samples.
@@ -304,7 +368,7 @@ class Clustergram:
 
         self.silhouette = pd.Series(name="silhouette_score", dtype="float64")
 
-        if self.backend == "sklearn":
+        if self.backend in ["sklearn", "scipy"]:
             for k in self.k_range:
                 if k > 1:
                     self.silhouette.loc[k] = metrics.silhouette_score(
@@ -375,7 +439,7 @@ class Clustergram:
             name="calinski_harabasz_score", dtype="float64"
         )
 
-        if self.backend == "sklearn":
+        if self.backend in ["sklearn", "scipy"]:
             for k in self.k_range:
                 if k > 1:
                     self.calinski_harabasz.loc[k] = metrics.calinski_harabasz_score(
@@ -441,7 +505,7 @@ class Clustergram:
 
         self.davies_bouldin = pd.Series(name="davies_bouldin_score", dtype="float64")
 
-        if self.backend == "sklearn":
+        if self.backend in ["sklearn", "scipy"]:
             for k in self.k_range:
                 if k > 1:
                     self.davies_bouldin.loc[k] = metrics.davies_bouldin_score(
@@ -569,13 +633,13 @@ class Clustergram:
             if self.plot_data_pca.empty:
                 pca_kwargs.pop("n_components", None)
 
-                if self.backend == "sklearn":
+                if self.backend in ["sklearn", "scipy"]:
                     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":
+                if self.backend in ["sklearn", "scipy"]:
                     self._compute_means_sklearn()
                 else:
                     self._compute_means_cuml()
@@ -612,7 +676,7 @@ class Clustergram:
         for i in k_range:
             cl = means[i].value_counts()
 
-            if self.backend == "sklearn":
+            if self.backend in ["sklearn", "scipy"]:
                 ax.scatter(
                     [i] * i,
                     [cl.index],
@@ -636,7 +700,7 @@ class Clustergram:
                 )
 
             try:
-                if self.backend == "sklearn":
+                if self.backend in ["sklearn", "scipy"]:
                     sub = means.groupby([i, i + 1]).count().reset_index()
                 else:
                     sub = means.groupby([i, i + 1]).count().reset_index().to_pandas()
diff --git a/doc/index.md b/doc/index.md
index 7d6fd98..77df37c 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, Mini Batch 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, plus hierarchical/agglomerative clustering using `SciPy`.
 
 ## Getting started
 
@@ -23,7 +23,7 @@ pip install clustergram
 ```
 
 In any case, you still need to install your selected backend
-(`scikit-learn` or `cuML`).
+(`scikit-learn` and `scipy` or `cuML`).
 
 The example of clustergram on Palmer penguins dataset:
 
@@ -87,9 +87,9 @@ cgram.plot(figsize=(12, 8), pca_weighted=False)
 ```
 ![Default clustergram](https://raw.githubusercontent.com/martinfleis/clustergram/master/doc/_static/pca_false.png)
 
-## Scikit-learn and RAPIDS cuML backends
+## Scikit-learn, SciPy and RAPIDS cuML backends
 
-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.
+Clustergram offers three backends for the computation - `scikit-learn` and `scipy` which use CPU and RAPIDS.AI `cuML`, which uses GPU. Note that all are optional dependencies but you will need at least one of them to generate clustergram.
 
 Using `scikit-learn` (default):
 
@@ -111,7 +111,7 @@ cgram.plot()
 
 ## Supported methods
 
-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.
+Clustergram currently supports K-Means, Mini Batch K-Means, Gaussian Mixture Model and SciPy's hierarchical clustering methods. Note tha GMM and Mini Batch K-Means are supported only for `scikit-learn` backend and hierarchical methods are supported only for `scipy` backend.
 
 Using K-Means (default):
 
@@ -137,6 +137,14 @@ cgram.fit(data)
 cgram.plot()
 ```
 
+Using Ward's hierarchical clustering:
+
+```python
+cgram = Clustergram(range(1, 8), method='hierarchical', linkage='ward')
+cgram.fit(data)
+cgram.plot()
+```
+
 ## Partial plot
 
 `Clustergram.plot()` can also plot only a part of the diagram, if you want to focus on a limited range of `k`.


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index a1095a2..09072fe 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -2,9 +2,9 @@ name: Tests
 
 on:
   push:
-    branches: [ master ]
+    branches: [ main ]
   pull_request:
-    branches: [ master ]
+    branches: [ main ]
   schedule:
     - cron:  '0 0 * * 5'
 
@@ -12,7 +12,7 @@ jobs:
   Test:
     name: Test clustergram
     runs-on: [ubuntu-latest]
-    
+
     steps:
       - uses: actions/checkout@v2
       - name: Install micromamba
@@ -28,7 +28,7 @@ jobs:
           export MAMBA_EXE=$(pwd)/micromamba
           . $MAMBA_ROOT_PREFIX/etc/profile.d/mamba.sh
           ./micromamba create -f ci/envs/latest.yaml -y
-      
+
       - name: Test clustergram
         shell: bash -l {0}
         run: |
diff --git a/clustergram/test_clustergram.py b/clustergram/test_clustergram.py
index 7ca807b..c391e72 100644
--- a/clustergram/test_clustergram.py
+++ b/clustergram/test_clustergram.py
@@ -156,7 +156,8 @@ def test_sklearn_gmm():
 
 
 @pytest.mark.skipif(
-    not RAPIDS, reason="RAPIDS not available.",
+    not RAPIDS,
+    reason="RAPIDS not available.",
 )
 def test_cuml_kmeans():
     n_samples = 10
@@ -254,14 +255,71 @@ def test_cuml_kmeans():
     )
 
 
+def test_hierarchical():
+    clustergram = Clustergram(range(1, 8), method="hierarchical")
+    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.344852974843885,
+        2.793993337861089,
+        2.827052593402202,
+        2.513999519667852,
+        2.3448529748438847,
+        2.059790669470653,
+        2.2716149466257645,
+    ]
+    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.2172638314564372, rel=1e-15
+    )
+    assert clustergram.plot_data.mean().mean() == pytest.approx(
+        2.3448529748438847, rel=1e-15
+    )
+
+
+def test_hierarchical_array():
+    clustergram = Clustergram(method="hierarchical")
+    clustergram.fit(data.values)
+
+    for i in range(1, 10):
+        assert clustergram.labels[i].nunique() == i
+    assert clustergram.labels.shape == (10, 9)
+    assert clustergram.labels.notna().all().all()
+
+
 def test_errors():
     with pytest.raises(ValueError):
         Clustergram(range(1, 3), backend="nonsense")
     with pytest.raises(ValueError):
         Clustergram(range(1, 3), method="nonsense")
+    with pytest.raises(ValueError):
+        Clustergram(range(1, 3), method="kmeans", backend="scipy")
+    with pytest.raises(ValueError):
+        Clustergram(range(1, 3), method="hieararchical", backend="sklearn")
+    with pytest.raises(ValueError):
+        Clustergram(range(1, 3), method="gmm", backend="cuML")
+    with pytest.raises(ValueError):
+        Clustergram()
 
 
-def test_repr_():
+def test_repr():
     expected = (
         "Clustergram(k_range=range(1, 30), backend='sklearn', "
         "method='kmeans', kwargs={'n_init': 10})"
@@ -294,7 +352,8 @@ def test_silhouette_score():
 
 
 @pytest.mark.skipif(
-    not RAPIDS, reason="RAPIDS not available.",
+    not RAPIDS,
+    reason="RAPIDS not available.",
 )
 def test_silhouette_score_cuml():
     n_samples = 10
@@ -390,7 +449,8 @@ def test_calinski_harabasz_score():
 
 
 @pytest.mark.skipif(
-    not RAPIDS, reason="RAPIDS not available.",
+    not RAPIDS,
+    reason="RAPIDS not available.",
 )
 def test_calinski_harabasz_score_cuml():
     n_samples = 10
@@ -472,7 +532,8 @@ def test_davies_bouldin_score():
 
 
 @pytest.mark.skipif(
-    not RAPIDS, reason="RAPIDS not available.",
+    not RAPIDS,
+    reason="RAPIDS not available.",
 )
 def test_davies_bouldin_score_cuml():
     n_samples = 10