[MRG+2] OPTICS: add extract_xi method

doc/modules/classes.rst

@@ -114,6 +114,7 @@ Functions
 
    cluster.affinity_propagation
    cluster.cluster_optics_dbscan
+   cluster.cluster_optics_xi
    cluster.compute_optics_graph
    cluster.dbscan
    cluster.estimate_bandwidth

doc/modules/clustering.rst

@@ -91,6 +91,12 @@ Overview of clustering methods
      - Non-flat geometry, uneven cluster sizes
      - Distances between nearest points
 
+   * - :ref:`OPTICS <optics>`
+     - minimum cluster membership
+     - Very large ``n_samples``, large ``n_clusters``
+     - Non-flat geometry, uneven cluster sizes, variable cluster density
+     - Distances between points
+
    * - :ref:`Gaussian mixtures <mixture>`
      - many
      - Not scalable
@@ -806,6 +812,11 @@ by black points below.
     be used (e.g., with sparse matrices). This matrix will consume n^2 floats.
     A couple of mechanisms for getting around this are:
 
+    - Use :ref:`OPTICS <optics>` clustering in conjunction with the
+      `extract_dbscan` method. OPTICS clustering also calculates the full
+      pairwise matrix, but only keeps one row in memory at a time (memory
+      complexity n).
+
     - A sparse radius neighborhood graph (where missing entries are presumed to
       be out of eps) can be precomputed in a memory-efficient way and dbscan
       can be run over this with ``metric='precomputed'``.  See
@@ -828,6 +839,92 @@ by black points below.
    Schubert, E., Sander, J., Ester, M., Kriegel, H. P., & Xu, X. (2017).
    In ACM Transactions on Database Systems (TODS), 42(3), 19.
 
+.. _optics:
+
+OPTICS
+======
+
+The :class:`OPTICS` algorithm shares many similarities with the :class:`DBSCAN`
+algorithm, and can be considered a generalization of DBSCAN that relaxes the
+``eps`` requirement from a single value to a value range. The key difference
+between DBSCAN and OPTICS is that the OPTICS algorithm builds a *reachability*
+graph, which assigns each sample both a ``reachability_`` distance, and a spot
+within the cluster ``ordering_`` attribute; these two attributes are assigned
+when the model is fitted, and are used to determine cluster membership. If
+OPTICS is run with the default value of *inf* set for ``max_eps``, then DBSCAN
+style cluster extraction can be performed repeatedly in linear time for any
+given ``eps`` value using the ``cluster_optics_dbscan`` method. Setting
+``max_eps`` to a lower value will result in shorter run times, and can be
+thought of as the maximum neighborhood radius from each point to find other
+potential reachable points.
+
+.. |optics_results| image:: ../auto_examples/cluster/images/sphx_glr_plot_optics_001.png
+        :target: ../auto_examples/cluster/plot_optics.html
+        :scale: 50
+
+.. centered:: |optics_results|
+
+The *reachability* distances generated by OPTICS allow for variable density
+extraction of clusters within a single data set. As shown in the above plot,
+combining *reachability* distances and data set ``ordering_`` produces a
+*reachability plot*, where point density is represented on the Y-axis, and
+points are ordered such that nearby points are adjacent. 'Cutting' the
+reachability plot at a single value produces DBSCAN like results; all points
+above the 'cut' are classified as noise, and each time that there is a break
+when reading from left to right signifies a new cluster. The default cluster
+extraction with OPTICS looks at the steep slopes within the graph to find
+clusters, and the user can define what counts as a steep slope using the
+parameter ``xi``. There are also other possibilities for analysis on the graph
+itself, such as generating hierarchical representations of the data through
+reachability-plot dendrograms, and the hierarchy of clusters detected by the
+algorithm can be accessed through the ``cluster_hierarchy_`` parameter. The
+plot above has been color-coded so that cluster colors in planar space match
+the linear segment clusters of the reachability plot. Note that the blue and
+red clusters are adjacent in the reachability plot, and can be hierarchically
+represented as children of a larger parent cluster.
+
+.. topic:: Examples:
+
+     * :ref:`sphx_glr_auto_examples_cluster_plot_optics.py`
+
+
+.. topic:: Comparison with DBSCAN
+
+    The results from OPTICS ``cluster_optics_dbscan`` method and DBSCAN are
+    very similar, but not always identical; specifically, labeling of periphery
+    and noise points. This is in part because the first samples of each dense
+    area processed by OPTICS have a large reachability value while being close
+    to other points in their area, and will thus sometimes be marked as noise
+    rather than periphery. This affects adjacent points when they are
+    considered as candidates for being marked as either periphery or noise.
+
+    Note that for any single value of ``eps``, DBSCAN will tend to have a
+    shorter run time than OPTICS; however, for repeated runs at varying ``eps``
+    values, a single run of OPTICS may require less cumulative runtime than
+    DBSCAN. It is also important to note that OPTICS' output is close to
+    DBSCAN's only if ``eps`` and ``max_eps`` are close.
+
+.. topic:: Computational Complexity
+
+    Spatial indexing trees are used to avoid calculating the full distance
+    matrix, and allow for efficient memory usage on large sets of samples.
+    Different distance metrics can be supplied via the ``metric`` keyword.
+
+    For large datasets, similar (but not identical) results can be obtained via
+    `HDBSCAN <https://hdbscan.readthedocs.io>`_. The HDBSCAN implementation is
+    multithreaded, and has better algorithmic runtime complexity than OPTICS,
+    at the cost of worse memory scaling. For extremely large datasets that
+    exhaust system memory using HDBSCAN, OPTICS will maintain *n* (as opposed
+    to *n^2*) memory scaling; however, tuning of the ``max_eps`` parameter
+    will likely need to be used to give a solution in a reasonable amount of
+    wall time.
+
+.. topic:: References:
+
+ *  "OPTICS: ordering points to identify the clustering structure."
+    Ankerst, Mihael, Markus M. Breunig, Hans-Peter Kriegel, and Jörg Sander.
+    In ACM Sigmod Record, vol. 28, no. 2, pp. 49-60. ACM, 1999.
+
 .. _birch:
 
 Birch

doc/whats_new/v0.21.rst

@@ -83,7 +83,8 @@ Support for Python 3.4 and below has been officially dropped.
 - |MajorFeature| A new clustering algorithm: :class:`cluster.OPTICS`: an
   algoritm related to :class:`cluster.DBSCAN`, that has hyperparameters easier
   to set and that scales better, by :user:`Shane <espg>`,
-  :user:`Adrin Jalali <adrinjalali>`, and :user:`Erich Schubert <kno10>`.
+  `Adrin Jalali`_, :user:`Erich Schubert <kno10>`, `Hanmin Qin`_, and
+  :user:`Assia Benbihi <assiaben>`.
 
 - |Fix| Fixed a bug where :class:`cluster.Birch` could occasionally raise an
   AttributeError. :issue:`13651` by `Joel Nothman`_.

examples/cluster/plot_cluster_comparison.py

@@ -74,14 +74,20 @@
                 'damping': .9,
                 'preference': -200,
                 'n_neighbors': 10,
-                'n_clusters': 3}
+                'n_clusters': 3,
+                'min_samples': 20,
+                'xi': 0.05,
+                'min_cluster_size': 0.1}
 
 datasets = [
     (noisy_circles, {'damping': .77, 'preference': -240,
-                     'quantile': .2, 'n_clusters': 2}),
+                     'quantile': .2, 'n_clusters': 2,
+                     'min_samples': 20, 'xi': 0.25}),
     (noisy_moons, {'damping': .75, 'preference': -220, 'n_clusters': 2}),
-    (varied, {'eps': .18, 'n_neighbors': 2}),
-    (aniso, {'eps': .15, 'n_neighbors': 2}),
+    (varied, {'eps': .18, 'n_neighbors': 2,
+              'min_samples': 5, 'xi': 0.035, 'min_cluster_size': .2}),
+    (aniso, {'eps': .15, 'n_neighbors': 2,
+             'min_samples': 20, 'xi': 0.1, 'min_cluster_size': .2}),
     (blobs, {}),
     (no_structure, {})]
 
@@ -116,6 +122,9 @@
         n_clusters=params['n_clusters'], eigen_solver='arpack',
         affinity="nearest_neighbors")
     dbscan = cluster.DBSCAN(eps=params['eps'])
+    optics = cluster.OPTICS(min_samples=params['min_samples'],
+                            xi=params['xi'],
+                            min_cluster_size=params['min_cluster_size'])
     affinity_propagation = cluster.AffinityPropagation(
         damping=params['damping'], preference=params['preference'])
     average_linkage = cluster.AgglomerativeClustering(
@@ -133,6 +142,7 @@
         ('Ward', ward),
         ('AgglomerativeClustering', average_linkage),
         ('DBSCAN', dbscan),
+        ('OPTICS', optics),
         ('Birch', birch),
         ('GaussianMixture', gmm)
     )

examples/cluster/plot_optics.py

@@ -0,0 +1,98 @@
+"""
+===================================
+Demo of OPTICS clustering algorithm
+===================================
+Finds core samples of high density and expands clusters from them.
+This example uses data that is generated so that the clusters have
+different densities.
+The :class:`sklearn.cluster.OPTICS` is first used with its Xi cluster detection
+method, and then setting specific thresholds on the reachability, which
+corresponds to :class:`sklearn.cluster.DBSCAN`. We can see that the different
+clusters of OPTICS's Xi method can be recovered with different choices of
+thresholds in DBSCAN.
+"""
+
+# Authors: Shane Grigsby <refuge@rocktalus.com>
+#          Adrin Jalali <adrin.jalali@gmail.com>
+# License: BSD 3 clause
+
+
+from sklearn.cluster import OPTICS, cluster_optics_dbscan
+import matplotlib.gridspec as gridspec
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Generate sample data
+
+np.random.seed(0)
+n_points_per_cluster = 250
+
+C1 = [-5, -2] + .8 * np.random.randn(n_points_per_cluster, 2)
+C2 = [4, -1] + .1 * np.random.randn(n_points_per_cluster, 2)
+C3 = [1, -2] + .2 * np.random.randn(n_points_per_cluster, 2)
+C4 = [-2, 3] + .3 * np.random.randn(n_points_per_cluster, 2)
+C5 = [3, -2] + 1.6 * np.random.randn(n_points_per_cluster, 2)
+C6 = [5, 6] + 2 * np.random.randn(n_points_per_cluster, 2)
+X = np.vstack((C1, C2, C3, C4, C5, C6))
+
+clust = OPTICS(min_samples=50, xi=.05, min_cluster_size=.05)
+
+# Run the fit
+clust.fit(X)
+
+labels_050 = cluster_optics_dbscan(reachability=clust.reachability_,
+                                   core_distances=clust.core_distances_,
+                                   ordering=clust.ordering_, eps=0.5)
+labels_200 = cluster_optics_dbscan(reachability=clust.reachability_,
+                                   core_distances=clust.core_distances_,
+                                   ordering=clust.ordering_, eps=2)
+
+space = np.arange(len(X))
+reachability = clust.reachability_[clust.ordering_]
+labels = clust.labels_[clust.ordering_]
+
+plt.figure(figsize=(10, 7))
+G = gridspec.GridSpec(2, 3)
+ax1 = plt.subplot(G[0, :])
+ax2 = plt.subplot(G[1, 0])
+ax3 = plt.subplot(G[1, 1])
+ax4 = plt.subplot(G[1, 2])
+
+# Reachability plot
+colors = ['g.', 'r.', 'b.', 'y.', 'c.']
+for klass, color in zip(range(0, 5), colors):
+    Xk = space[labels == klass]
+    Rk = reachability[labels == klass]
+    ax1.plot(Xk, Rk, color, alpha=0.3)
+ax1.plot(space[labels == -1], reachability[labels == -1], 'k.', alpha=0.3)
+ax1.plot(space, np.full_like(space, 2., dtype=float), 'k-', alpha=0.5)
+ax1.plot(space, np.full_like(space, 0.5, dtype=float), 'k-.', alpha=0.5)
+ax1.set_ylabel('Reachability (epsilon distance)')
+ax1.set_title('Reachability Plot')
+
+# OPTICS
+colors = ['g.', 'r.', 'b.', 'y.', 'c.']
+for klass, color in zip(range(0, 5), colors):
+    Xk = X[clust.labels_ == klass]
+    ax2.plot(Xk[:, 0], Xk[:, 1], color, alpha=0.3)
+ax2.plot(X[clust.labels_ == -1, 0], X[clust.labels_ == -1, 1], 'k+', alpha=0.1)
+ax2.set_title('Automatic Clustering\nOPTICS')
+
+# DBSCAN at 0.5
+colors = ['g', 'greenyellow', 'olive', 'r', 'b', 'c']
+for klass, color in zip(range(0, 6), colors):
+    Xk = X[labels_050 == klass]
+    ax3.plot(Xk[:, 0], Xk[:, 1], color, alpha=0.3, marker='.')
+ax3.plot(X[labels_050 == -1, 0], X[labels_050 == -1, 1], 'k+', alpha=0.1)
+ax3.set_title('Clustering at 0.5 epsilon cut\nDBSCAN')
+
+# DBSCAN at 2.
+colors = ['g.', 'm.', 'y.', 'c.']
+for klass, color in zip(range(0, 4), colors):
+    Xk = X[labels_200 == klass]
+    ax4.plot(Xk[:, 0], Xk[:, 1], color, alpha=0.3)
+ax4.plot(X[labels_200 == -1, 0], X[labels_200 == -1, 1], 'k+', alpha=0.1)
+ax4.set_title('Clustering at 2.0 epsilon cut\nDBSCAN')
+
+plt.tight_layout()
+plt.show()

sklearn/cluster/__init__.py

@@ -11,7 +11,8 @@
                            FeatureAgglomeration)
 from .k_means_ import k_means, KMeans, MiniBatchKMeans
 from .dbscan_ import dbscan, DBSCAN
-from .optics_ import OPTICS, cluster_optics_dbscan, compute_optics_graph
+from .optics_ import (OPTICS, cluster_optics_dbscan, compute_optics_graph,
+                      cluster_optics_xi)
 from .bicluster import SpectralBiclustering, SpectralCoclustering
 from .birch import Birch
 
@@ -21,6 +22,7 @@
            'DBSCAN',
            'OPTICS',
            'cluster_optics_dbscan',
+           'cluster_optics_xi',
            'compute_optics_graph',
            'KMeans',
            'FeatureAgglomeration',