Repository: scikit-learn/scikit-learn
Pull request: Added distance_threshold parameter to hierarchical clustering
#### Reference Issue 
Fixes #3796 

#### What does this implement/fix? Explain your changes.

Hierarchical clustering now has a `distance_threshold` parameter which can be set instead of `n_clusters` and this is used to determine the number of clusters to cut the tree at.

*Note: that this works only when compute_full_tree=True and n_clusters=None.*

Changes
* Either `distance_threshold` or `n_clusters` is accepted as parameter. When `distance_threshold` is set  `n_clusters` needs to be set to None as the default of 2 clusters has been retained.
* When building the tree the return_distance is set to True if the `distance_threshold` has been set.
The distances returned is then used to calculate the number of clusters when cutting the tree.

Tests
* Updated distance_threshold tests to pass with `n_clusters=None` as parameter as the default for n_clusters is 2.
* Test to raise error when neither n_cluster nor distance_threshold is passed in.
* Test agglomerative clustering with distance_threshold passed in and compare with the different number of clusters produced with and without connectivity.

Documentation
* Updated the User Guide to show `distance_threshold` parameter in the table.
* Doc string updated with `distance_threshold` explanation

#### Any other comments?

Is this implementation the right behaviour of `distance_threshold` or should the construction of the tree stop when the distance is reached?

--------
Apologies for the constant stream of commits, I was trying to sort out the errors myself.

Problem statement:
Hierarchical clustering: distance threshold
So, the output of hierarchichal clustering results can be determined either by number of clusters, or by the a distance thereshold to cut the tree at that threshold. However, scikit learn only supports one way! 

class sklearn.cluster.AgglomerativeClustering(n_clusters=2, ...

As a suggestion, is it possible to add the other option, to give distance_threshold as input argument and get as many clusters are created in the output! So, one of this argument should be given, either n_clusters or the threshold, not both!

Hints from issue discussion:
sure it is possible !

code is open and we always welcome good contributions...

actually there are more options for the stopping/threshold criteria:
https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.cluster.hierarchy.fcluster.html#scipy.cluster.hierarchy.fcluster
http://www.mathworks.com/help/stats/cluster.html

the inconsistency coefficient seems interesting and is the default in matlab and scipy...

You might also review some of the literate on graph modularity as a more robust measure of cluster quality.  The igraph package (with python bindings) has a large number of these functions available for reference.

Some discussion here also re: clustering and "communitiy detection", agglomerative/divisive/spectrall approaches, etc.:
http://arxiv.org/abs/0906.0612
http://bommaritollc.com/2012/06/summary-community-detection-algorithms-igraph-0-6/

Do you mind if I pick this up? I have implemented clustering in C++ years ago so I'll brush up and get back with questions.
I think you'd be welcome to

On 11 May 2017 at 20:12, Vathsala Achar <notifications@github.com> wrote:

> Do you mind if I pick this up? I have implemented clustering in C++ years
> ago so I'll brush up and get back with questions.
>
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/scikit-learn/scikit-learn/issues/3796#issuecomment-300746178>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AAEz66oNne-dsdHaN8gyrMBFploImN5bks5r4t8ngaJpZM4CyRYi>
> .
>

Thanks @jnothman 
I saw that @vmichel has made the necessary changes for a distance threshold (and no PR?), but I was interested in the `criterion` that `scipy.fcluster` uses (mention by @eyaler) and was wondering if that would be useful along with the distance threshold. Or is it a feature of it's own?

designing an api that is extensible to alternate stopping criteria may bea
good idea, but it's often good to implement one at a time

On 12 May 2017 7:45 pm, "Vathsala Achar" <notifications@github.com> wrote:

> Thanks @jnothman <https://github.com/jnothman>
> I saw that @vmichel <https://github.com/vmichel> has made the necessary
> changes for a distance threshold (and no PR?), but I was interested in the
> criterion that scipy.fcluster uses (mention by @eyaler
> <https://github.com/eyaler>) and was wondering if that would be useful
> along with the distance threshold. Or is it a feature of it's own?
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/scikit-learn/scikit-learn/issues/3796#issuecomment-301033531>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AAEz668aVTFBtfzK_lPIIEOJUhCcNc0Vks5r5CpNgaJpZM4CyRYi>
> .
>

Commit message: Changes based on review

* Doc string updates for clear information

* Removed redundant attribute n_clusters_

* Fixed tests

* Changes to FeatureAgglomeration to include distance threshold

Patch to review:
diff --git a/doc/modules/clustering.rst b/doc/modules/clustering.rst
index 9dfb0d08eaa41..388bfb8200158 100644
--- a/doc/modules/clustering.rst
+++ b/doc/modules/clustering.rst
@@ -73,13 +73,13 @@ Overview of clustering methods
      - Graph distance (e.g. nearest-neighbor graph)
 
    * - :ref:`Ward hierarchical clustering <hierarchical_clustering>`
-     - number of clusters
+     - number of clusters or distance threshold
      - Large ``n_samples`` and ``n_clusters``
      - Many clusters, possibly connectivity constraints
      - Distances between points
 
    * - :ref:`Agglomerative clustering <hierarchical_clustering>`
-     - number of clusters, linkage type, distance
+     - number of clusters or distance threshold, linkage type, distance
      - Large ``n_samples`` and ``n_clusters``
      - Many clusters, possibly connectivity constraints, non Euclidean
        distances
diff --git a/sklearn/cluster/hierarchical.py b/sklearn/cluster/hierarchical.py
index deb0bb5b0c23c..55d20cc5b5f3a 100644
--- a/sklearn/cluster/hierarchical.py
+++ b/sklearn/cluster/hierarchical.py
@@ -541,7 +541,8 @@ def _average_linkage(*args, **kwargs):
 ###############################################################################
 # Functions for cutting  hierarchical clustering tree
 
-def _hc_cut(n_clusters, children, n_leaves):
+def _hc_cut(n_clusters, children, n_leaves,
+            distance_threshold=None, distances=None):
     """Function cutting the ward tree for a given number of clusters.
 
     Parameters
@@ -560,12 +561,29 @@ def _hc_cut(n_clusters, children, n_leaves):
     n_leaves : int
         Number of leaves of the tree.
 
+    distance_threshold : float (optional)
+        The distance threshold to cluster at.
+        If distance_threshold is set then n_clusters will be ignored. Instead,
+        the number of clusters will be determined by when distances between
+        clusters first exceed the distance_threshold during agglomeration.
+        NOTE: You should set either ``n_clusters`` or ``distance_threshold``,
+        NOT both. If the ``distance_threshold`` is set then ``n_clusters`` is
+        ignored.
+
+    distances : ndarray, shape (n_nodes-1,) (optional)
+        Only used when distance_threshold is set.
+        distances[i] refers to the distance between children[i][0] and
+        children[i][1] when they are merged.
+
     Returns
     -------
     labels : array [n_samples]
         cluster labels for each point
 
     """
+    # cut the tree when distance threshold is reached
+    if distance_threshold:
+        n_clusters = np.count_nonzero(distances >= distance_threshold) + 1
     if n_clusters > n_leaves:
         raise ValueError('Cannot extract more clusters than samples: '
                          '%s clusters where given for a tree with %s leaves.'
@@ -604,6 +622,12 @@ class AgglomerativeClustering(BaseEstimator, ClusterMixin):
     n_clusters : int, default=2
         The number of clusters to find.
 
+    distance_threshold : float (optional)
+        The distance threshold to cluster at.
+        NOTE: You should set either ``n_clusters`` or ``distance_threshold``,
+        NOT both. If the ``distance_threshold`` is set then ``n_clusters`` is
+        ignored.
+
     affinity : string or callable, default: "euclidean"
         Metric used to compute the linkage. Can be "euclidean", "l1", "l2",
         "manhattan", "cosine", or 'precomputed'.
@@ -669,11 +693,13 @@ class AgglomerativeClustering(BaseEstimator, ClusterMixin):
 
     """
 
-    def __init__(self, n_clusters=2, affinity="euclidean",
+    def __init__(self, n_clusters=2, distance_threshold=None,
+                 affinity="euclidean",
                  memory=None,
                  connectivity=None, compute_full_tree='auto',
                  linkage='ward', pooling_func='deprecated'):
         self.n_clusters = n_clusters
+        self.distance_threshold = distance_threshold
         self.memory = memory
         self.connectivity = connectivity
         self.compute_full_tree = compute_full_tree
@@ -703,10 +729,16 @@ def fit(self, X, y=None):
         X = check_array(X, ensure_min_samples=2, estimator=self)
         memory = check_memory(self.memory)
 
-        if self.n_clusters <= 0:
+        if self.n_clusters is not None and self.n_clusters <= 0:
             raise ValueError("n_clusters should be an integer greater than 0."
                              " %s was provided." % str(self.n_clusters))
 
+        if self.n_clusters is None and self.distance_threshold is None:
+            raise ValueError("Either n_clusters (>0) or distance_threshold "
+                             "needs to be set, got n_clusters={} and "
+                             "distance_threshold={} instead.".format(
+                                self.n_clusters, self.distance_threshold))
+
         if self.linkage == "ward" and self.affinity != "euclidean":
             raise ValueError("%s was provided as affinity. Ward can only "
                              "work with euclidean distances." %
@@ -729,7 +761,7 @@ def fit(self, X, y=None):
         compute_full_tree = self.compute_full_tree
         if self.connectivity is None:
             compute_full_tree = True
-        if compute_full_tree == 'auto':
+        if compute_full_tree == 'auto' and self.distance_threshold is None:
             # Early stopping is likely to give a speed up only for
             # a large number of clusters. The actual threshold
             # implemented here is heuristic
@@ -743,14 +775,30 @@ def fit(self, X, y=None):
         if self.linkage != 'ward':
             kwargs['linkage'] = self.linkage
             kwargs['affinity'] = self.affinity
-        self.children_, self.n_components_, self.n_leaves_, parents = \
-            memory.cache(tree_builder)(X, connectivity,
-                                       n_clusters=n_clusters,
-                                       **kwargs)
+        distance_threshold = self.distance_threshold
+        # if distance_threshold is set then distances is returned
+        if distance_threshold is not None:
+            ch, n_comps, n_lvs, parents, distances = \
+                memory.cache(tree_builder)(X, connectivity,
+                                           n_clusters=n_clusters,
+                                           return_distance=True,
+                                           **kwargs)
+        else:
+            ch, n_comps, n_lvs, parents = \
+                memory.cache(tree_builder)(X, connectivity,
+                                           n_clusters=n_clusters,
+                                           **kwargs)
+        self.children_, self.n_components_, self.n_leaves_ = ch, n_comps, n_lvs
+
         # Cut the tree
         if compute_full_tree:
-            self.labels_ = _hc_cut(self.n_clusters, self.children_,
-                                   self.n_leaves_)
+            if distance_threshold is not None:
+                self.labels_ = _hc_cut(self.n_clusters, self.children_,
+                                       self.n_leaves_,
+                                       distance_threshold, distances)
+            else:
+                self.labels_ = _hc_cut(self.n_clusters, self.children_,
+                                       self.n_leaves_)
         else:
             labels = _hierarchical.hc_get_heads(parents, copy=False)
             # copy to avoid holding a reference on the original array
@@ -773,6 +821,12 @@ class FeatureAgglomeration(AgglomerativeClustering, AgglomerationTransform):
     n_clusters : int, default 2
         The number of clusters to find.
 
+    distance_threshold : float (optional)
+        The distance threshold to cluster at.
+        NOTE: You should set either ``n_clusters`` or ``distance_threshold``,
+        NOT both. If the ``distance_threshold`` is set then ``n_clusters`` is
+        ignored.
+
     affinity : string or callable, default "euclidean"
         Metric used to compute the linkage. Can be "euclidean", "l1", "l2",
         "manhattan", "cosine", or 'precomputed'.
@@ -835,7 +889,8 @@ class FeatureAgglomeration(AgglomerativeClustering, AgglomerationTransform):
         are merged to form node `n_features + i`
     """
 
-    def __init__(self, n_clusters=2, affinity="euclidean",
+    def __init__(self, n_clusters=2, distance_threshold=None,
+                 affinity="euclidean",
                  memory=None,
                  connectivity=None, compute_full_tree='auto',
                  linkage='ward', pooling_func=np.mean):
diff --git a/sklearn/cluster/tests/test_hierarchical.py b/sklearn/cluster/tests/test_hierarchical.py
index c4534663236b0..ffc992c85bf1c 100644
--- a/sklearn/cluster/tests/test_hierarchical.py
+++ b/sklearn/cluster/tests/test_hierarchical.py
@@ -13,6 +13,7 @@
 from scipy import sparse
 from scipy.cluster import hierarchy
 
+from sklearn.metrics.cluster.supervised import adjusted_rand_score
 from sklearn.utils.testing import assert_true
 from sklearn.utils.testing import assert_raises
 from sklearn.utils.testing import assert_equal
@@ -465,7 +466,8 @@ def test_connectivity_callable():
     connectivity = kneighbors_graph(X, 3, include_self=False)
     aglc1 = AgglomerativeClustering(connectivity=connectivity)
     aglc2 = AgglomerativeClustering(
-        connectivity=partial(kneighbors_graph, n_neighbors=3, include_self=False))
+        connectivity=partial(kneighbors_graph, n_neighbors=3,
+                             include_self=False))
     aglc1.fit(X)
     aglc2.fit(X)
     assert_array_equal(aglc1.labels_, aglc2.labels_)
@@ -527,11 +529,44 @@ def test_agg_n_clusters():
 
     rng = np.random.RandomState(0)
     X = rng.rand(20, 10)
-    for n_clus in [-1, 0]:
-        agc = AgglomerativeClustering(n_clusters=n_clus)
-        msg = ("n_clusters should be an integer greater than 0."
-               " %s was provided." % str(agc.n_clusters))
-        assert_raise_message(ValueError, msg, agc.fit, X)
+    n_clus = -1
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("n_clusters should be an integer greater than 0."
+           " %s was provided." % str(agc.n_clusters))
+    assert_raise_message(ValueError, msg, agc.fit, X)
+    n_clus = 0
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("n_clusters should be an integer greater than 0."
+           " %s was provided." % str(agc.n_clusters))
+    assert_raise_message(ValueError, msg, agc.fit, X)
+
+
+def test_agg_n_cluster_or_distance_threshold():
+    # Test that an error is raised when n_clusters
+    # and distance_threshold are None
+
+    n_clus = None
+    rng = np.random.RandomState(0)
+    X = rng.rand(20, 10)
+    # distance_threshold is None by default
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("Either n_clusters (>0) or distance_threshold "
+           "needs to be set, got n_clusters={} and "
+           "distance_threshold=None instead.".format(agc.n_clusters))
+    assert_raise_message(ValueError, msg, agc.fit, X)
+
+
+def test_agg_n_cluster_and_distance_threshold():
+    # Test that when distance_threshold is set n_clusters_ is unchanged
+
+    n_clus, dist_thresh = 10, 10
+    rng = np.random.RandomState(0)
+    X = rng.rand(20, 10)
+    agc = AgglomerativeClustering(n_clusters=n_clus,
+                                  distance_threshold=dist_thresh)
+    agc.fit(X)
+    # Expecting no errors here
+    assert_true(agc.n_clusters == n_clus)
 
 
 def test_affinity_passed_to_fix_connectivity():
@@ -559,3 +594,51 @@ def increment(self, *args, **kwargs):
     linkage_tree(X, connectivity=connectivity, affinity=fa.increment)
 
     assert_equal(fa.counter, 3)
+
+
+def test_agglomerative_clustering_with_distance_threshold():
+    # Check that we obtain the correct number of clusters with
+    # agglomerative clustering with distance_threshold.
+
+    rng = np.random.RandomState(0)
+    mask = np.ones([10, 10], dtype=np.bool)
+    n_samples = 100
+    X = rng.randn(n_samples, 50)
+    connectivity = grid_to_graph(*mask.shape)
+    # test when distance threshold is set to 10
+    distance_threshold = 10
+    for linkage in ("ward", "complete", "average"):
+        for conn in [None, connectivity]:
+            clustering = AgglomerativeClustering(
+                distance_threshold=distance_threshold,
+                connectivity=conn, linkage=linkage)
+            clustering.fit(X)
+            clusters_produced = clustering.labels_
+            num_clusters_produced = len(np.unique(clustering.labels_))
+            # test if the clusters produced match the point in the linkage tree
+            # where the distance exceeds the threshold
+            tree_builder = _TREE_BUILDERS[linkage]
+            children, n_components, n_leaves, parent, distances = \
+                tree_builder(X, connectivity=conn, n_clusters=None,
+                             return_distance=True)
+            num_clusters_at_threshold = np.count_nonzero(
+                distances >= distance_threshold) + 1
+            # test number of clusters produced
+            assert_true(num_clusters_at_threshold == num_clusters_produced)
+            # test clusters produced
+            clusters_at_threshold = _hc_cut(n_clusters=num_clusters_produced,
+                                            children=children,
+                                            n_leaves=n_leaves)
+            assert_true(np.array_equiv(clusters_produced,
+                                       clusters_at_threshold))
+
+
+def test_agglomerative_clustering_with_distance_threshold_edge_case():
+    # test boundary case of distance_threshold matching the distance
+    X = [[0], [1]]
+    for linkage in ("ward", "complete", "average"):
+        for threshold, y_true in [(0.5, [1, 0]), (1.0, [1, 0]), (1.5, [0, 0])]:
+            clusterer = AgglomerativeClustering(distance_threshold=threshold,
+                                                linkage=linkage)
+            y_pred = clusterer.fit_predict(X)
+            assert_equal(1, adjusted_rand_score(y_true, y_pred))


Code review comments:
diff_hunk: @@ -527,11 +529,44 @@ def test_agg_n_clusters():
 
     rng = np.random.RandomState(0)
     X = rng.rand(20, 10)
-    for n_clus in [-1, 0]:
-        agc = AgglomerativeClustering(n_clusters=n_clus)
-        msg = ("n_clusters should be an integer greater than 0."
-               " %s was provided." % str(agc.n_clusters))
-        assert_raise_message(ValueError, msg, agc.fit, X)
+    n_clus = -1
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("n_clusters should be an integer greater than 0."
+           " %s was provided." % str(agc.n_clusters))
+    assert_raise_message(ValueError, msg, agc.fit, X)
+    n_clus = 0
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("n_clusters should be an integer greater than 0."
+           " %s was provided." % str(agc.n_clusters))
+    assert_raise_message(ValueError, msg, agc.fit, X)
+
+
+def test_agg_n_cluster_or_distance_threshold():
+    # Test that an error is raised when n_clusters
+    # and distance_threshold are None
+
+    n_clus = None
+    rng = np.random.RandomState(0)
+    X = rng.rand(20, 10)
+    # distance_threshold is None by default
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("Either n_clusters (>0) or distance_threshold "
+           "needs to be set, got n_clusters={} and "
+           "distance_threshold=None instead.".format(agc.n_clusters))
+    assert_raise_message(ValueError, msg, agc.fit, X)
+
+
+def test_agg_n_cluster_and_distance_threshold():
+    # Test that when distance_threshold is set n_clusters_ is unchanged
path: sklearn/cluster/tests/test_hierarchical.py
text: @user1:
This isn't the right comment anymore

diff_hunk: @@ -743,14 +775,30 @@ def fit(self, X, y=None):
         if self.linkage != 'ward':
             kwargs['linkage'] = self.linkage
             kwargs['affinity'] = self.affinity
-        self.children_, self.n_components_, self.n_leaves_, parents = \
-            memory.cache(tree_builder)(X, connectivity,
-                                       n_clusters=n_clusters,
-                                       **kwargs)
+        distance_threshold = self.distance_threshold
+        # if distance_threshold is set then distances is returned
+        if distance_threshold is not None:
+            ch, n_comps, n_lvs, parents, distances = \
+                memory.cache(tree_builder)(X, connectivity,
+                                           n_clusters=n_clusters,
+                                           return_distance=True,
+                                           **kwargs)
+        else:
+            ch, n_comps, n_lvs, parents = \
+                memory.cache(tree_builder)(X, connectivity,
+                                           n_clusters=n_clusters,
+                                           **kwargs)
+        self.children_, self.n_components_, self.n_leaves_ = ch, n_comps, n_lvs
+
         # Cut the tree
         if compute_full_tree:
-            self.labels_ = _hc_cut(self.n_clusters, self.children_,
-                                   self.n_leaves_)
+            if distance_threshold is not None:
+                self.labels_ = _hc_cut(self.n_clusters, self.children_,
path: sklearn/cluster/hierarchical.py
text: @user1:
why don't you just calculate n_clusters from distances here and avoid modifying `_hc_cut` (adding two parameters for the sake of 1 line of logic)?

diff_hunk: @@ -527,11 +529,44 @@ def test_agg_n_clusters():
 
     rng = np.random.RandomState(0)
     X = rng.rand(20, 10)
-    for n_clus in [-1, 0]:
-        agc = AgglomerativeClustering(n_clusters=n_clus)
-        msg = ("n_clusters should be an integer greater than 0."
-               " %s was provided." % str(agc.n_clusters))
-        assert_raise_message(ValueError, msg, agc.fit, X)
+    n_clus = -1
+    agc = AgglomerativeClustering(n_clusters=n_clus)
+    msg = ("n_clusters should be an integer greater than 0."
+           " %s was provided." % str(agc.n_clusters))
path: sklearn/cluster/tests/test_hierarchical.py
text: @user1:
Either use a loop to test both -1 and 0, or put -1 directly in the expected error message